Primes Race (Mathematica Efficiency)

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!

• Just do addition each step, no need to reselect (which makes an O(1) computation into an O(n) computation). – Daniel Lichtblau Mar 1 '16 at 23:30
• @Dops "It looks as though Team 1 always wins however there are certain cases this is not true." The observed effect of team 1 winning is called the Chebychev bias. Hardy and Littlewood have shown, however, that the "lead" changes infinitely often. More information here: Weisstein, Eric W. "Modular Prime Counting Function." From MathWorld--A Wolfram Web Resource. mathworld.wolfram.com/ModularPrimeCountingFunction.html – Dr. Wolfgang Hintze Mar 9 '16 at 10:43

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.

• Hm very interesting. I've never used MapThread and therefore don't really know what's going on however I'll try and pick it apart. Thank you for the feedback. – Dops Mar 1 '16 at 23:52
• @Dops MapThread is like Map but applies the function in its first argument to corresponding elements of multiple lists. – bbgodfrey Mar 1 '16 at 23:55
• Could someone try and explain the syntax of what exactly is going? Also do you mean "Thus, there are 1034 instances among the first million primes in which a list1 element is smaller than the corresponding list2 element." I'm looking for when list 1 has less values than list 2. – Dops Mar 2 '16 at 3:36
• @Dops Does the additional material in my answer adequately address your comment just above? My original answer looked at x equal to the primes in list1, whereas the more recent answer looks at all x up to M. Although the qualitative behavior is similar, my newer answer is sufficiently fast that it ought to be the one used. – bbgodfrey Mar 2 '16 at 20:45

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.
• This is a very clever approach (+1), yielding a solution equivalent to my "Earlier Solution" but somewhat faster. Do you see a similar approach that could reproduce my "Improved Solution" faster? – bbgodfrey Mar 2 '16 at 16:09
• @bbgodfrey I do not think it will be faster since checking for every integers is actually much less efficient. – njpipeorgan Mar 3 '16 at 0:40