Now we come to the Python program that amazed me.

The problem with version 4 was that my choice of setting the search limit via a floating point value choked the program prematurely — Python long integers, it turns out, are as technically beyond Python floats as a spaceship is beyond a potato cannon.

I changed the code so that the limit is not the actual Mersenne prime value, but the exponent of the Mersenne function — in other words, for Mn = 2^n - 1, I’m limiting n, not Mn.

Here’s the revised code for mersenne5.py (in “more” tag hyperspace, of course):

mersenne5.py:

#!/usr/bin/env python

__metaclass__ = type
import time
from math import sqrt, log

def is_prime(n):
    """Determines if a positive number is prime or not."""
    if n == 2:
        return True
    if n % 2 == 0:
        return False
    max = int(sqrt(n))
    factor = 1
    while factor < max:
        factor += 2
        if n % factor == 0:
            return False
    return True

def lucas_lehmer_test(prime, mersenne_prime):
    """
    The actual Lucas-Lehmer test only works for odd exponents.

    http://en.wikipedia.org/wiki/Lucas-Lehmer_test_for_Mersenne_numbers
    """
    s = 4
    for x in range(prime-2):
        s = ((s*s) - 2) % mersenne_prime
    if s == 0:
        return True
    else:
        return False

def mersenne(limit):
    """
    Returns a list of Mersenne Primes up to a given limit.
    This new version uses the Lucas-Lehmer test:

    http://en.wikipedia.org/wiki/Lucas-Lehmer_test_for_Mersenne_numbers

    Because it only works for odd exponents, there is one Mersenne prime
    which would habitually be left off -- M(2) = 3.  Therefore, I'm starting
    with 3 as a given.  An alternative would be to add ten lines of code
    to catch a math problem that could probably be worked out by a third
    grader using pencil and paper.

    On the other hand, the Lucas-Lehmer test is so efficient I went beyond
    the maximum floating point limit and had to switch my definition for
    the limiting value.  Now, mersenne(limit) defines the limit of the base
    prime, not the Mersenne Prime.  Therefore mersenne(1000) would search
    values up to 2^1000.
    """
    cache = [3]
    for x in xrange(2,limit):
        if is_prime(x):
            candidate = (2 ** x) - 1
            if lucas_lehmer_test(x, candidate):
                cache.append(candidate)
    return cache

if __name__ == '__main__':
    print "Start:", time.ctime()
    print mersenne(1300)
    print "Finish:", time.ctime()

OK, now the test run. I was a bit paranoid at first, thinking this was entirely too good to be true, so I started with a limit of 1300 — good enough to find the 15th Mersenne Prime:

Start: Mon May 21 19:48:49 2007
[3, 7, 31,
127, 8191, 131071, 524287, 2147483647L,
2305843009213693951L, 618970019642690137449562111L,
162259276829213363391578010288127L,
170141183460469231731687303715884105727L,
686479766013060971498190079908139321726943530014330540939446
345918554318339765605212255964066145455497729631139148085803
7121987999716643812574028291115057151L,
531137992816767098689588206552468627329593117727031923199444
138200403559860852242739162502265229285668889329486246501015
346579337652707239409519978766587351943831270835393219031728
127L,
104079321946643990819252403273640855386152622472667048053191
123504036080596733602980122394417323241848424216139542810077
913835662483234649081399066056773207629241295093892203457731
833496615835504729594205476898112116936771475484788669625013
844382602917323488853111608285384165850282556046662248318909
188018470682222031405210266984354887329580288780508697361869
00714720710555703168729087L]
Finish: Mon May 21 19:48:55 2007

real    0m6.567s
user    0m6.489s
sys     0m0.040s

And there it is, in only 6.5 seconds. How high can this go?

I bumped the exponent limit up to 4500 — we’re talking reliable math on thousand-plus digit numbers, now. The result:

