This article originally attributed the last Mersenne prime found by hand to Leonhard Euler. It was Édouard Lucas. So, that’s what the title is about. We also only realised after publication that there’s a wordplay in the title that only works if you say Euler wrong. The proofreaders apologise.

There is an ongoing quest in mathematics to find larger and larger prime numbers. Since the primes go on forever, there is no single largest one, but every day around the world computers are churning to find a new record holder.

This project was inspired by a Numberphile video about how mathematicians find enormous prime numbers. The goal is to find a large prime number ourselves, one that took mathematicians over 2,000 years to find. The record for the largest known prime has since been broken since that video was released; the largest known prime as of May 2021 is

\[ M_{82,589,933} = 2^{82,589,933} – 1 \]

If you were to write the number down it would be longer than the entire text of the Harry Potter series. We won’t find any prime numbers nearly that big, but we can try to find one too big to fit in a single tweet.

2,000 years of math in 10 lines of Python

Mathematicians have been trying to find progressively larger prime numbers for several thousand years – the University of Tennessee at Martin has a great introduction to the topic. While the current record holders are found with computers, these techniques were once done by hand. The last time a human set the record for calculating the largest known prime was done by Édouard Luca in 1876 when he proved that $M_{127}$ is prime. Lucas’ record stood for 76 years until 1952 when a computer finally found a larger prime number. Since then all records have been set with the help of computers. In 10 lines of Python code, we will prove Lucas was right and find a bigger prime, surpassing more than 2,000 years of mathematical history.

With pen and paper, Lucas showed that $M_{127}$ is prime.

\[ M_{127} = 2^{127} – 1 = 170,141,183,460,469,231,731,687,303,715,884,105,727 \]

There’s a lot of Python code in this post. If you don’t have Python on your computer, here’s an online notebook you can use to run the code in your browser.

How do they set world records?

$M_{82,589,933}$ was found by a worldwide computer search called the Greater Internet Mersenne Prime Search, known as GIMPS. GIMPS works by picking a number that might be prime and testing it using a massive network of volunteered computers. By breaking the rather difficult question “Is this number prime?” into smaller pieces each machine in the network can work on a chunk of the problem they can handle, then report its results. This massively parallelized structure allows their network of more than two million computers to tackle prime numbers far larger than any single machine could handle.

GIMPS runs several tests on each candidate number before running the most conclusive, and slowest, test to see if a candidate number is indeed prime. These initial screens, which include trial division and Pollard’s P-1 method, can quickly show if a candidate is composite, so they can eliminate candidates, but cannot prove a number is prime. Once a candidate passes all the early rounds it’s time for the conclusive test. Until 2018, their conclusive test was the Lucas-Lehmer Primality test. The Lucas-Lehmer test is a definitive method to prove whether or not a number is prime and the one we’ll be focusing on. In 2018, GIMPs did change their final test to the Fermat Probable Prime test since it could be implemented with a lower chance of hardware error. They improved the method again in 2020 when new work by Krzysztof Pietrzak eliminated the need to double-check results. It’s exciting to see new techniques still being created to solve such an old problem, but we’ll stick with the Lucas-Lehmer test for now.

The Lucas-Lehmer Primality Test

Computers are helpful but the Lucas-Lehmer test was originally developed and done by hand. The method is so easy to write down it fits in 10 lines of Python code. Proving it works is not as straightforward, so we’ll understand a related piece of mathematics to get a feel for how the Lucas-Lehmer test works. Let’s focus on Fermat’s Little Theorem as a way to understand how the behavior of prime numbers can help us build a test to find them. We won’t prove anything here but I hope this provides a sense of how mathematicians can turn a property of a number into a way to detect it.

Fermat’s Little Theorem is a statement that shows one way in which prime numbers tie together exponents and modular arithmetic. It gives us an exponent that can be used to make any base equal to $1 \mod p$.

Fermat’s Little Theorem:

