Primality Tests

Math 242 Modern Computational Math, spring 2021

Primality testing by trial division

Our task from Monday was to write a function isPrime which tests whether a number is prime.

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

Here is one possible implementation:

isPrime[n_Integer]:=Module[{},
Do[
If[Divisible[n,k],Return[False, Module]], 
{k,2,Sqrt[n]}];
Return[True]
]

Note that Return works differently in Mathematica than in many other languages. The second argument Module, causes the return to quit not only the loop but also the module.

Try it out:

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

Another implementation:

This implementation does not use a Return statement. Instead, it uses a variable to remember whether or not a divisor has been found.

isPrime2[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:

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

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.

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

Try it out:

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

This seems to be a bit faster than the other implementations:

Timing[isPrime[79009213931]]

Timing[isPrime2[79009213931]]

Timing[isPrime3[79009213931]]

Listing the primes

Next we want to write a function listPrimes that will create a list of all primes up to some value n.

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

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

Here is the module framework. We just need to fill in the details.

listPrimes[n_]:=Module[{list={}},
Do[
If[isPrime[k], AppendTo[list,k]]
,{k,2,n}];
list (* return the list of primes *)
]

listPrimes[50]

Timing[listPrimes[100000]][[1]]

Make a function that returns the time required to find a list of primes up to n:

primetime[n_]:=Timing[listPrimes[n]][[1]];

runtimes=Table[{m,primetime[m]},{m,10000,100000,10000}]

ListPlot[runtimes]