Primality Testing and Sieves

Math 242 Modern Computational Math

Primality testing by trial division

Yesterday we wrote a function, myPrimeQ, which tests whether a number is prime.

myPrimeQ:
	Argument: an integer n
	Return: True if n is prime, False if n is not prime

For example:
	myPrimeQ[11] ⟶ True
	myPrimeQ[12] ⟶ False

Here is one implementation from yesterday:

myPrimeQ[n_]:=Module[{foundDivisor=False},
Do[
If[Divisible[n,k],foundDivisor=True],
{k,2,Sqrt[n]}
];

Return[!foundDivisor]  (* return value *)
];

Try it out:

myPrimeQ/@{2,3,19,22,25,29}

This implementation isn’t very efficient, since the loop doesn’t stop when it finds a divisor.

A faster implementation:

Here, we use a while loop that stops when a divisor is found.

myPrimeQ2[n_]:=Module[ 
{foundDivisor=False,
d=2,
m=Sqrt[n]},

While[!foundDivisor && d≤m,
If[Divisible[n,d],foundDivisor=True];
d++
];

!foundDivisor (* return value *)
]

Try it out:

myPrimeQ2/@{2,3,19,22,25,29}

Note that this implementation is much faster than the previous when n has a small factor:

Timing[myPrimeQ[44466501336143]]

Timing[myPrimeQ2[44466501336143]]

A sophisticated implementation:

This implementation uses Throw to stop a Do loop when a divisor is found; in this case, the enclosing Catch takes the value False. If no divisor is found, then the Catch statement takes the value True when the loop finishes. The value of the Catch statement is also the value that the Module returns.

myPrimeQ3[n_Integer]:=Module[{m=Floor[Sqrt[n]]},
Catch[
Do[ 
If[Divisible[n,d],Throw[False]],
{d,2,m}
];
True
] (* the Module returns the "caught" value *)
]

Try it out:

myPrimeQ3/@{2,3,19,22,25,29}

This is a bit faster than the previous implementation:

Timing[myPrimeQ2[79009213931]]

Timing[myPrimeQ3[79009213931]]

Timing[myPrimeQ3[44466501336143]]

Listing the Primes

Your next task is to use your function, myPrimeQ, to build a function myPrimes that will create a list of all primes up to some value n.

myPrimes:
	Argument: an integer n
	Return: a list of primes less than or equal to n

For example:  myPrimes[20] ⟶ {2, 3, 5, 7, 11, 13, 17, 19}

myPrimes[n_]:=Module[{list={}},
Do[
If[myPrimeQ3[k],AppendTo[list,k]]
,{k,2,n}];
list
]

myPrimes[100]

Runtime

How efficient is your function? Find the runtime of your function for several choices of n. Make a plot of the time to compute this list. Potential values to test include n=(2^k) for k=1,2,3,.... Hint: Use ListPlot or ListLinePlot.

Based on your data, about how long would your function take to compute all the primes less than 1,000,000,000?

Timing[myPrimes[1000]][[1]]

Timing[myPrimes[10000]][[1]]

Timing[myPrimes[100000]][[1]]

Timing[myPrimes[1000000]][[1]]

Sieve of Eratosthenes

The Sieve of Eratosthenes starts with a big list of numbers 2,3,4,...,n. The algorithm marks 2 as prime, and then marks all multiples of 2 as composite. Then the algorithm iteratively marks the next unmarked number as prime, and all of its multiples as composite.

Implement the Sieve of Eratosthenes.

sieveErat:
	Argument: an integer n
	Return: a list of primes less than or equal to n

For example:  sieveErat[20] ⟶ {2, 3, 5, 7, 11, 13, 17, 19}

sieveErat[n_]:=Module[

]

Try it out:

sieveErat[1000]

Runtime

Analyze the time required to compute primes up to n using sieveErat. How does it compare with your previous method?

Sieve of Sundaram

Start with the list of integers 1,2,...,n. For every pair of integers (i,j) satisfying 1≤i≤j and i+j+2ij≤n, remove the number i+j+2ij from the list. Take all the numbers still in the list: double them and add 1 to each. This results in a list of all the odd primes below 2n+2.

Implement the Sieve of Sundaram.

sieveSund:
	Argument: an integer n
	Return: a list of primes less than or equal to 2n+2

For example:  sieveSund[20] ⟶ {2, 3, 5, 7, 11, 13, 17, 19, 23, 29, 31, 37, 41}