\[ a^{p-1} \equiv 1 \mod p \]

where $a$ is an integer and $p$ is a prime.

This feels very abstract, so let’s try the example from Wikipedia where $a=2$ and $p=7$.

\[ a^{p-1} \equiv 1 \mod p \implies a^{p-1} – 1 \equiv 0 \mod p \]

Keep in mind that $\equiv 0 \mod p$ is the same thing as saying a number is a multiple of $p$.

\[ 2^{7-1} – 1= 2^6 -1 = 64 -1= 63 \text{ is a multiple of 7, so } 63 \equiv 0 \mod 7 \]

The takeaway from this is to understand that raising any base to the power $p-1$ and having the result be equivalent to $1 \mod p$ is a property of every prime number, no matter how big. This means we know something about the prime numbers that have not yet been discovered. We use this knowledge to our advantage by changing how we use the equation. If we have a number that we hope is prime but aren’t sure, we can use Fermat’s Little Theorem to see if the candidate behaves like a prime.

Let’s consider a candidate, $n$, and see how we can use Fermat’s Little Theorem to test if the candidate behaves like a prime. If we have a candidate, $n$, and check that $a^{n-1} \equiv 1 \mod n$ for many values of $a$, then $n$ is behaving like a prime number. We have to be careful here because the converse of Fermat’s Little Theorem here is not true. In our example, $n$ is behaving like a prime number but we don’t know for sure if $n$ is a prime number. The Carmichael Numbers live in the murky gap of behaving like something without being that thing, but we won’t venture there today.

We introduced Fermat’s Little Theorem to show how true statements about prime numbers can be used to detect them. We took a statement that is true for all primes and then used that relationship to test if an integer is behaving the same way prime numbers do. Unlike our approach with Fermat’s Little Theorem, the Lucas-Lehmer test is definitive and does not have the same “walks like a duck, quacks like a duck, but isn’t a duck” problem that the Carmichael Numbers create. Lucas-Lehmer’s proof closes that gap in two steps. First, Fermat’s Little Theorem is true for all primes, while the Lucas-Lehmer test only applies to primes of a very specific form. Second, the Lucas-Lehmer test uses iteration alongside exponentiation and modular arithmetic to create a sequence of numbers that end in zero if and only if the candidate is prime. This iteration is also why the test can be slow and why GIMPS runs it last.

With the background that the Lucas-Lehmer Test uses iteration to definitively prove a number prime, we can finally write it down. The entire test fits in 10 lines of Python. This code isn’t good enough to set new records, but it’ll prove, and beat, Lucas!

def LucasLehmer(p):
    '''Returns True if the number 2^(p - 1) is prime, otherwise returns False'''
    s = 4
    M = (2**p) - 1
    for i in range(p - 2):
        s = (s**2 - 2) % M
    if s == 0:
        return True
    else:
        return False

p = 127
lucas_prime = (2 ** p) - 1
if LucasLehmer(p):
    print(lucas_prime, 'is prime!')
else:
    print(lucas_prime, 'is not prime.')

Running the code above gives us the output below.

170141183460469231731687303715884105727 is prime!

The number that took the record from Lucas, and from humans once and for all, was $M_{521} = 2^{521} – 1$. We can test that one too.

p = 521
M_p = (2 ** p) - 1
if LucasLehmer(p):
    print(M_p, 'is prime!')
else: 
    print(M_p, 'is not prime.')

6864797660130609714981900799081393217269435300143305409394463459185543183397656052122559640661454554977296311391480858037121987999716643812574028291115057151 is prime!

We’ve beaten Lucas! To be fair, we did not really find a number bigger than Lucas’ – we looked one up and checked it. Next, we show how to discover $M_{521}$ for ourselves and then find even bigger primes. We use the same ideas GIMPS does when deciding which numbers to test next.

How do we know which numbers to test?