Start: Mon May 21 19:51:45 2007
[3, 7, 31,
127, 8191, 131071, 524287, 2147483647L,
2305843009213693951L, 618970019642690137449562111L,
162259276829213363391578010288127L,
170141183460469231731687303715884105727L,
686479766013060971498190079908139321726943530014330540939446
345918554318339765605212255964066145455497729631139148085803
7121987999716643812574028291115057151L,
531137992816767098689588206552468627329593117727031923199444
138200403559860852242739162502265229285668889329486246501015
346579337652707239409519978766587351943831270835393219031728
127L,
104079321946643990819252403273640855386152622472667048053191
123504036080596733602980122394417323241848424216139542810077
913835662483234649081399066056773207629241295093892203457731
833496615835504729594205476898112116936771475484788669625013
844382602917323488853111608285384165850282556046662248318909
188018470682222031405210266984354887329580288780508697361869
00714720710555703168729087L,
147597991521418023508489862273738173631206614533316977514777
121647857029787807894937740733704938928938274850753149648047
728126483876025919181446336533026954049696120111343015690239
609398909022625932693502528140961498349938822283144859860183
431853623092377264139020949023183644689960821079548296376309
423663094541083279376990539998245718632294472963641889062337
217172374210563644036821845964963294853869690587265048691443
463745750728044182367681351785209934866084717257940842231667
809767022401199028017047489448742692474210882353680848507250
224051945258754287534997655857267022963396257521263747789778
550155264652260998886991401354048380986568125041949768669777
1007L,
446087557183758429571151706402101809886208632412859901111991
219963404685792820473369112545269003989026153245931124316702
395758705693679364790903497461147071065254193353938124978226
307947312410798874869040070279328428810311754844108094878252
494866760969586998128982645877596028979171536962503068429617
331702184750324583009171832104916050157628886606372145501702
225925125224076829605427173573964812995250569412480720738476
855293681666712844831190877620606786663862190240118570736831
901886479225810414714078935386562497968178729127629594924411
960961386713946279899275006954917139758796061223803393537381
034666494402951052059047968693255388647930440925104186817009
640171764133172418132836351L,
259117086013202627776246767922441530941818887553125427303974
923161874019266586362086201209516800483406550695241733194177
441689509238807017410377709597512042313066624082916353517952
311186154862265604547691127595848775610568757931191017711408
826252153849035830401185072116424747461823031471398340229288
074545677907941037288235820705892351068433882986888616658650
280927692080339605869308790500409503709875902119018371991620
994002568935113136548829739112656797303241986517250116412703
509705427773477972349821676443446668383119322540099648994051
790241624056519054483690809616061625743042361721863339415852
426431208737266591962061753535748892894599629195183082621860
853400937932839420261866586142503251450773096274235376822938
649407127700846077124211823080804139298087057504713825264571
448379371125032081826126566649084251699453951887789613650248
405739378594599444335231188280123660406262468609212150349937
584782292237144339628858485938215738821232393687046160677362
909315071L,
190797007524439073807468042969529173669356994749940177394741
882673528979787005053706368049835514900244303495954950709725
762186311224148828811920216904542206960744666169364221195289
538436845390250168663932838805192055137154390912666527533007
309292687539092257043362517857366624699975402375462954490293
259233303137330643531556539739921926201438606439020075174723
029056838272505051571967594608350063404495977660656269020823
960825567012344189908927956646011998057988548630107637380993
519826582389781888135705408653045219655801758081251164080554
609057468028203308718724654081055323215860189611391296030471
108443146745671967766308925858547271507311563765171008318248
647110097614890313562856541784154881743146033909602737947385
055355960331855614540900081456378659068370317267696980001187
750995491090350108417050917991562167972281070161305972518044
872048331306383715094854938415738549894606070722584737978176
686422134354526989443028353644037187375385397838259511833166
416134323695660367676897722287918773420968982326089026150031
515424165462111337527431154890666327374921446276833564519776
797633875503548665093914556482031482248883127023777039667707
976559857333357013727342079099064400455741830654320379350833
236245819348824064783585692924881021978332974949906122664421
376034687815350484991L,
285542542228279613901563566102164008326164238644702889199247
456602284400390600653875954571505539843239754513915896150297
878399377056071435169747221107988791198200988477531339214282
772016059009904586686254989084815735422480409022344297588352
526004383890632616124076317387416881148592486188361873904175
783145696016919574390765598280188599035578448591077683677175
520434074287726578006266759615970759521327828555662781678385
691581844436444812511562428136742490459363212810180276096088
111401003377570363545725120924073646921576797146199387619296
560302680261790118132925012323046444438622308877924609373773
012481681672424493674474488537770155783006880852648161513067
144814790288366664062257274665275787127374649231096375001170
901890786263324619578795731425693805073056119677580338084333
381987500902968831935913095269821311141322393356490178488728
982288156282600813831296143663845945431144043753821542871277
745606447858564159213328443580206422714694913091762716447041
689678070096773590429808909616750452927258000843500344831628
297089902728649981994387647234574276263729694848304750917174
186181130688518792748622612293341368928056634384466646326572
476167275660839105650528975713899320211121495795311427946254
553305387067821067601768750977866100460014602138408448021225
053689054793742003095722096732954750721718115531871310231057
902608580607L]
Finish: Mon May 21 20:04:26 2007

real    12m40.597s
user    10m48.031s
sys     0m4.486s

Amazing. I was actually profoundly moved by this. In just 12 minutes, it found the first 20 Mersenne primes — the 20th being a 1,332 digit number.

How high can this go?

It’s starting to look like a very interesting question.

Scripting languages get a bad rap for poor performance; and yet, here, I’ve seen a rapid development of a module capable of doing insanely complex mathematical operations. Had I tried to write this in a “faster” language like C — well, I’d still probably be setting up Makefiles, honestly, or programming code for handling arbitrarily long integers. Python has it already in hand. I wasn’t expecting such a deep and persistent level of quality. My hat is off to Python!

Software examples on this page are licensed under the CC-GNU GPL.