Fastest way to generate prime numbers

Questions
,
1

Has anyone invesitgated what is the fastest way to generate prime numbers sequentially in Mathematica, and whether this is effectively parallelizable?

I’m trying to do some statistics on prime numbers and right now I am generating them in batches of size ‘chunk’, something like this:

(* chunk is a large number *)
Prime@Range[r, r + chunk]
(* do stuff *)
r += chunk

The performance characteristics of the prime-related functions are rather opaque and they seem to rely in some sort of internal caching (for example for as long as I’m generating all primes, using Prime[] instead of NextPrime[] didn’t seem to carry any great penalties).

I was wondering if someone has already investigated the fastest way to generate primes in chunks, sequentially. Because of the caching that’s going on it’s not clear how to Parallelize[] this effectively, without duplicating a lot of the work among subkernels.

2

Pointing out a somewhat related discussion - perhaps there are some tidbits there:

How to find the first 100 prime twins without a do loop or iteration

I wonder if MSE has some related discussions - but have not found any yet.

3

Depends on the size. For 10 digits I would expect your method to be fine, although using NextPrime is faster. Each of the below were done from a clean start so as to remove cache effects.

In[1]:= Timing[p1 = NextPrime[252097800623, Range[100]];]

(* Out[1]= {0.004000, Null} *)

In[1]:= Timing[p2 = Prime[Range[10^10 + 1, 10^10 + 100]];]

(* Out[1]= {0.752000, Null} *)

Double the size of the starting point and you won’t have a choice of methods anymore.

4

Is it feasible to parallelize listing primes sequentially using Mathematica?

Generally, Prime won’t compute all primes up to n. Prime[10^6] is much faster than Prime@Range[10^6]. If this weren’t true, then parallelization wouldn’t be feasible as all the work would be duplicated among subkernels.

This does indeed work:

AbsoluteTiming[
 Join @@ ParallelTable[Prime[Range[10^6] + k 10^6], {k, 0, 3}];
]

And it does indeed speed things up considrably compared to the non-parallel version: Prime@Range[4*10^6] (I got a speedup of about 2.5 for this example using 4 subkernels).

The above also applies to NextPrime when used as NextPrime[start, Range[…]]. Unfortunately this naive approach to parallelization breaks down (stops being efficient) if used for very large primes, say starting the calculation from Prime[10^11] == 2760727302517. This is in fact the magnitude I’m interested in. I computed primes up to just above 10^11, but it wasn’t a serious calculation (I just didn’t have a reason to stop it, as the server was free, so I let it run for a long time and watched what happens using the magic of Dynamic …). Now I’m interested in repeating the calculation again up to 10^10 or 10^11, and I’m hoping to be able to reduce the calculation time by parallelization.

I expect that Prime and NextPrime do some smart things internally and automatically choose the best method depending on the magnitude of the numbers. This makes it difficult to guess how to parallilize best, e.g. what size of chunks to use instead of the 10^6 above, or even if the chunk size should be constant. Prime is a black box to the user.

If I understand Daniel’s comment correctly, then NexrPrime switches methods above 10^10 or so, so perhaps for very large numbers it will not be possible to get any speedup with this type of parallelization.

5

I do not recall whether or where NextPrime switches methods.

As for parallelization, if you are in a range where Prime works then you might try a ParallelTable. Below is more or less what you did and it seems to give a notable speed boost in the range in question.

p1 = Prime[10^10];
n = 10^5;
first = Prime[Range[10^10, 10^10 + 4 n - 1, n]];
AbsoluteTiming[
 t1 = Flatten[
    ParallelTable[NextPrime[j, Range[0, n - 1]], {j, first}]];]

Compare to:

AbsoluteTiming[t2 = NextPrime[p1, Range[0, 4*(n + 1) - 1]];]

There may be better ways to parallelize this as far as speed goes.

6

Actually one can do quite well in regions above the working zone of Prime. Say we want to get the first n primes larger than m, and we intend to parallelize over nproc processors. We can get reasonable approximations of where the n/nproc groupings begin, overshoot each by some fraction, and take Union at the end. Also in case we underdo the total we may need a final run to make up the deficit. The code below does what I have described here.

nPrimesAboveM[m_, n_, nproc_] := Module[
  {markers, res, primes, frac = Ceiling[105/100*n/nproc] + 5},
  markers =
   Round[NestList[# + NIntegrate[Log[x], {x, #, # + n/nproc}] &, m,
     nproc - 1]];
  primes =
   ParallelTable[NextPrime[markers[[j]], Range[frac]], {j, nproc}];
  res = Partition[primes, 2, 1];
  Map[If[Last[#[[1]]] < First[#[[2]]], Print["sinister gap"]] &,
   res];
  res = Union @@ primes;
  If[Length[res] < n,
   res = Join[res, NextPrime[Last[res], Range[n - Length[res]]]]];
  Take[res, n]
  ]

Here is an example, finding a million primes after 10^20.

AbsoluteTiming[plist1 = nPrimesAboveM[10^20, 10^6, 4];]

Out[47]= {111.080568, Null}

It is around 3x faster than its nonparallelized counterpart.

AbsoluteTiming[plist2 = NextPrime[10^20, Range[10^6]];]

Out[10]= {366.245392, Null}

plist2 === plist1 (* check result *)

Out[48]= True

One important caveat: it is not flawless. I have not put in a patch for the case where we get a gap, and this case can arise.