The Lucas-Lehmer test we’re using is powerful and easy to write down, but like many good things, these advantages come with a trade-off. The trade-off is our test only works in very specific settings. This test can only determine if the number $M_p = 2^p – 1$ is prime. It is not guaranteed that $M_p$ is prime, but $p$ must be prime for $M_p$ to even stand a chance. In short, we need a small prime number, like 2, 3, 127, or 521, to try to make a bigger one.

The proof that $p$ must be prime is challenging but we can try a few values to see how this idea could work.

$p=2$ is prime and $2^2 – 1 = 4 – 1 = 3$ is also prime.

$p=3$ is prime and $2^3 – 1 = 8 – 1 = 7$ is also prime.

$p=4$ is not prime and $2^4 – 1 = 16 – 1 = 15 = 3 \times 5$ is also not prime.

Every time I focus on abstractions for too long I’m always a bit shocked when examples like this work out. For completeness, we’ll see that if $p$ is prime $M_p$ doesn’t have to be prime, so one last example.

$p=11$ is prime and $2^{11} – 1 = 2048 – 1 = 2047 = 23 \times 89$ is not prime.

The point is that before we find a big prime we need a lot of small ones. There are several ways to find small primes, but an easy method that is fast enough for us is the Sieve of Eratosthenes. This wonderful animation from Wikipedia shows how it works, but you don’t need to understand the details to run the code. 

def sieveOfEratosthenes(n):
    '''Returns a list of all prime numbers up to n. Fairly quick for n up to 1,000,000.'''
    candidates = {i:True for i in range(2, n)}
    for i in candidates:
        if candidates[i]:
            for j in range(2*i, n, i):
                candidates[j] = False
    primes = [i for i in candidates if candidates[i]]
    
    return primes

With this function we can generate lots of small primes. Let’s get all the primes up to 1,000 and beat Lucas without cheating.

N = 1000
primes = sieveOfEratosthenes(N)
print(f'Testing the {len(primes)} primes less than {N}. The largest prime less than {N} is {max(primes)}.')

mersenne_primes = []
for p in primes:
    if LucasLehmer(p):
        print(f'2^{p} - 1 = {2**p -1} is prime')
        mersenne_primes.append(p)

print(f'\nFound {len(mersenne_primes)} Mersenne Primes - {mersenne_primes}')
print(f'The largest is 2^{max(mersenne_primes)} - 1 and it has {len(str(2**max(mersenne_primes)-1))} digits.')

The block above will generate a list of small prime numbers, each less than 1,000, and test each Mersenne Prime $M_p$ using the Lucas-Lehmer Primality Test. The code prints out each Mersenne Prime as it finds it and then tells us how big the longest one is. Its output is below.

Testing the 168 primes less than 1000. The largest prime less than 1000 is 997.
2^3 - 1 = 7 is prime
2^5 - 1 = 31 is prime
2^7 - 1 = 127 is prime
2^13 - 1 = 8191 is prime
2^17 - 1 = 131071 is prime
2^19 - 1 = 524287 is prime
2^31 - 1 = 2147483647 is prime
2^61 - 1 = 2305843009213693951 is prime
2^89 - 1 = 618970019642690137449562111 is prime
2^107 - 1 = 162259276829213363391578010288127 is prime
2^127 - 1 = 170141183460469231731687303715884105727 is prime
2^521 - 1 = 6864797660130609714981900799081393217269435300143305409394463459185543183397656052122559640661454554977296311391480858037121987999716643812574028291115057151 is prime
2^607 - 1 = 531137992816767098689588206552468627329593117727031923199444138200403559860852242739162502265229285668889329486246501015346579337652707239409519978766587351943831270835393219031728127 is prime

Found 13 Mersenne Primes - [3, 5, 7, 13, 17, 19, 31, 61, 89, 107, 127, 521, 607]
The largest is 2^607 - 1 and it has 183 digits.

