Why does iterating Prime in reverse order require much more time?

Question

Say I would like to display the $10$ greatest primes that are less than $10^5$. I could do the following:

AbsoluteTiming[
 M = 10^5; m = PrimePi[M];
 prms = Prime[#] & /@ Range[1, m];
 prms[[#]] & /@ Range[-1, -10, -1]
 ]

And the result comes out :

{0.0156250, {99991, 99989, 99971, 99961, 99929, 99923, 99907, 99901, 99881, 99877}}

But if I tried to do in in reverse,

AbsoluteTiming[
 M = 10^5; m = PrimePi[M];
 prms = Prime[#] & /@ Range[m, 1, -1];
 prms[[#]] & /@ Range[1, 10]
 ]

the process takes a whole lot longer:

{0.6250000, {99991, 99989, 99971, 99961, 99929, 99923, 99907, 99901, 
  99881, 99877}}

Using the second method, I can't even increase M to $10^6$, as the program takes extremely long to execute. Can anybody offer some insight into this ? $\;$ Am I essentially not doing the same thing in both cases ?

Answer 1

Given a large n, to find k largest primes below n (as well as above) the best approach uses NextPrime (it has been added to Mathematica 6) :

NextPrime[n] gives the next prime above n.

NextPrime[n,k] gives the k-th prime above n. If k is negative it gives k-th largest prime below n.

k need not be a single number but it may be a list of integers, so if we are looking for k consecutive primes we can take advanted of Range, e.g. :

NextPrime[ 100000, Range[-10, -1]]

{99877, 99881, 99901, 99907, 99923, 99929, 99961, 99971, 99989, 99991}

The issue with Prime and PrimePi is that they are internally related however their documentation pages are not very informative. There are certain limitations of these functions (look at a related question : What is so special about Prime? ). Prime calls PrimePi (e.g. this comment by Oleksandr R.) if Prime[n] < 25 10^13. One can guess what is going on from Some Notes on Internal Implementation where it says:

Prime and PrimePi use sparse caching and sieving. For large $n$, the Lagarias-Miller-Odlyzko algorithm for PrimePi is used, based on asymptotic estimates of the density of primes, and is inverted to give Prime.

So if one has found a large prime, generically the system definitely has found some close primes too (sparse caching and sieving) and of course internal algorithms are not symmetric around a large $n$, i.e. finding closest $k$ primes below and above $n$ is not symmetric (basically it is implied by decreasing density of primes (globally) but directly it is determined by the Lagarias-Miller-Odlyzko method ). For more information take a look at this crucial reference : Computing $ \pi(x)$: the Meissel-Lehmer method. If you want to find really large primes a fast algorithm should use PrimeQ however it is known to be correct only for $n < 10^{16}$. Another algorithm which is correct for all natural n is much slower, one can find it in PrimalityProving package .

Answer 2

Addendum

If you just want the greatest 10 primes less than M, you can start from Prime[PrimePi[M]-9]. By doing so, you gain a speed increase of 2 orders of magnitude when M = 100000.

M = 100000;
m = PrimePi[M]
AbsoluteTiming[Table[Prime[k], {k, m - 9, m}]]

9592
{0.000171, {99877, 99881, 99901, 99907, 99923, 99929, 99961, 99971, 99989, 99991}}

Now backwards (after closing and re-opening Mathematica and defining M and m.). This is 10 x slower than the above, but still faster than starting from Prime[1].

AbsoluteTiming[Table[Prime[k], {k, m, m - 9, -1}]]

{0.001467, {99991, 99989, 99971, 99961, 99929, 99923, 99907, 99901, 99881, 99877}}

One might expect that NextPrime would provide an advantage. We restart Mathematica and find out that it is only marginally better than our last attempt.

AbsoluteTiming[NestList[NextPrime, Prime[m - 9], 9]]

{0.001193, {99877, 99881, 99901, 99907, 99923, 99929, 99961, 99971, 99989, 99991}}

---

My Earlier answers about the order of the iterator.

I conducted the following tests starting at Prime[1]. Even though this is a very inefficient approach, the results will be easier to compare to those of the OP.

We want to rule out Range so we take it out of the timing. Likewise, we set M and m first, before conducting the timing.

M = 100000;
m = PrimePi[M];
r1 = Range[1, m];
r2 = Range[m, 1, -1];

Backwards request of primes much slower than forwards request (the original observation)

r1 == Reverse[r2]
AbsoluteTiming[prms = Prime[#] & /@ r1;]
AbsoluteTiming[prms = Prime[#] & /@ r2;]

True
{0.009920, Null}
{0.354147, Null}

Puzzling. Yet absolutely consistent with what you observed.

It appears that r1 may be stored in a form that can be processed more quickly than r2.

---

Half backwards, Half forwards: Only half as bad as all backwards.

Let's look a little deeper.

r3 is the sub-list of odd numbers, followed by the even numbers.

r5 is like r3 but the odd numbers were reversed.

r3 = Transpose[Partition[r1, 2]] // Flatten;
r4 = Transpose[Partition[r1, 2]];
r5 = Join[Reverse[r4[[1]]], r4[[2]]];

Let's time the processing of r3 and r5:

AbsoluteTiming[prms = Prime[#] & /@ r3;]
AbsoluteTiming[prms = Prime[#] & /@ r5;]

{0.012765, Null}
{0.205590, Null}

Half the numbers in r5 were ordered from larger to smaller. The delay was about 1/2 the delay we saw when all numbers were ordered from larger to smaller.

Provisional conclusion: Mathematica "counts" primes forwards faster than backwards(?)

---

Random Ordering

r6 = RandomSample[r1, 9592];
AbsoluteTiming[prms = Prime[#] & /@ r6;]

{0.190770, Null}

The (pseudo)random ordering of numbers runs as fast as the backwards-forwards list!

---

Lingering Question: Why should a random ordering of numbers be faster to process than a "backwards request" for primes? I would have guessed that the random case would be at least as slow, or slower than the backward case.

P.S. @Szabolcs Please look at the results above and see whether they are consistent with your hunch. (I'm uncertain whether they are.) Perhaps Mathematica only stores the results of the final outcome. In this case the backwards request would receive no benefits.

---

Szabolcs' suggestion

Here are some orderings recommended by Szabolcs. Notice how the orderings increasingly move toward canonical ordering from least to greatest.

Table[{k, First@Timing[
 Prime /@ Flatten@Reverse@Partition[Range[1, 20], k]]}, {k, 1, 10}];
Table[{k, Prime /@ Flatten@Reverse@Partition[Range[1, 20], k]}, {k, 1,10}]

And their respective timings:

Table[{k, First@Timing[Prime /@ Flatten@Reverse@Partition[Range[1, 20000], k];]}, {k, 1,10}]

{{1, 2.540567}, {2, 1.275608}, {3, 0.858899}, {4, 0.647391}, {5, 0.517372}, {6, 0.431639}, {7, 0.370851}, {8, 0.325734}, {9, 0.291674}, {10, 0.262256}}

Here's the canonical ordering:

First@Timing[Prime /@ Flatten@Partition[Range[1, 20000], 1];]

0.011894