Counting Primes

Math 242 Modern Computational Mathematics

Counting Primes

The prime-counting function π(x) counts the number of primes less than or equal to a real number x. For example π(8)=4 because there are four primes less than or equal to 8. We will compute values of π(x) using a list of primes found by the Sieve of Eratosthenes. Here is a fast version of the sieve:

sieveErat[n_Integer]:=Module[{nums=Range[n],m=Floor[Sqrt[n]]},
nums[[1]]=0;
Do[
If[nums[[i]]≠0,nums[[2i;;;;i]]=0],
{i,2,m}];
Select[nums,#≠0 &]
]

One way to count the primes less than or equal to n is simply to find the length of the list of primes returned by sieveErat[n]. For example:

Length[sieveErat[17]]

This works, since π(17)=7. However, we would like to compute lots of values of π(x). Running the sieve of Eratosthenes for each value of x would be inefficient. Here is a more efficient method, implemented as a module countPrimes with the following specification.

Module: countPrimes
Arguments: positive integer n
Return: list containing values π(x), for x=1,2,3,...,n

The module will do the sieve of Eratosthenes only once, finding all primes up to n. Then the module will loop over the integers i from 2 through n, keeping track of how many primes it has found up to i.

countPrimes[n_Integer]:=Module[{primeList,c,i,counts=Table[0,n]},
primeList=sieveErat[n];
c=0; (* count of primes found so far *)
Do[
If[c<Length[primeList],
If[i==primeList[[c+1]], (* then i has reached the next prime *)
 c++
];
];
counts[[i]]=c  (* there are c primes ≤ i *)
,{i,2,n}];
counts (* return the list of counts *)
]

This is quite fast:

Timing[countPrimes[100000];]

Exercise 1

Verify that countPrimes correctly computes π(x) for a few small values of x. Then use ListPlot to plot π(x) for x from 1 to some positive integer m that you choose. What is the shape of the graph? Could you approximate the graph by a smooth function?

Density of Primes

In about 1793, Gauss conjectured that the density of primes near x is approximately 1/ln(x). Gauss had printed tables of prime numbers to work with. You have a computer and fast sieve algorithms. Can you provide evidence for or against Gauss’ conjecture?

To do this, you will need to approximate the density of primes using values of π(x). One way to do this is to fix an interval length d, and then compute the number of primes between x and x+d, divided by the length of the interval. That is, the density of primes near x is approximately:

((π(x + d) - π(x))/d).

Exercise 2

Plot an approximation of the prime density function. Can you provide evidence in support of Gauss’ conjecture about the density of primes? It may help to experiment with different interval sizes d: a larger value of d will help smooth out the density approximation.

Note that our function primeCounts only accepts integer input. If x is not an integer, then you might want to use primeCounts[Floor[x]].

Approximating the prime counting function

The prime counting function is not continuous. However, various continuous functions can be used to approximate the prime counting function. A simple approximation is given by:

(x/(ln(x)))

Another approximation is given by the integral of Gauss’ density estimate. This is known as the logarithmic integral:

li(x)=(∫_0^x)(dt/(ln t))

In Mathematica, the logarithmic integral is implemented by LogIntegral.

Exercise 3

Plot the prime counting function, together with the approximations given above, on different intervals. Do you think these are good approximation? What is the percent error for each approximation at various values of x? Does the percent error decrease when x gets large?

Twin Primes

A pair of primes that differ by 2 are called twin primes. For example, 11 and 13 are twin primes, as are 29 and 31.

The twin prime counting function is denoted (Subscript[π, 2])(x). Specifically, (Subscript[π, 2])(x) is the number of primes p≤x such that p+2 is also prime. The following function computes (Subscript[π, 2])(x).

Module: countTwinPrimes
Arguments: positive integer n
Return: list containing values (Subscript[π, 2])(x), for x=1,2,3,...,n

countTwinPrimes[n_Integer]:=Module[{primeList,j,c,counts=Table[0,n]},
primeList=sieveErat[n+2];
j=1; (* tracks position in primeList *)
c=0; (* counts twin primes found so far *)
Do[
If[j+1≤Length[primeList], (* avoid running off the end of the list *)
If[i==primeList[[j]], (* then i has reached the next prime *)
If[i+2==primeList[[j+1]],c++]; (* then (i, i+2) are twin primes *)
 j++
];
];
counts[[i]]=c  (* there are c twin primes (p, p+2) such that p ≤ i *)
,{i,2,n}];
counts  (* return the list of counts *)
]

Exercise 4

Make a graph of (Subscript[π, 2])(x). What do you notice? How does the graph of (Subscript[π, 2])(x) compare with the graph of π(x)?

Other Pairs of Primes

A pair of primes that differ by 4 are called cousin primes. A pair of primes that differ by 6 are called sexy primes (a pun on the Latin word for six). Make a copy of countTwinPrimes and modify it to count cousin primes or sexy primes.

Exercise 5

How do the counts of cousin primes and sexy primes compare to the counts of twin primes? Is one type of prime pair more common than the others?