The largest prime we found is $M_{607}$, which is even bigger than the 1952 recorder holder $M_{527}$. That means your computer has caught up with more than 2,000 years of mathematicians. Your computer just beat every human who has ever tried to find a large prime number by hand and you know something the famous Lucas did not. Try turning up the $N$ from $1,000$ to $2,000$ and see if your computer can find a prime so big it can’t fit in a tweet – that is, find a prime with more than 280 digits.

Going further

My laptop was able to check all numbers $2^p -1 $ for $p < 10,000$ and found the following $p$ generate Mersenne Primes.

\[ 3, 5, 7, 13, 17, 19, 31, 61, 89, 107, 127, 521, 607, 1279, 2203, 2281, 3217, 4253, 4423, 9689, 9941 \]

That matches the OEIS list and Wikipedia page for Mersenne primes with $p < 10,000$. The biggest of these, $M_{9941} = 2^{9941} – 1$, held the record in 1963 and is mind-bogglingly big.

M_9941 = (2 ** 9941) - 1
print(f'M_9941 has {len(str(M_9941))} digits!')
print(M_9941, 'is prime.')

M_9941 has 2993 digits!
34608828249085121524296039576741331672262866890023854779048928344500622080983411446436437554415370753366448674763505018641470709332373970608376690404229265789647993709760358469552319045484910050304149809818540283507159683562232941968059762281334544739720849260904855192770626054911793590389060795981163838721432994278763633095377438194844866471124967685798888172212033000821469684464956146997194126921284336206463313859537577200462442029064681326087558257488470489384243989270236884978643063093004422939603370010546595386302009073043944482202559097406700597330570799507832963130938739885080198416258635194522913042562936679859587495721031173747796418895060701941717506001937152430032363631934265798516236047451209089864707430780362298307038193445486493756647991804258775574973833903315735082891029392359352758617185019942554834671861074548772439880729606244911940066680112823824095816458261761861746604034802056466823143718255492784779380991749580255263323326536457743894150848953969902818530057870876229329803338285735419228259022169602665532210834789602051686546011466737981306056247480055071718250333737502267307344178512950738594330684340802698228963986562732597175372087295649072830289749771358330867951508710859216743218522918811670637448496498549094430541277444079407989539857469452772132166580885754360477408842913327292948696897496141614919739845432835894324473601387609643750514699215032683744527071718684091832170948369396280061184593746143589068811190253101873595319156107319196071150598488070027088705842749605203063194191166922106176157609367241948160625989032127984748081075324382632093913796444665700601391278360323002267434295194325607280661260119378719405151497555187549252134264394645963853964913309697776533329401822158003182889278072368602128982710306618115118964131893657845400296860012420391376964670183983594954112484565597312460737798777092071706710824503707457220155015899591766244957768006802482976673920392995410164224776445671222149803657927708412925555542817045572430846389988129960519227313987291200902060882060733762075892299473666405897427035811786879875694315078654420055603469625309399653955932310466430039146465805452965014040019423897552675534768248624631951431493188170905972588780111850281190559073677771187432814088678674286302108275149258477101296451833651979717375170900505673645964696355331369819296000267389583289299126738345726980325998955997501176664201042888546085699446442834195232948787488410595750197438786353119204210855804692460582533832967771946911459901921324984968810021189968284941331573164056304725480868921823442538199590383852412786840833479611419970101792978355653650755329138298654246225346827207503606740745956958127383748717825918527473164970582095181312905519242710280573023145554793628499010509296055849712377978984921839997037415897674154830708629145484724536724572622450131479992681684310464449439022250504859250834761894788889552527898400988196200014868575640233136509145628127191354858275083907891469979019426224883789463551 is prime.

To catch up with the current record holder you would need a much more powerful computer and an improved version of our test. Currently, GIMPS takes about a month to test a single number and you can contribute directly to that project by running their code on your computer. You can learn more about the Greater Internet Mersenne Prime Search on their “How GIMPS Works” page.

Thank you for reading!