Shameless Mathematica - Some Fun with Prime Numbers
Hareem Naveed

March 18, 2016

Credit: This project is inspired by some of the work we did with Fibonacci numbers. It mainly draws on the theory listed here: https://primes.utm.edu/prove/merged.html, and relies on “A Classical Introduction to Modern Number Theory” as a reference.

In this project, we implement a test to check for primality of large numbers, specifically, we check the primality of Fermat numbers. We then implement a heuristic search method to find large prime numbers. 

Fermat Numbers 
A Fermat Number is a number of the form (2^((2n)))+1, where n is a positive integer. First, let us write a short function and output the first few Fermat numbers:

fermatnumber[n_]:=2^(2^n)+1

fermatnumber[2]

Grid[
Prepend[
Table[{n,fermatnumber[n],PrimeQ[fermatnumber[n]]},{n,0,10}],{"n","F_n","Prime"}],Frame -> All]

Only the first four Fermat numbers are prime, and as we can see, Fermat numbers grow very rapidly! At this point, it is noted that using current computational methods, it is unlikely that other prime Fermat numbers will be found. Note from the above table that should another prime Fermat number be found, it would be a very large number and it would be very difficult to check it’s primality.

For example, Mathematica has trouble computing the factors for even the 10th Fermat number, and times out of the computation. Note that the it only takes 1.06 seconds to compute the factors of the 8th Fermat number:

FactorInteger[fermatnumber[8]] // Timing

In this project, we will explore both Fermat numbers and work with some large primes. 

First, we will implement Pépin’s test, which is a test to check the primality of a Fermat number. Second, we will work with Mathematica’s inbuilt “Prime[n]” function which returns the n’th prime, and implement a simple methods to find large prime numbers. 

First, let us update the above table with timing for each of the evaluations:

Grid[Prepend[Table[{n,fermatnumber[n],Timing[PrimeQ[fermatnumber[n]]]},{n,0,13}],{"n","F_n","Time, Prime"}],Frame-> All]

In another format, note that when we plot the time for computation, it grows exponentially. We want an algorithm that accomplishes the same task in polynomial time.

ListLinePlot[Table[First[Timing[PrimeQ[fermatnumber[n]]]],{n,0,15}],PlotStyle->PointSize[Large]]

First[Timing[PrimeQ[fermatnumber[15]]]]

Pépin’s Test:
If we let (Subscript[F, n]) be the n’th Fermat number, (Subscript[F, n]) is prime iff (3^(((Subscript[F, n])-1)/2)) ≡ -1 (mod (Subscript[F, n])). We can implement this using the repeated-squaring technique in a “Do” function. Note that the base can similarly be 5 or 7.

a = 3;
Do[a = Mod[a^2,fermatnumber[3]],{n,1,2^3-1}]
Print[{a,fermatnumber[3],a-fermatnumber[3]}] // Timing

PrimeQ[fermatnumber[3]] // Timing

a = 7;
Do[a = Mod[a^2,fermatnumber[8]],{n,1,2^8-1}]
Print[{a,fermatnumber[8],a-fermatnumber[8]}] // Timing

PrimeQ[fermatnumber[8]] // Timing

Note that this method is only faster than the native “PrimeQ” function when working with larger Fermat numbers. Let us plot the output of calling this function for the different Fermat numbers and compare with the “PrimeQ” plot above. To do this, first, let us write a function that implements Pépin’s test.

pepinstest[n_]:=Module[{m,b},
m=fermatnumber[n];
b=3;
Do[b= Mod[b^2,m],{m,1,2^n-1}];If[Mod[b,m]==m-1,Return["prime"],Return["composite"]]]

pepinstest[8] // Timing

pepinstest[10] // Timing

PrimeQ[fermatnumber[10]]//Timing

The new method is slightly faster than the in-built “PrimeQ” function up to n = 10, but does much better than the “FactorInteger” function which times out before giving an answer.

pepintimes = Table[First[Timing[pepinstest[n]]],{n,0,15}]

plot1 = ListLinePlot[pepintimes,PlotStyle->Red,ImagePadding -> 25,Frame->{False,False,False,True},FrameTicks->{None,None,None,All},FrameStyle->{Automatic,Automatic,Automatic,Red}]

primetimes = Table[First[Timing[PrimeQ[fermatnumber[n]]]],{n,0,15}];

plot2 = ListLinePlot[primetimes,PlotStyle->Blue,ImagePadding -> 25,Axes-> False,Frame->{True,True,False,False},FrameTicks->{All,None,None,None},FrameStyle->{Automatic,Blue,Automatic,Automatic}]

Overlay[{plot1,plot2}]

pepintimes

primetimes

Note from the above computations that after (Subscript[F, 12]), the native “PrimeQ” function is much slower than Pépin’s test for evaluting the primality of the Fermat Number. For (Subscript[F, 15]), “PrimeQ” is 200 times slower than Pépin’s test. We can even evaluate the primality of (Subscript[F, 20]) relatively quickly using this method!

plot3 = ListLinePlot[Table[First[Timing[pepinstest[n]]],{n,0,20}]]

pepinstest[20]

Of course, we know that it will be a composite number, but the purpose of this exercise is to demonstrate the implementation of a trick (Pépin’s test) that is faster than the inbuilt computational methods in Mathematica for determining the primality of a large number.

Finding Large Prime Numbers

In Mathematica, you can specify the nth prime number with a simple command:

Prime[100]

Prime[10^8] // Timing

The list of largest known primes includes Mersenne primes, as the method used for testing primality is based on factorization of n+1 or n-1. For Mersenne primes, the factorization of n+1 is trivial as it is simply a power of 2. 

One way to find a large prime number is to start with a random prime number, then multiply it by another large random integer and add 1.

a= Prime[k=Random[Integer,{1,10^9}]]

b = Random[Integer,{1,10^8}]

possibleprime= a*b+1;

possibleprime

PowerMod[2,possibleprime-1,possibleprime]

This number is clearly not prime. We can write a loop that will check for the different values of b that will give us a prime number.

While[PowerMod[2,possibleprime-1,possibleprime]≠ 1 &&
PowerMod[3,possibleprime-1,possibleprime]≠ 1 &&
PowerMod[5,possibleprime-1,possibleprime]≠ 1,
b=Random[Integer,{1,10^15}];
possibleprime=a*b+1]

possibleprime

PrimeQ[possibleprime]

By changing the parameters of a and b, and the search space for the “Random” functions, we can specify larger and larger primes.