Primes Race (Mathematica Efficiency)

Question

I am currently working with a paper that deals with this concept of Prime Races. You essentially create a large list of prime numbers and then split that list into two teams. You are assigned to team 1 if your remainder is 1 after being divided by 4. You are assigned to team 2 if your remainder is 3. You then choose an x value and see how many primes fall under that x value for each team. It looks as though Team 1 always wins however there are certain cases this is not true. For instance at x=26,861, Team 1= 1472 while Team 2 = 1473. I'm trying to expand x out to 1 million and find where this occurs at other places in a computational efficient way. However the method I'm using takes forever. Are there any other ideas or concepts that may work?

Here is my Mathematica code:

M = 100000;
listPrimes = Table[Prime[i], {i, 1, M}];
list1 = Select[listPrimes, Mod[#1, 4] == 3 &];
list2 = Select[listPrimes, Mod[#1, 4] == 1 &];
x = 26862;
list1cnt = Length[Select[list1, # <= x &]]
list2cnt = Length[Select[list2, # <= x &]]

mchk = 30000;
For[x = 26850, x <= mchk, x++,
  list1cnt = Length[Select[list1, # <= x &]];
  list2cnt = Length[Select[list2, # <= x &]];
  If[list2cnt > list1cnt, Print["Found One:", x]]
  ];

I'm testing it between 26850 and 30000 and this takes a couple minutes so 1 to 1 million would take far too long.

Thanks for any advice or suggestions!

Answer 1

Improved Solution

The challenge posed in the question is determine for any x up to M = 1000000 whether list1 or list2 has more elements less than x. The question suggests that M primes are needed to make that determination. In fact, only the substantially smaller (by a factor of about 0.08) set of primes less than or equal to M is needed. It is computed by

M = 1000000;
listPrimes = Table[Prime[i], {i, 1, PrimePi[M]}];
list1 = Select[listPrimes, Mod[#1, 4] == 3 &];
list2 = Select[listPrimes, Mod[#1, 4] == 1 &];

Next, we count the number of primes less than any integer from 1 to M.

listm = Append[Prepend[list1, 1], M + 1];
count1 = Flatten@Table[ConstantArray[i - 1, listm[[i + 1]] - listm[[i]]], 
    {i, Length[listm] - 1}];
listm = Append[Prepend[list2, 1], M + 1];
count2 = Flatten@Table[ConstantArray[i - 1, listm[[i + 1]] - listm[[i]]], 
    {i, Length[listm] - 1}];

and Tally the number of instances that one or the other list is ahead.

Tally@MapThread[Sign[#1 - #2] &, {count1, count2}]
(* {{0, 1352}, {1, 995242}, {-1, 3406}} *)

list1 is ahead for 995242 values of integer x between 1 and M, list2 ahead for 3406 values of x, and the lists are tied for 1352 values of x. The first ten x, for instance, at which list2 is ahead is given by

Flatten@Position[MapThread[Sign[#1 - #2] &, {count1, count2}], -1][[1 ;; 10]]
(* {26861, 26862, 616841, 616842, 616849, 616850, 616851, 616852, 616853, 616854} *)

The AbsoluteTiming for this entire calculation is about 1.9 sec on my PC.

Solution for M = 100000000

My 8 GB PC can just barely handle M = 100000000, requiring about 200 sec and over 8.2 GB (some on disk). Thus, run time varies linearly with M, as expected. (Larger M could be handled by breaking the calculation into parts to reduce memory usage.) The Tally results are

(* {{0, 3866}, {1, 99965510}, {-1, 30624}} *)

The list2 "wins" are clustered into three groups, much as shown in the earlier plot below.

Earlier Solution

This earlier solution was derived for 1000000 primes. It compares the magnitude of the nth prime of list1 with that of list2. Whichever list has the smaller prime is ahead at that point. In effect, it involves sampling but nonetheless gives a good qualitative picture of the behavior of the race.

Count[MapThread[#1 < #2 &, {list1[[1 ;; Length[list2]]], list2}], False]
(* 1034 *)

Thus, there are 1034 instances among the first million primes in which a list1 element is larger than the corresponding list2 element. The first one is readily found with

Position[MapThread[#1 < #2 &, {list1[[1 ;; Length[list2]]], list2}][[1 ;; 10000]], False]
(* {{1473}} *)

{list1[[%[[1, 1]]]], list2[[%[[1, 1]]]]}
(* {26863, 26861} *)

The distribution of rare cases can be plotted as follows.

ListPlot[{#, list1[[#]]} & /@ Flatten@Position[
    MapThread[#1 < #2 &, {list1[[1 ;; Length[list2]]], list2}], False], 
    PlotRange -> {{1, Length[list2]}, {1, list2[[-1]]}}]

So, the rare cases where a list1 element is larger than the corresponding list2 element are clustered into just three groups, the first barely visible near the origin. Closer examination reveals that there is 1 case in the first cluster, about 150 cases in the second, and about 1330 in the third.

Answer 2

The solution is fairly simple,

Prime[Sign@Accumulate[Mod[Prime@Range@PrimePi[1*^6],4]-2]~Position~-1]

where

• Prime@Range@PrimePi[1*^6] gives the list of all primes up to 1*^6;
• Mod[...,4]-2 automatically threads over the list above, giving 1 or -1 for two teams respectively. Note that Mod[2,4]-2 equals to 0.
• Sign@Accumulate[...] gives a list of results about which team is leading so far.
• Prime[...~Position~-1] gives all primes team 2 is leading when up to which.

---

It takes 0.056s on my laptop to get the final result,

{{26861},{616841},{616849},...,{633653},{633667},{633797}}

---

Other techniques

Here are some techniques to improve performance, but be less elegant.

• Mod[#,4]& --> BitAnd[#,3]&
• Position[#,-1]& --> Pick[Range@Length@#,#,-1]; they are expressing the same idea in this case.