From: bruck@mtha.usc.edu (Ronald Bruck) Newsgroups: sci.math Subject: Re: nPrimes Date: 3 Feb 1995 14:55:28 -0800 In article <3gu6cd$fil@ixnews3.ix.netcom.com> HJSmith@ix.netcom.com (Harry Smith) writes: >In <3gsuen$fvc@kitten.umdc.umu.se> Niklas Frykholm writes: > >> >>schneken@Informatik.TU-Muenchen.DE (Thomas Schnekenburger) wrote: >> >>> I put the exercise to compute nPrimes(n), >>> that is the number of primes smaller than n. >>> For example nPrimes(10)=4, nPrimes(100)=25, nPrimes(1000)=168, ... >> >>> Can someone give me any hints about fast algorithms for computing nPrimes() ? >> >>How about using the seive of Erastotenes (or however it is spelled). ... > >Try my program: > [LONG sieve program omitted] Let me begin by saying I'm not a number theorist, and what I'm about to say isn't of my own personal knowledge, but what I've been told by number theorists. Apparently the function \pi(x)--the number of primes <= x--can be computed very efficiently by using an asymptotic expansion. My understanding is that this is how Mathematica probably computes PrimePi[x] (as it calls the function). Math- ematica's response to Timing[PrimePi[1734871252]] is {0.1 second, 85803055}, and comes too fast to be done by sieving! By keeping a few precomputed values of PrimePi[x], and using a bisection algorithm, you can then find Prime[n] rather efficiently. This, at least, is what I've been told. But I, and I'm sure others, would appreciate hearing this from the horse's mouth, or at least an expert's newsreader response. --Ron Bruck ============================================================================== From: rusin@washington.math.niu.edu (Dave Rusin) Newsgroups: sci.math Subject: Re: nPrimes Date: 4 Feb 1995 11:48:32 GMT First schneken@Informatik.TU-Muenchen.DE (Thomas Schnekenburger) wrote: > I put the exercise to compute nPrimes(n), > that is the number of primes smaller than n. > For example nPrimes(10)=4, nPrimes(100)=25, nPrimes(1000)=168, ... > > Can someone give me any hints about fast algorithms for computing nPrimes() ? ...a perfectly reasonable question. Then In <3gsuen$fvc@kitten.umdc.umu.se> Niklas Frykholm writes: >How about using the seive of Erastotenes (or however it is spelled). ...and in a followup post, Harry Smith wrote: >Try my program: > >/* Start of file Prime8By.h ---------------------------------------------- */ >/* > Name = 'Prime8By - Generates Primes'; Both these responses are fine if you want to have all the primes, but if you just want the _number_ of primes this is way too slow. I don't have my references here at home (nor, at this hour, can they be fetched) but there have been numerous papers computing this function, typically (in recent years) in Math. of Computation. But the idea goes back to Chebyshev: Of the numbers less than N, half are even and so excluded; one-third are multiples of 3 and so excluded; one-sixth have been counted twice. So the number not divisible by 2 or 3 is about N-N/2-N/3+N/6. Continue in this way to account for successively more primes. This will give an approximation to (I'll use standard notation now) pi(N). If this looks like N.Prod(1-1/p), you're right; thinking of this as the reciprocal of the zeta function at (OK, near) s=1, you'll not be surprised to learn that greater control of the behaviour of pi(N) comes from a study of the zeta function and its zeros. This is part of why the Riemann hypothesis is considered so important. dave ============================================================================== From: cet1@cus.cam.ac.uk (Chris Thompson) Newsgroups: sci.math Subject: Re: nPrimes Date: 6 Feb 1995 03:05:38 GMT In article <3gvpig$sjq@mp.cs.niu.edu>, rusin@washington.math.niu.edu (Dave Rusin) writes: |> |> Both these responses are fine if you want to have all the primes, but if |> you just want the _number_ of primes this is way too slow. I don't have |> my references here at home (nor, at this hour, can they be fetched) but |> there have been numerous papers computing this function, typically (in recent |> years) in Math. of Computation. "Computing \pi(x): The Meissel-Lehmer Method" J.C. Lagarias, V.S. Miller, and A.M. Odlyzko Math. Comp. 44 (1985) 537-560 Required reading, and with a useful bibliography of previous work. But I think you are wrong to imply that much of the latter had appeared in Mathematics of Computation. Indeed, was there anything there between David Mapes' 1963 paper and the above? Or anything since? |> But the idea goes back to Chebyshev: [description of inclusion-exclusion omitted] Or even to Legendre. Chris Thompson Internet: cet1@cus.cam.ac.uk JANET: cet1@uk.ac.cam.cus ============================================================================== From: cet1@cus.cam.ac.uk (Chris Thompson) Newsgroups: sci.math Subject: Re: nPrimes Date: 8 Feb 1995 02:35:56 GMT In article <3guc90$7g3@mtha.usc.edu>, bruck@mtha.usc.edu (Ronald Bruck) writes: |> [snip] |> |> Apparently the function \pi(x)--the number of primes <= x--can be computed very |> efficiently by using an asymptotic expansion. There is a conditionally convergent series for \pi(x) in terms of the zeroes of the Riemann zeta-function. Andrew Odlyzko has turned this into an algorithm for computing \pi(x) in time O(x^{1/2+\epsilon}), which is theoretically the best known for large enough x. But it would be rather misleading to call it "efficient": the constants are not good, as far as I know the method has never actually been used to calculate exact values of \pi(x) for interestingly large values of x. The Meissel-Lehmer method mentioned in other postings gives an algorithm for computing \pi(x) in time O(x^{2/3+\epsilon}). If I have time, I will compose and post a sketch of how this works later this week. |> My understanding is that this is |> how Mathematica probably computes PrimePi[x] (as it calls the function). Math- |> ematica's response to Timing[PrimePi[1734871252]] is {0.1 second, 85803055}, |> and comes too fast to be done by sieving! |> |> By keeping a few precomputed values of PrimePi[x], and using a bisection |> algorithm, you can then find Prime[n] rather efficiently. I very much doubt that Mathematica is using either of the previously mentioned methods. More likely, as you imply, it is holding pre-calculated values of \pi(x), perhaps at intervals of something like 10^6 around your 1.73*10^9. Then to find \pi(1734871252)-\pi(1735000000) it would only need to sieve a bitmap for 128748 numbers against the 4358 primes <= 41651, which should be easy enough to do in 0.1 seconds on a modern workstation. But this is just speculation. Does anyone have any inside knowledge of how Mathematica actually computes PrimePi ? Chris Thompson Internet: cet1@cus.cam.ac.uk JANET: cet1@uk.ac.cam.cus ============================================================================== From: HJSmith@ix.netcom.com (Harry J. Smith) Newsgroups: sci.math Subject: Re: nPrimes Date: 10 Feb 1995 13:25:00 GMT In <3heka9$hc7@www.interramp.com> pp001000@interramp.com (Kurt Schneckenburger) writes: > >In article <3guc90$7g3@mtha.usc.edu>, bruck@mtha.usc.edu says... >> >>In article <3gu6cd$fil@ixnews3.ix.netcom.com> HJSmith@ix.netcom.com (Harry >Smith) writes: >>>In <3gsuen$fvc@kitten.umdc.umu.se> Niklas Frykholm >writes: >>> >>>> >>>>schneken@Informatik.TU-Muenchen.DE (Thomas Schnekenburger) wrote: >>>> >>>>> I put the exercise to compute nPrimes(n), >>>>> that is the number of primes smaller than n. >>>>> For example nPrimes(10)=4, nPrimes(100)=25, nPrimes(1000)=168, .. >>>> >>>>> Can someone give me any hints about fast algorithms for computing >nPrimes() ? >>>> >>>>How about using the seive of Erastotenes (or however it is spelled). > >> > >Along the lines of the above: > >The book "Introduction to Algorithms" by Thomas H. Cormen, Charles E. >Leiserson and Ronald L. Rivest (1993) on page 837 says there are 50,847,478 >prime numbers below 10^9 (1,000,000,000). > >The book "Projects in Scientific Computation" by Richard E.Crandall (1994) >on page 121 says that the number of primes for 10^4 is 1,229, 10^8 is >5,761,455, 10^12 is 37,607,912,018 etc. > >My program matches for 10^4 and 10^8. For 10^9 my program computed >50,847,533 using a basic sieve approach (stores all the primes into a >virtual array), a difference of 55 (additional). Anyone know of any >references citing similar statistics: how many primes below 10^6, 10^7, 10^9 >etc.? How about the number of primes for numbers between 10^8 and 10^9. >Any chance 50,847,478 in "Introduction to Algroithms" is a typo? > >- Kurt Schneckenburger > > Here or EMail at pp001000@interramp.com > > > The statement that there is 50,847,478 prime numbers below 10^9 (1,000,000,000) is a well known error that has appeared in many books. The correct number is 50,847,534 given in Paulo Ribenboim's book "The Little Book of Big Primes", Springer Verlag, New York. 1991. x pi(x) --------------------------------- 10^8 5,761,455 10^9 50,847,534 10^10 455,052,511 10^11 4,118,054,813 10^12 37,607,912,018 10^13 346,065,536,839 10^14 3,204,941,750,802 10^15 29,844,570,422,669 10^16 279,238,341,033,925 2 * 10^16 547,863,431,950,008 4 * 10^16 1,075,292,778,753,150 -- | Harry J. Smith | 19628 Via Monte Dr., Saratoga, CA 95070-4522, USA | Home Phone: 408 741-0406, Work Phone: 408 235-5088 (Voice Mail) | EMail: HJSmith@ix.netcom.com on the Internet via Netcom NetCruiser -- ============================================================================== From: cet1@cus.cam.ac.uk (Chris Thompson) Newsgroups: sci.math Subject: How the Meissel-Lehmer method computes \pi(x) in time x^{2/3+} Date: 13 Feb 1995 06:20:18 GMT Keywords: Meissel Lehmer Lagarias Miller Odlyzko prime pi complexity I promised to write a sketch for sci.math of how the Meissel-Lehmer method enables one to compute \pi(x), the number of primes <= x, in time O(x^{2/3+epsilon}) and store O(x^{1/3+epsilon}). First, note that this all comes from the paper Computing \pi(x): The Meissel-Lehmer Method J.C. Lagarias, V.S. Miller, & A.M. Odlyzko Math. Comp. 44 (1985) 537-560 This was the first application of the Meissel-Lehmer method in such a way as to get a time better than O(x^(1-epsilon)). The x^epsilon's in the time and store estimates can be replaced by explicit powers of (log x), and to see how to optimise these, and how to improve the implied constants, you should consult the paper. (Victor Miller tells me that an extra (log x) improvement has been made since that time.) The version presented here is a bare bones affair to show how the exponents of x arise. Secondly, note that the time and store estimates are precisely those for an Eratosthenes sieve running up to x^{2/3}. The blocks in such a sieve can be (in order of magnitude) the square root of the largest numbers being sieved, i.e. of the largest primes whose multiples are being sieved out, but not smaller. If the number of multiples of a given prime in the block became o(1), then the block start/stop overheads would dominate. Define f(x,n) to be the number of positive integers <= x which are not divisible by any of the first n primes p_1, p_2, ..., p_n. We intend to calculate f(x,A) where A=\pi(x^{1/3}). This is the number of positive integers <=x which are the products of primes > x^{1/3}; hence of exactly 0, 1 or 2 of them: f(x,A) = 1 /* don't forget to count 1 ! */ + (\pi(x) - A) + \sum_{A p_a. To see this, bear in mind that the first argument of f(.,.) must just have changed from >= x^{2/3} to < x^{2/3}. There would be O(x^{2/3}) terms in this sum even without the restriction on (a,n), so we are still on target for our time estimates. Now we are going to run a modified Eratosthenes sieve up to x^{2/3} to compute these type (a) terms. We use some of our O(x^{1/3+epsilon}) store to hold tables of \mu(n) and g(n) for n <= x^{1/3}, where g(n) is the i such that p_i is the smallest prime divisor of n, as well as a table of the p_i (1<=i<=A) themselves. We must perform the sieve for blocks y_min < y <= y_max, with y_max - y_min about x^{1/3}, to keep within our store limits. Yet another table of this size holds the f(y_min,n) for 1<=n<=A and will be updated to f(y_max,n) as part of processing the block. After sieving the block by the first a-1 primes, we will need to compute the values of f(x/(n*p_a),a-1) for n such that x/(y_max*p_a) <= n < x/(y_min*p_a), g(n) > a, and \mu(n) is non-zero. For this we need to be able to find f(y,a-1) for y anywhere in the range of the block. We can use a binary tree whose level m entries give the count of unsieved numbers in sub-blocks of length 2^m; then f(y,a-1) is computed as f(y_min,a-1) plus at most one count from each level of the tree, and so takes only time O(log x). To maintain this structure, sieving out a number must decrement one count at each level of the tree, and also takes time O(log x). The increment f(y_max,a-1)-f(y_min,a-1) also comes for free from this structure. It should now be clear how to proceed. For each block, we alternately sieve out multiples of a prime, and then compute the contributions of type (a) to f(x,A) with p_a the *next* prime to be sieved out. The sieving takes the same time as for a regular Eratosthenes sieve, apart from a logarithmic factor. The number of contributions to f(x,A), even counting the n's we skip over because they fail the tests on g(n) or \mu(n), is O(x^{2/3}). The start/stop overheads on the the ranges of n (i.e. the fact that we have to compute the ends of range even if it contains no integers) occur just as often as the start/stop overheads of the sieve (because they alternate) and will therefore also come to O(x^{2/3+epsilon}) provided the blocks are no smaller than x^{1/3}. I think it is time to go away and sleep after that... Chris Thompson Internet: cet1@cus.cam.ac.uk JANET: cet1@uk.ac.cam.cus ============================================================================== From: desnogue@unice.fr (Desnogues) Newsgroups: sci.math Subject: Re: How the Meissel-Lehmer method computes \pi(x) in time x^{2/3+} Date: 17 Feb 1995 14:35:32 GMT > (Victor Miller tells > me that an extra (log x) improvement has been made since that time.) This improvement, that yields an O(x^{2/3}/\log^2(x)) time complexity and an O(x^{1/3}\log^3(x)\log\log(x)) space complexity, is due to Marc Deleglise and Joel Rivat from the French University of Lyon 1. Using their improvement they computed \pi(x) upto 10^18. -- Laurent Desnogues