Take the 2-minute tour ×
Mathematica Stack Exchange is a question and answer site for users of Mathematica. It's 100% free, no registration required.

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 ?

share|improve this question
5  
Probably Prime[n] makes use of results for Prime[m], m < n, but it caches only a limited number of results. But this is just a guess and you have surely though of it already. –  Szabolcs Dec 5 '12 at 22:24
    
So with a clear memory, Prime[n] runs in O(n) time.. –  cartonn Dec 5 '12 at 22:34
    
Prime and PrimePi are Listable, they are automatically mapped over a list, thus you can simply write e.g. Prime @ Range @ m instead of Prime[#] & /@ Range[1, m]. –  Artes Jan 5 '13 at 18:03
add comment

2 Answers

up vote 6 down vote accepted

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 .

share|improve this answer
    
this helps a lot, thank you –  cartonn Dec 7 '12 at 3:50
add comment

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}]

ordering

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

share|improve this answer
    
r1 and r2 I think you meant.. –  cartonn Dec 5 '12 at 22:30
    
Yes, thanks. I corrected it. –  David Carraher Dec 5 '12 at 22:56
    
This kind of ordering is an interesting experiment too: Table[ {k, First@ Timing[Prime /@ Flatten@Reverse@Partition[Range[1, 20000], k];]}, {k, 1, 10} ] –  Szabolcs Dec 6 '12 at 0:12
    
Yes! Very nice. I'm comparing that output to First@Timing[Prime /@ Flatten@Partition[Range[1, 20000], 1];] –  David Carraher Dec 6 '12 at 0:22
    
One would expect the random ordering to be the slowest, very interesting.. –  cartonn Dec 6 '12 at 16:05
add comment

Your Answer

 
discard

By posting your answer, you agree to the privacy policy and terms of service.

Not the answer you're looking for? Browse other questions tagged or ask your own question.