The famous mathematician Bernhard Riemann did some novel and exciting work in 1859 on the old question of how many primes there are below a given number, which is usually written π n This function is equivalent to asking, "what is the value of the k-th prime for any given k?" The latter is actually a more specific question since an arbitrary n will in general fall between primes. Indeed, the notation is most curious because π refers to the number of primes up to n, not to any prime per se, while n actually refers to the values of the primes (since π changes by unity every time n passes through a prime value. Thus if n is a prime, it is the π-th prime, pretty much the inverse of defining a function p n where p is the value of the n-th prime. We will encounter this inversion as a practical matter later. However, it is a testiment to Riemann's genius that he derived closed expressions that gave the exact number of primes. For example, the 4th prime is 7, but if you chose n8,9,10, Riemann's formula will give the same answer (4) for each of these values for n. For this reason, the notation seems reversed, but actually makes sense. The PNT as a theorem asserts that the following interpolation formulas, which attempt to smoothly fit onto the trends in the distribution of the primes, are asymptotically exact. An early estimate was that π n was approximately n n This estimate systematically underestimates the numbers of primes, except possibly for staggeringly large values of n. As a theorem it would state that in the limit of n going to infinity, the ratio of π n to n n is unity. The theorem is of particular interest to mathematicians in that it follows from an as-yet unproven conjecture by Riemann on the distribution of zeros of the zeta function, a complicated topic which is not required for our simple analysis here. An improved approximation is Li n where Li, the log integral function is defined to be t 2 n 1 t This estimate is substantially better and systematically underestimates the number of primes (again until one gets to staggeringly large values). Notice that the first approximation results from considering the log term to be slowly varying and taking it out of the integral to get the (crude) estimate of n n . The log integral function diverges at n 1 , but this has no relevance to its applications at large n. Some large values (from Derbyshire, p. 116). n π Li π 10 8 5,761,455 754 10 9 50,847,534 1,701 10 10 455,052,511 3,104 10 11 4,118,054,813 11,588 10 12 37,607,912,018 38,263 10 13 346,065,536,839 108,971 10 14 3,204,941,750,802 314,890

Anyone used to examining errors for trends would see that the steady increase in the difference Li π by about a factor of 3 per decade in n suggests that there is a scaling correction that needs to be applied.

Calculating the primes The method of choice to calculate a list of primes is to use the "sieve" of Eratosthenes. Here one starts with a list the natural numbers. Then starting with 2, one removes every even number. The next number is 3, so now one removes every third number (half of these have already been removed by the 2). Here it is easiest to simply replace the values with zero, rather than literally removing the numbers. Now the next (non-zero) number is 5, and we set every fifth number to zero. Then 7 and then 11. By eleven, the first 121 non-zero numbers will all be primes (i.e., 11 squared). The alternative of testing each successive number to see if a smaller prime divided it would be hopelessly inefficient. Here is a simple MATLAB program (easily translated into the program of your choice):

Here one simply chooses the value of N, and after executing the program, the vector "primes" will contain the list: 23571113...97 (if N100). Note that the "sieve" itself never removes "1" and this number has all of the properties of a prime (not divisible by any other number besides 1 and itself). But it has no use in the unique factoring of natural (whole) numbers into primes, and is so excluded. This program has been checked against the list given in Abramowitz and Stegun (the primes up to 99,991, which are roughly the first 10,000 primes).

Derivation of a PNT The sieve suggests a simple way of estimating the distribution of primes. After any given prime, there will be a density of non-zero numbers remaining. For example, after 2, half the remaining numbers will be nonzero. After the 3, one removes 13 but 12 are already gone, so we get 12 16 13 left. Since there are an infinite number of numbers left, it is easier to think in terms of remaining blocks of numbers; after removing 2 and 3, we have blocks of 2 3 6 , and in any arbitrary block of 6 consecutive numbers beyond the 3, there will be exactly 2 non-zero numbers. If we go to 5, the blocks become 2 3 5 30 , and there will remain exactly 8 non-zero numbers in any block of 30 numbers beyond 5. It is easy to verify that this density evolves quite systematically. After each successive prime p, the average density per block drops by exactly a factor of p 1 p This density then corresponds to a mean distance between primes, since the next prime is chosen from the first non-zero number in the following block, and we can guess this distance from the mean density. In particular, after each successive prime p, the density decreases by p 1 p and the mean distance correspondingly increases by p p 1 . Let us call this mean distance after the k-th prime Δ k , while the k-th prime itself we would call P k . It follows then that the next prime will be approximately P k + 1 P k Δ k but given this new prime, the mean step size increases to Δ k + 1 Δ k P k + 1 P k + 1 1 But now we are all done. If we chose the first prime ( P 1 ) and the distance to the next prime ( Δ 1 ), we can simply iterate these relationships indefinetely. Indeed, we know these values: P 1 2 and Δ 1 1 A MATLAB program would simply be

We chose N 9500 because this available on the tabulated lists, or if you created the list of primes, you should get primes 9500 98,947 . In contrast, the above program gives 102,014, which is too large by only about 3%. This result is rather astonishingly good given that the initial primes and steps between primes are anything but regular.

Obtaining the PNT Results It is easy to rearrange the equations as difference equations and then write them as differential equations. First we write P k + 1 P k Δ k and Δ k + 1 Δ k Δ k 1 P k + 1 1 For large k, we can approximate these as k P k Δ k and k Δ k Δ k P k One might worry about simplifying P k + 1 1 to P k , but the "corrections" are tiny, and unimportant as we will show. The ratio is Δ P P which integrates to give P P 1 Δ Δ 1 which then gives k 1 N P P 1 P 1 Δ 1 P P 1 Notice that this is in fact the log integral function of the improved PNT. However, because of the curious choice of notation, the natural symbols for π (here N) and n (here P) are reversed. Notice that the Δ 1 in the denominator can be absorbed into the P P 1 and the rescaling P moves these terms to the limits on the integral. Effectively, this simply rescales Li n . Accordingly, we can chose a Δ 1 such that the curve is through the primes instead of being always off to one side. Such a rescaling makes perfect sense given that the PNT is just a fit to the primes and such a fit has to be global. Thus the fit is not obliged to either pass exactly through the prime 2 nor have the initial step size be exactly 1. Somehow the PNT derivation of Li n has effectively slipped in an implicit assumption that Δ 1 P 1 (e.g., Δ 1 2 0.6931 ...). In fact, choosing Δ 1 0.625 provides a better fit that at large k gives an estimate for the value of the primes that oscillates about the correct values, being off a few percent to either side. This oscillation shows that there is little point in trying to do "better" (e.g., trying to solve a more exact differential equation above; in fact we just iterated instead). A significant improvement could be gotten by letting Δ 1 be a weak function of k. Notice that we have not proven the PNT, we have just derived the same asymptotic functions cited in the proofs. One can create a number of interesting relationships using these approximate expressions and the PNT itself. For example, if we integrate the log integral function between two (large) successive primes, the result should be just unity (i.e., one new prime). Since the P k denominator will hardly change, the integral itself will be approximately 1 P k + 1 P k P k which can be rewritten as the recursion relation P k + 1 P k P k which must be one of the simplest possible recursion relations for primes. If we start it with the obviously simplest choice P 1 2 , we get .

The first thousand or so primes, plotting P k against k. The uppermost (ragged) line are the actual primes, the dashed line just below it is Li k , and the solid line is the simple recursion relation P k + 1 P k P k starting with P 1 2 . Another step is to replace P k + 1 P k with Δ k , namely the mean distance between primes. This gives us k k P k 1 1 P 1 1 1 P 2 ... 1 1 P k const and the result for the first million primes is plotted in . The constant is actually γ where gamma is Euler's constant, 0.5772157... and the exponential is 0.56145946... and we can see the convergence to this value. This result was derived in 1874 and is known as Merten's theorem (but the proof was based on Riemann's (1859) paper and not on a causal idea of a "mean distance between primes"). The above product form is a favorite of number theorists because its inverse is a special case of the zeta function, namely the value of that function at unity, which happens to be an infinity that just cancels the infinity in the limit of P k (in the sense of the above limit).

Plot of P k divided by Δ k . If we were to plot this function on ordinary linear scale, it would look like a step function starting at 0.346 and almost immediately jumping to what looks like the asymptotic value, and then being constant out to a value of almost 100,000. So we have expanded the vertical scale to show the trend and plotted the horizontal scale as a log plot to show how the "fluctuations" die off with increasing size of the primes. Here a linear plot suppresses that trend because almost all of the points to the right of the origin would now correspond to "large" primes. In Merten's theorem we have used the "natural" definition for Δ 1 of 12, whereas when we used it in the recursion relation for primes we instead defined Δ 1 1 , since that was the actual spacing between the starting prime 2 and the next prime, 3. Alternatively, to agree with the log integral expression we would have had to define Δ 1 2 . In recursion relations there is often such a freedom of choice, as well a "natural" choice.