Skip to main content

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



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



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 .


Comments