How to Find the Next Largest Known Prime

The Electronic Frontier Foundation (EFF) is offering a $150,000 cash prize “to the first individual or group who discovers a prime number with at least 100,000,000 decimal digits” and a $250,000 cash prize “to the first individual or group who discovers a prime number with at least 1,000,000,000 decimal digits” (https://www.eff.org/awards/coop).  The current largest known prime number, 2^74,207,281 – 1, is just over 22 million digits long and was found on January 7, 2016 by mathematics professor Curtis Cooper using software from GIMPS (the Great International Mersenne Prime Search, http://www.mersenne.org/), an organization where private users can help calculations by donating computer processing time during off-hours via the internet.

A prime number is a number that is only divisible by 1 and itself, and a number that is not prime is called a composite number.  For example, the number 17 is a prime number because it is only divisible by 1 and 17 and no other numbers.  The first few prime numbers are 2, 3, 5, 7, 11, 13, 17, 19, and so on. 

In fact, back in ancient times Euclid proved that the list of prime numbers never ends.  If you suppose that there is a final prime number P, then there exists a number q which equals all the primes multiplied together plus one (q = p₁·p₂·p₃·…·P + 1).  Now, q cannot have a factor because if it did it would divide into q (by definition) and p₁·p₂·p₃·…·P (which has all the primes) and therefore also into 1 (since q = p₁·p₂·p₃·…·P + 1), which is impossible.  So q must be a prime number larger than P, but this contradicts the original assumption that P is the largest prime.  Therefore, there can be no final prime number P, which means the list of prime numbers never ends.  It also means that no matter how big the largest known prime is, there is always a bigger one waiting to be discovered.

There are also no known patterns for the distribution of prime numbers, but there are some tests.  The most basic way to test if a number n is a prime number is to use trial division, which divides the number n by all prime numbers p between 2 and √n.  If any p divides into the number n evenly, then the number n is not a prime number; otherwise if no p divides into the number n evenly, then the number n is a prime number.  For example, if we were to test n = 17 to be prime, we would divide 17 by all the prime numbers between 2 and √17 ≈ 4.1, or by 2 and by 3.  Since neither 2 nor 3 divide evenly into 17, 17 is a prime number. 

Unfortunately, trial division can become cumbersome very quickly with large numbers, because you have to know (and divide by) a list of all the prime numbers between 2 and √n.  Fortunately, there is a more efficient method called the Lucas-Lehmer test which can check whether a number is a Mersenne prime or not (and is also the method used by GIMPS).  A Mersenne prime, named after the French mathematician Marin Mersenne, is a prime number that can be expressed in the form of 2^p – 1.  For example, the number 7 is a Mersenne prime for p = 3, because 2³ – 1 = 7.  The first few Mersenne primes are 3 (p = 2), 7 (p = 3), 31 (p = 5), 127 (p = 7), 8191 (p = 13), and so on.  (It may be tempting to conclude that a Mersenne prime is formed whenever p is a prime number, but this is not true for the prime p = 11, because 2¹¹ – 1 = 2047 = 23·89 which is not prime.  It is true, however, that p must be a prime number for all Mersenne primes.)

The Lucas-Lehmer test uses the recursive sequence s₀ = 4 and s_n = s_n-1² – 2, and states that for all p > 2, s_p-2 is divisible by 2^p – 1 if and only if 2^p – 1 is a prime number.  For example, if we were to test p = 5, where 2^p – 1 = 2⁵ – 1 = 31, we would need to find s_p-2 = s_5-2 = s₃ and see if it is divisible by 31.  So s₀ = 4, s₁ = 4² – 2 = 14, s₂ = 14² – 2 = 194, and s₃ = 194² – 2 = 37634, which is divisible by 31, so 2⁵ – 1 = 31 is a prime number.  Since we are testing divisibility, the algorithm can be made even more efficient by using a modulus function, where s_n = (s_n-1² – 2) mod 2^p – 1, which means the Lucas-Lehmer test can be restated as s_p-2 ≡ 0 (mod 2^p – 1) if and only if 2^p – 1 is a prime number.  So to test p = 5 again, where 2^p – 1 = 31, the sequence would be s₀ = 4, s₁ = (4² – 2) mod 31 = 14, s₂ = (14² – 2) mod 31 = 194 mod 31 = 8, and s₃ = (8² – 2) mod 31 = 62 mod 31 = 0, so 2⁵ – 1 = 31 is a prime number.

Although the Lucas-Lehmer test is difficulty to prove (see here for the proof), it is very easy to program in a computer:

01	# Lucas-Lehmer Mersenne Prime Tester
02	# Python 2.7.3
03	# After running, type "M(p)" where p is the Mersenne Prime exponent in 2^p - 1.
04	#  For example, if you would like to test Mp = 2^7 - 1 = 127, type "M(7)".
05
06	# import python libraries
07	from time import time, strftime
08	import datetime
09
10	# M function
11	def M(p):
12	# make sure p > 2, a restriction for the Lucas-Lehmer algorithm
13	if p > 2:
14	# start timer
15	timestart = time()
16	# find Mp
17	Mp = 2 ** p - 1
18	# Lucas-Lehmer algorithm
19	s = 4
20	for x in range(0, p - 2):
21	s = (s * s - 2) % Mp
22	# display results
23	if s == 0:
24	print "2^" + str(p) + " - 1 is PRIME"
25	else:
26	print "2^" + str(p) + " - 1 is COMPOSITE"
27	#display elapsed time
28	timeelapsedint = round(time() - timestart, 2)
29	timeelapsedstr = str(datetime.timedelta(seconds = round(
30	timeelapsedint, 0)))
31	print "runtime: " + timeelapsedstr + " or " + str(
32	timeelapsedint) + " seconds."
33	else:
34	# p <= 2, so display a message to ask for p > 2
35	print "Please enter a number larger than 2."

Here’s what happens when you use this program to test a somewhat large Mersenne prime, 2^11,213 – 1, which has approximately 3,000 digits, on a dual-core 1.67 GHz processor (a mediocre-speed laptop for 2016):

>> M(11213) 2^11213 - 1 is PRIME runtime: 0:00:55 or 54.64 seconds.

So with less than 40 lines of code, you can use your computer to verify that a 3,000 digit number is prime in less than one minute!

Not bad, but this is still a far cry from the current record of 22 million digits, and an even further cry from the future cash-prize-awarding prime of 100 million digits.  In any case, finding the next largest prime number will require an educated guess, a fast computer (or more), and lots and lots of patience.

Proof for the Lucas-Lehmer Primality Test

The Lucas-Lehmer primality test states that for the recursive sequence s₀ = 4 and s_n = s_n-1² – 2, and for all p > 2, s_p-2 is divisible by 2^p – 1 if and only if 2^p – 1 is a prime number.  (The following is loosely based on the proof found here with some improvements.)

Definitions

s₀ = 4

Let a = ½s₀ = ½(4) = 2

Let b = a² – 1 = 2² – 1 = 3

Let u = a + √b

Let ū = a – √b

Then uū = (a + √b)(a – √b) = a² – b = a² – (a² – 1) = a² – a² + 1 = 1

Then s_n = u^{2^n} + ū^{2^n} by induction:

s₀ = u^{2^0} + ū^{2^0} = u¹ + ū¹ = u + ū = a + √b + a – √b = 2a = 2(½s₀) = s₀

Assuming s_n = u^{2^n} + ū^{2^n}:

s_n+1 = s_n² – 2

s_n+1 = (u^{2^n} + ū^{2^n})² – 2

s_n+1 = u^{2(2^n)} + 2u^{2^n}ū^{2^n} + ū^{2(2^n)} – 2

s_n+1 = u^{2^(n+1)} + 2(uū)^{2^n} + ū^{2^(n+1)} – 2

s_n+1 = u^{2^(n+1)} + 2(1)^{2^n} + ū^{2^(n+1)} – 2

s_n+1 = u^{2^(n+1)} + 2 + ū^{2^(n+1)} – 2

s_n+1 = u^{2^(n+1)} + ū^{2^(n+1)}

Proof of Sufficiency

Prove: If s_p-2 is divisible by 2^p – 1, then 2^p – 1 is a prime number.

Assume s_p-2 is divisible by 2^p – 1:

Let M = 2^p – 1

s_p–2 = kM

u^{2^(p-2)} + ū^{2^(p-2)} = kM

u^{2^(p-2)} = kM – ū^{2^(p-2)}

u^{2^(p-2)}u^{2^(p-2)} = kMu^{2^(p-2)} – ū^{2^(p-2)}u^{2^(p-2)}

(u^{2^(p-2)})² = kMu^{2^(p-2)} – (uū)^{2^(p-2)}

u^{2^(p-1)} = kMu^{2^(p-2)} – 1

Assume that 2^p – 1 is not a prime number:

Let q be the smallest factor of M = 2^p – 1

M ≡ 0 (mod q)

u^{2^(p-1)} ≡ k(0)u^{2^(p-2)} – 1 (mod q)

u^{2^(p-1)} ≡ -1 (mod q)

(u^{2^(p-1)})² ≡ (-1)² (mod q)

u^{2^p} ≡ 1 (mod q)

Since u^{2^(p-1)} ≠ 1, the order does not divide 2p – 1, so the order is 2^p

Let X = {n + m√b | n, m ε Z_q} and X* = {all elements of X except 0}

|X| = q² and |X*| = q² – 1 (as long as b = a² – 1 is not a perfect square)

u is in X*, and since the order of an element is at most the order of its size,

so 2^p ≤ q² – 1 < q²

Since q is the smallest factor of M, q² < M = 2^p – 1

But then 2^p < 2^p – 1 which is a contradiction

So 2^{p – 1} must be a prime number.

Proof of Necessity

Prove: If 2^p – 1 is a prime number, then s_p-2 is divisible by 2^p – 1.

Let M = 2^p – 1

Assume p is not a prime number

Then p has two factors m ≥ 2 and n ≥ 2 such that p = mn

But then 2^p – 1 = 2^mn – 1 which is divisible by 2^m – 1 and 2ⁿ – 1

This contradicts the assumption that 2^p – 1 is a prime number

So p must be a prime number

(-1)^(M-1)/2 ≡ -1 (mod M)

the exponent (M – 1)/2 = ½(2^p – 1 – 1) = ½(2^p – 2) = 2^p-1 – 1 must be odd

-1 to an odd exponent is -1

so (-1)^(M-1)/2 ≡ -1 (mod M)

2^(M-1)/2 ≡ 1 (mod M)

by Fermat’s Little Theorem, 2^p ≡ 2 (mod p), so 2^p – 2 = np

2^p – 2 must be even, and p (as a prime larger than 2) must be odd, so n = 2k and 2^p – 2 = 2kp

so the exponent (M – 1)/2 = ½(2^p – 1 – 1) = ½(2^p – 2) = kp

now 2^p ≡ 1 (mod 2^p – 1) or 2^p ≡ 1 (mod M)

so 2^(M-1)/2 = 2^kp = (2^p)^k ≡ 1^k (mod M) = 1 (mod M)

so 2^(M-1)/2 ≡ 1 (mod M)

3^(M-1)/2 ≡ -1 (mod M)

2²ⁿ⁺¹ – 1 ≡ 7 (mod 12) for all positive integers of n by induction:

2²⁽¹⁾⁺¹ – 1 = 2³ – 1 = 7 ≡ 7 (mod 12)

Assuming 2²ⁿ⁺¹ – 1 ≡ 7 (mod 12):

then 2²ⁿ⁺¹ ≡ 8 (mod 12)

2^{2(n + 1) +1} – 1 = 2^{2n + 3} – 1 = 2^{2n + 1 + 2} – 1 = 2^{2n + 1}2² – 1 ≡ 8(4) – 1 (mod 12) ≡ 7 (mod 12)

Since p is a prime larger than 2, p = 2n + 1 for all positive integers,

so M = 2^p – 1 = 2²ⁿ⁺¹ – 1 ≡ 7 (mod 12)

Since M ≡ 7 (mod 12), M = 12n + 7 for some integer n

By the law of quadratic reciprocity, (3 | 12n + 7)(12n + 7 | 3)

= (-1)^{½(3 – 1)½(12n + 7 – 1)} = (-1)^{6n + 3} = -1

Now, (12n + 7 | 3) = (3·4n + 3·2 + 1 | 3) ≡ (1 | 3) = 1

So (3 | 12n + 7) = -1

Or (3 | M) = -1, which means 3 is a quadratic nonresidue of M

So by Euler’s criterion, 3^(M-1)/2 ≡ -1 (mod M)

u^(M+1)/2 ≡ -1 (mod M)

u = a + √b = ^{(a + √b)2(a + 1)}/_{2(a + 1)} = ^{(a + 1 + √b)^2}/_{2(a + 1)}

because the numerator (a + √b)2(a + 1) = 2(a + 1)a + 2(a + 1)√b

= 2a² + 2a + 2(a + 1)√b = a² + a² + 2a + 1 – 1 + 2(a + 1)√b

= a² + 2a + 1 + 2(a + 1)√b + a² – 1 = (a + 1)² + 2(a + 1)√b + b = (a + 1 + √b)²

(recall that b = a² – 1)

so u^(M+1)/2 = [^{(a + 1 + √b)^2}/_{2(a + 1)}]^(M+1)/2 ≡ ^{2(a + 1)}/_{-2(a + 1)} (mod M) = -1 (mod M)

the numerator (a + 1 + √b)^2(M+1)/2 = (a + 1 + √b)^(M+1)

= (a + 1 + √b)(a + 1 + √b)^M ≡ -2(a + 1) (mod M)

using the binomial theorem in a finite field and prime exponent M,

(a + 1 + √b)^M ≡ (a + 1)^M + √b^M (mod M)

and using Fermat’s Little Theorem, (a + 1)^M ≡ a + 1 (mod M)

and √b^M = b^M/2 = b^{M/2 – ½ + ½} = b^{(M – 1)/2 + ½} = b^{(M – 1)/2}b^½ = (a² – 1)^{(M – 1)/2}b^½

= ((a – 1)(a + 1))^{(M – 1)/2}b^½ = (a – 1)^{(M – 1)/2}(a + 1)^{(M – 1)/2}b^½

since a = 2, (a – 1)^{(M – 1)/2} = (2 – 1)^{(M – 1)/2} = 1^{(M – 1)/2} ≡ 1 (mod M), and (a + 1)^{(M – 1)/2} = (2 + 1)^{(M – 1)/2} = 3^{(M – 1)/2} ≡ -1 (mod M)

so √b^M = (a – 1)^{(M – 1)/2}(a + 1)^{(M – 1)/2}b^½ ≡ (1)(-1)b½ (mod M) = -√b (mod M)

therefore, (a + 1 + √b)^M ≡ (a + 1 – √b) (mod M)

so (a + 1 + √b)^2(M+1)/2 = (a + 1 + √b)(a + 1 + √b)^M ≡ (a + 1 + √b)(a + 1 – √b) (mod M)

= (a + 1)² – b (mod M) = (a + 1)² – (a² – 1) (mod M) = a² + 2a + 1 – a² + 1 (mod M)

= 2a + 2 (mod M) = 2(a + 1) (mod M)

the denominator (2(a + 1))^(M+1)/2 = (2(a + 1))^{(M–1)/2 + 1} = (2(a + 1))^(M–1)/2(2(a + 1))

= 2^(M–1)/2(a + 1)^(M–1)/2(2(a + 1)) ≡ (1)(-1)(2(a + 1)) (mod M) = -2(a + 1) (mod M)

so u^(M+1)/2 ≡ -1 (mod M)

Therefore,

u^(M+1)/2 ≡ -1 (mod M)

u^{(M+1)/4+(M+1)/4} ≡ -1 (mod M)

u^(M+1)/4u^(M+1)/4 ≡ -1 (mod M)

u^(M+1)/4u^(M+1)/4ū^(M+1)/4 ≡ -ū^(M+1)/4 (mod M)

u^(M+1)/4(uū)^(M+1)/4 ≡ -ū^(M+1)/4 (mod M)

u^(M+1)/4(1)^(M+1)/4 ≡ -ū^(M+1)/4 (mod M)

u^(M+1)/4 ≡ -ū^(M+1)/4 (mod M)

u^(M+1)/4 + ū^(M+1)/4 ≡ 0 (mod M)

u^{(2^p–1+1)/4} + ū^{(2^p–1+1)/4} ≡ 0 (mod M)

u^{(2^p)/4} + ū^{(2^p)/4} ≡ 0 (mod M)

u^{2^(p–2)} + ū^{2^(p–2)} ≡ 0 (mod M)

s_p-2 ≡ 0 (mod M)

which means s_p-2 is divisible by 2^p – 1

Conclusion

We have proved the conditional (if s_p-2 is divisible by 2^p – 1, then 2^p – 1 is a prime number) and its converse (if 2^p – 1 is a prime number, then s_p-2 is divisible by 2^p – 1), which means the biconditional (s_p-2 is divisible by 2^p – 1 if and only if 2^p – 1 is a prime number) is also true.

Interesting Notes

Other s₀ values can be chosen, as long as the following criteria are met:

● b = a² – 1 is not a perfect square

● (a – 1)^{(M – 1)/2} ≡ 1 (mod M)

● (a + 1)^{(M – 1)/2} ≡ -1 (mod M)

The traditional value for s₀ is s₀ = 4 (a = 2, a – 1 = 1, a + 1 = 3, and b = 3), but other values that fit the above criteria are s₀ = 10 (a = 5, a – 1 = 4, a + 1 = 6, and b = 24) and s₀ = ²/₃ (a = ¹/₃, a – 1 = -²/₃, a + 1 = ⁴/₃, and b = -⁸/₉).

Using Your Computer to Calculate Pi

On March 21, 2015, Rajveer Meena recited 70,000 digits of pi from memory in just under 10 hours in Vellore, India, for the Guinness World Record of most digits of pi memorized.  Apart from utter amazement by this feat of memory, have you ever wondered how mathematicians can even calculate such a long list numbers for pi?

There are actually several methods for calculating pi.  You might recall from your math classes that pi is the ratio between the circumference and diameter of any circle, and so one way to calculate pi is to construct a circle, directly measure its circumference and diameter, and divide the two.  Unfortunately, this is only precise enough to obtain a few decimal places of pi.  Another way to calculate pi is to find the perimeter or areas of inscribed or circumscribed polygons of circles.  This was the technique favored by mathematicians before the invention of calculators and computers, but its precision is still limited to “just” a few hundred digits of pi.

inscribed octagon

Today, computers can be used to calculate millions of digits of pi by using infinite series formulas.  There are several different ways to do this, but one of the most efficient methods that is relatively straightforward is to combine Machin’s formula:

with the arc-cotangent power series formula:

The computer program for this is as follows:

01	# Pi Calculator
02	# Python 2.7.3
03	# After running, type "pi(n)" where n is the number of decimals for pi.  For
04	#  example, if you would like to calculate 100 decimals for pi, type "pi(100)".
05
06	# import python libraries
07	from decimal import Decimal, getcontext
08	from time import time, strftime
09	import datetime
10
11	# arccot function using power formula arccot = 1/x - 1/(3x^3) + 1/(5x^5) ...
12	def arccot(x, digits):
13	# set precision and starting values
14	getcontext().prec = digits
15	total = 0
16	n = 1
17	# loop while new term is large enough to actually change the total
18	while Decimal((2 * n - 1) * x ** (2 * n - 1)) < Decimal(10 ** digits):
19	# find value of new term
20	term = ((-1)**(n - 1)) * 1 / Decimal((2 * n - 1) * x ** (2 * n - 1))
21	# add the new term to the total
22	total += term
23	# next n
24	n += 1
25	# return the sum
26	return total
27
28	# pi function
29	def pi(decimals):
30	# start timer
31	timestart = time()
32
33	# find pi using Machin's Formula pi = 4 * (4 * arccot(5) - arccot(239))
34	#  and the power formula for arccot (see arccot function above)
35	print "pi = " + str(Decimal(4 * (4 * arccot(5, decimals + 3) - arccot(239,
36	decimals + 3))).quantize(Decimal(10) ** (-decimals)))
37
38	# display elapsed time
39	timeelapsedint = round(time() - timestart, 2)
40	timeelapsedstr = str(datetime.timedelta(seconds = round(
41	timeelapsedint, 0)))
42	print "runtime: " + timeelapsedstr + " or " + str(
43	timeelapsedint) + " seconds."

Here’s what happens when you use the program to find the first 1,000 decimal places for pi on a dual-core 1.67 GHz processor (a mediocre-speed laptop for 2016):

>> pi(1000) pi = 3.14159265358979323846264338327950288419716939937510 582097494459230781640628620899862803482534211706798214808 651328230664709384460955058223172535940812848111745028410 270193852110555964462294895493038196442881097566593344612 847564823378678316527120190914564856692346034861045432664 821339360726024914127372458700660631558817488152092096282 925409171536436789259036001133053054882046652138414695194 151160943305727036575959195309218611738193261179310511854 807446237996274956735188575272489122793818301194912983367 336244065664308602139494639522473719070217986094370277053 921717629317675238467481846766940513200056812714526356082 778577134275778960917363717872146844090122495343014654958 537105079227968925892354201995611212902196086403441815981 362977477130996051870721134999999837297804995105973173281 609631859502445945534690830264252230825334468503526193118 817101000313783875288658753320838142061717766914730359825 349042875546873115956286388235378759375195778185778053217 12268066130019278766111959092164201989 runtime: 0:00:04 or 3.73 seconds.

  

So with less than 50 lines of code, you can have your computer calculate 1,000 decimal places of pi in less than 4 seconds!  At this rate, you will have 1,000,000 decimal places of pi in just over an hour!  Not bad, considering the very first computer, the ENIAC, took 70 hours to calculate pi to 2,037 decimal places in 1949.  Computers have come a long way since then!

With better computers and more efficient (but more complicated) infinite series formulas, pi has been calculated to over 13.3 trillion digits!  This number is unfathomable to most of us.  If you were to recite this many digits of pi at the same speed as the 2015 world record holder Meena (about 2 digits per second), it would take you over 200,000 years to finish!  So why bother to find so many digits of pi?   Most mathematicians will give the same reason that mountain climbers use for climbing Mount Everest: because it’s there.