# Lucas-Lehmer test

Lucas-Lehmer test
You are encouraged to solve this task according to the task description, using any language you may know.

Lucas-Lehmer Test: for p an odd prime, the Mersenne number 2p − 1 is prime if and only if 2p − 1 divides S(p − 1) where S(n + 1) = (S(n))2 − 2, and S(1) = 4.

The following programs calculate all Mersenne primes up to the implementation's maximum precision, or the 47th Mersenne prime. (Which ever comes first).

## Contents

with Ada.Text_Io; use Ada.Text_Io;with Ada.Integer_Text_Io; use Ada.Integer_Text_Io; procedure Lucas_Lehmer_Test is   type Ull is mod 2**64;   function Mersenne(Item : Integer) return Boolean is      S : Ull := 4;      MP : Ull := 2**Item - 1;   begin      if Item = 2 then         return True;      else         for I in 3..Item loop            S := (S * S - 2) mod MP;         end loop;         return S = 0;      end if;   end Mersenne;   Upper_Bound : constant Integer := 64;begin   Put_Line(" Mersenne primes:");   for P in 2..Upper_Bound loop      if Mersenne(P) then         Put(" M");         Put(Item => P, Width => 1);      end if;   end loop;end Lucas_Lehmer_Test;

Output:

Mersenne primes:
M2 M3 M5 M7 M13 M17 M19 M31


## Agena

Because of the very large numbers computed, the mapm binding is used to calculate with arbitrary precision.

readlib 'mapm'; mapm.xdigits(100); mersenne := proc(p::number) is   local s, m;   s := 4;   m := mapm.xnumber(2^p) - 1;   if p = 2 then      return true   else      for i from 3 to p do         s := (mapm.xnumber(s)^2 - 2) % m      od;      return mapm.xtoNumber(s) = 0   fiend; for i from 3 to 64 do   if mersenne(i) then      write('M' & i & ' ')   fiod;

produces:

M3 M5 M7 M13 M17 M19 M31

## ALGOL 68

Works with: ALGOL 68 version Revision 1 - no extensions to language used
Works with: ALGOL 68G version Any - tested with release 1.18.0-9h.tiny
PRAGMAT stack=1M precision=20000 PRAGMAT PROC is prime = ( INT p )BOOL:  IF p = 2 THEN TRUE  ELIF p <= 1 OR p MOD 2 = 0 THEN FALSE  ELSE    BOOL prime := TRUE;    FOR i FROM 3 BY 2 TO ENTIER sqrt(p)      WHILE prime := p MOD i /= 0 DO SKIP OD;    prime  FI; PROC is mersenne prime = ( INT p )BOOL:  IF p = 2 THEN TRUE  ELSE    LONG LONG INT m p :=  LONG LONG 2 ** p - 1, s := 4;    FROM 3 TO p DO      s := (s ** 2 - 2) MOD m p    OD;    s = 0  FI; test:(  INT upb prime = ( long long bits width - 1 ) OVER 2; # no unsigned #  INT upb count = 45; # find 45 mprimes if INT has enough bits #   printf(($" Finding Mersenne primes in M[2.."g(0)"]: "l$,upb prime));   INT count:=0;  FOR p FROM 2 TO upb prime WHILE    IF is prime(p) THEN      IF is mersenne prime(p) THEN        printf (($" M"g(0)$,p));        count +:= 1      FI    FI;    count <= upb count  DO SKIP OD)

Output:

Finding Mersenne primes in M[2..33252]:
M2 M3 M5 M7 M13 M17 M19 M31 M61 M89 M107 M127 M521 M607 M1279 M2203 M2281 M3217 M4253 M4423 M9689 M9941 M11213 M19937 M21701 M23209


## AWK

 # syntax: GAWK -f LUCAS-LEHMER_TEST.AWK# converted from PascalBEGIN {    printf("Mersenne primes:")    n = 1    for (exponent=2; exponent<=32; exponent++) {      s = (exponent == 2) ? 0 : 4      n = (n+1)*2-1      for (i=1; i<=exponent-2; i++) {        s = (s*s-2)%n      }      if (s == 0) {        printf(" M%s",exponent)      }    }    printf("\n")    exit(0)}

output:

Mersenne primes: M2 M3 M5 M7 M13 M17 M19


## BBC BASIC

Using its native arithmetic BBC BASIC can only test up to M23.

      *FLOAT 64      PRINT "Mersenne Primes:"      FOR p% = 2 TO 23        IF FNlucas_lehmer(p%) PRINT "M" ; p%      NEXT      END       DEF FNlucas_lehmer(p%)      LOCAL i%, mp, sn      IF p% = 2 THEN = TRUE      IF (p% AND 1) = 0 THEN = FALSE      mp = 2^p% - 1      sn = 4      FOR i% = 3 TO p%        sn = sn^2 - 2        sn -= (mp * INT(sn / mp))      NEXT      = (sn = 0)

Output:

Mersenne Primes:
M2
M3
M5
M7
M13
M17
M19


## Bracmat

Only exponents that are prime are tried. The primality test of these numbers uses a side effect of Bracmat's attempt at computing a root of a small enough number. ('small enough' meaning that the number must fit in a computer word, normally 32 or 64 bits.) To do that, Bracmat first creates a list of factors of the number and then takes the root of each factor. For example, to compute 54^2/3, Bracmat first creates the expression (2*3^3)^2/3 and then 2^2/3*3^(3*2/3), which becomes 2^2/3*9. If a number cannot be factorized, (either because it is prime or because it is to great to fit in a computer word) the root expression doesn't change much. For example, the expression 13^(13^-1) becomes 13^1/13, and this matches the pattern 13^%.

  ( clk$:?t0:?now & ( time = ( print = . put$ ( str              $( div$(!arg,1)                  ","                  (   div$(mod$(!arg*100,100),1):?arg                    & !arg:<10                    & 0                  |                   )                  !arg                  " "                )              )        )      & -1*!now+(clk$:?now):?SEC & print$!SEC      & print$(!now+-1*!t0) & put$"s: "    )  & 3:?exponent  &   whl    ' ( !exponent:~>12000      & (   !exponent^(!exponent^-1):!exponent^%          & 4:?s          & 2^!exponent+-1:?n          & 0:?i          &   whl            ' ( 1+!i:?i              & !exponent+-2:~<!i              & mod$(!s^2+-2.!n):?s ) & ( !s:0 & !time & out$(str$(M !exponent " is PRIME!")) | ) | ) & 1+!exponent:?exponent ) & done ); Output (after 4.5 hours): 0,00 0,00 s: M3 is PRIME! 0,00 0,00 s: M5 is PRIME! 0,00 0,00 s: M7 is PRIME! 0,00 0,00 s: M13 is PRIME! 0,00 0,00 s: M17 is PRIME! 0,00 0,01 s: M19 is PRIME! 0,00 0,01 s: M31 is PRIME! 0,00 0,01 s: M61 is PRIME! 0,01 0,02 s: M89 is PRIME! 0,01 0,03 s: M107 is PRIME! 0,00 0,04 s: M127 is PRIME! 0,50 0,54 s: M521 is PRIME! 0,29 0,84 s: M607 is PRIME! 6,81 7,65 s: M1279 is PRIME! 38,35 46,01 s: M2203 is PRIME! 6,32 52,33 s: M2281 is PRIME! 116,01 168,34 s: M3217 is PRIME! 293,09 461,44 s: M4253 is PRIME! 64,61 526,05 s: M4423 is PRIME! 8863,90 9389,95 s: M9689 is PRIME! 1101,12 10491,08 s: M9941 is PRIME! 5618,45 16109,53 s: M11213 is PRIME! ## C Works with: gcc version 4.1.2 20070925 (Red Hat 4.1.2-27) Works with: C99 Compiler options: gcc -std=c99 -lm Lucas-Lehmer_test.c -o Lucas-Lehmer_test #include <math.h>#include <stdio.h>#include <limits.h>#pragma precision=log10l(ULLONG_MAX)/2 typedef enum { FALSE=0, TRUE=1 } BOOL; BOOL is_prime( int p ){ if( p == 2 ) return TRUE; else if( p <= 1 || p % 2 == 0 ) return FALSE; else { BOOL prime = TRUE; const int to = sqrt(p); int i; for(i = 3; i <= to; i+=2) if (!(prime = p % i))break; return prime; }} BOOL is_mersenne_prime( int p ){ if( p == 2 ) return TRUE; else { const long long unsigned m_p = ( 1LLU << p ) - 1; long long unsigned s = 4; int i; for (i = 3; i <= p; i++){ s = (s * s - 2) % m_p; } return s == 0; }} int main(int argc, char **argv){ const int upb = log2l(ULLONG_MAX)/2; int p; printf(" Mersenne primes:\n"); for( p = 2; p <= upb; p += 1 ){ if( is_prime(p) && is_mersenne_prime(p) ){ printf (" M%u",p); } } printf("\n");} Output: Mersenne primes: M2 M3 M5 M7 M13 M17 M19 M31  ## C++ Library: GMP #include <iostream>#include <gmpxx.h> mpz_class pow2(mpz_class exp);bool is_mersenne_prime(mpz_class p); int main(){ mpz_class maxcount(45); mpz_class found(0); mpz_class check(0); for( mpz_nextprime(check.get_mpz_t(), check.get_mpz_t()); found < maxcount; mpz_nextprime(check.get_mpz_t(), check.get_mpz_t())) { //std::cout << "P" << check << " " << std::flush; if( is_mersenne_prime(check) ) { ++found; std::cout << "M" << check << " " << std::flush; } }} bool is_mersenne_prime(mpz_class p){ if( 2 == p ) return true; else { mpz_class div = pow2(p) - mpz_class(1); mpz_class s(4); mpz_class s(4); for( mpz_class i(3); i <= p; ++i ) { s = (s * s - mpz_class(2)) % div ; } return ( s == mpz_class(0) ); }} mpz_class pow2(mpz_class exp){ // Unfortunately, GMP doesn't have a left-shift method. // It also doesn't have a pow() equivalent that takes arbitrary-precision exponents. // So we have to do it the hard (and presumably slow) way. mpz_class ret(2); mpz_class ret(2); for(mpz_class i(1); i < exp; ++i) ret *= mpz_class(2); ret *= mpz_class(2); //std::cout << "pow2( " << exp << " ) = " << ret << std::endl; return ret; } Output: (Incomplete; It takes a long time.) M2 M3 M5 M7 M13 M17 M19 M31 M61 M89 M107 M127 M521 M607 M1279 M2203 M2281 M3217 M4253 M4423 M9689 M9941 M11213 M19937 M21701 M23209 M44497  ## C# Works with: Visual Studio version 2010 Works with: .NET Framework version 4.0 using System;using System.Collections.Generic;using System.Numerics;using System.Threading.Tasks; namespace LucasLehmerTestForRosettaCode{ public class LucasLehmerTest { static BigInteger ZERO = new BigInteger(0); static BigInteger ONE = new BigInteger(1); static BigInteger TWO = new BigInteger(2); static BigInteger FOUR = new BigInteger(4); private static bool isMersennePrime(int p) { if (p % 2 == 0) return (p == 2); else { for (int i = 3; i <= (int)Math.Sqrt(p); i += 2) if (p % i == 0) return false; //not prime BigInteger m_p = BigInteger.Pow(TWO, p) - ONE; BigInteger s = FOUR; for (int i = 3; i <= p; i++) s = (s * s - TWO) % m_p; return s == ZERO; } } public static int[] GetMersennePrimeNumbers(int upTo) { List<int> response = new List<int>(); Parallel.For(2, upTo, i => { if (isMersennePrime(i)) response.Add(i); }); response.Sort(); return response.ToArray(); } static void Main(string[] args) { int[] mersennePrimes = LucasLehmerTest.GetMersennePrimeNumbers(11213); foreach (int mp in mersennePrimes) Console.Write("M" + mp+" "); Console.ReadLine(); } }} Output: (Run only to 11213) M2 M3 M5 M7 M13 M17 M19 M31 M61 M89 M107 M127 M521 M607 M1279 M2203 M2281 M3217 M4253 M4423 M9689 M9941 M11213  ## Clojure (defn prime? [i] (cond (< i 4) (>= i 2) (zero? (rem i 2)) false :else (not-any? #(zero? (rem i %)) (range 3 (inc (Math/sqrt i)))))))) (defn mersenne? [p] (or (= p 2) (let [mp (dec (bit-shift-left 1 p))] (loop [n 3 s 4] (if (> n p) (zero? s) (recur (inc n) (rem (- (* s s) 2) mp))))))) (filter mersenne? (filter prime? (iterate inc 1))) Output: Infinite list of Mersenne primes: (2 3 5 7 13 17 19 31 61 89 107 127 521 607 1279 2203 2281 3217 4253...  ## D Translation of: Python import std.stdio, std.math, std.bigint; bool isPrime(in int p) pure nothrow { if (p < 2 || p % 2 == 0) return p == 2; foreach (immutable i; 3 .. cast(uint)real(p).sqrt + 1) if (p % i == 0) return false; return true;} bool isMersennePrime(in int p) pure /*nothrow*/ { if (!p.isPrime) return false; if (p == 2) return true; immutable mp = (1.BigInt << p) - 1; auto s = 4.BigInt; foreach (immutable _; 3 .. p + 1) s = (s ^^ 2 - 2) % mp; return s == 0;} void main() { foreach (immutable p; 2 .. 2_300) if (p.isMersennePrime) { write('M', p, ' '); stdout.flush; }} Output: M2 M3 M5 M7 M13 M17 M19 M31 M61 M89 M107 M127 M521 M607 M1279 M2203 M2281  With p up to 10_000 it prints: M2 M3 M5 M7 M13 M17 M19 M31 M61 M89 M107 M127 M521 M607 M1279 M2203 M2281 M3217 M4253 M4423 M9941  ## DWScript Using Integer type, which is 64bit, limits the search to M31. function IsMersennePrime(p : Integer) : Boolean;var i, s, m_p : Integer;begin if p=2 then Result:=True else begin m_p := (1 shl p)-1; s := 4; for i:=3 to p do s:=(s*s-2) mod m_p; Result:=(s=0); end;end; const upperBound = Round(Log2(High(Integer))/2); PrintLn('Finding Mersenne primes in M[2..' + IntToStr(upperBound) + ']: ');Print('M2'); var p : Integer;for p:=3 to upperBound step 2 do begin if IsMersennePrime(p) then Print(' M'+IntToStr(p));end;PrintLn(''); Output:  M2 M3 M5 M7 M13 M17 M19 M31 ## Erlang -module(mp).-export([main/0]). main() -> [ io:format("M~p ", [P]) || P <- lists:seq(2,700), (P == 2) orelse (s((1 bsl P) - 1, P-1) == 0) ]. s(MP,1) -> 4 rem MP;s(MP,N) -> X=s(MP,N-1), (X*X - 2) rem MP. In 3 seconds will print M2 M3 M5 M7 M13 M17 M19 M31 M61 M89 M107 M127 M521 M607  Testing larger numbers (i.e. 5000) is possible but will take few minutes. ## F# Simple arbitrary-precision version: let rec s mp n = if n = 1 then 4I % mp else ((s mp (n - 1)) ** 2 - 2I) % mp [ for p in 2..47 do if p = 2 || s ((1I <<< p) - 1I) (p - 1) = 0I then yield p ] Tail-recursive version: let IsMersennePrime exponent = if exponent <= 1 then failwith "Exponent must be >= 2" let prime = 2I ** exponent - 1I; let rec LucasLehmer i acc = match i with | x when x = exponent - 2 -> acc | x -> LucasLehmer (x + 1) ((acc*acc - 2I) % prime) LucasLehmer 0 4I = 0I  Version using library folding function (way shorter and faster than the above): let IsMersennePrime exponent = if exponent <= 1 then failwith "Exponent must be >= 2" let prime = 2I ** exponent - 1I; let LucasLehmer = [| 1 .. exponent-2 |] |> Array.fold (fun acc _ -> (acc*acc - 2I) % prime) 4I LucasLehmer = 0I  ## Fortran Works with: Fortran version 90 and later Only Mersenne number with prime exponent can be themselves prime but for the small numbers used in this example it was not worth the effort to include this check. As the size of the exponent increases this becomes more important. PROGRAM LUCAS_LEHMER IMPLICIT NONE INTEGER, PARAMETER :: i64 = SELECTED_INT_KIND(18) INTEGER(i64) :: s, n INTEGER :: i, exponent DO exponent = 2, 31 IF (exponent == 2) THEN s = 0 ELSE s = 4 END IF n = 2_i64**exponent - 1 DO i = 1, exponent-2 s = MOD(s*s - 2, n) END DO IF (s==0) WRITE(*,"(A,I0,A)") "M", exponent, " is PRIME" END DO END PROGRAM LUCAS_LEHMER ## GAP LucasLehmer := function(n) local i, m, s; if n = 2 then return true; elif not IsPrime(n) then return false; else m := 2^n - 1; s := 4; for i in [3 .. n] do s := RemInt(s*s, m) - 2; od; return s = 0; fi;end; Filtered([1 .. 2000], LucasLehmer);[2, 3, 5, 7, 13, 17, 19, 31, 61, 89, 107, 127, 521, 607, 1279] ## Go Processing the first list indicates that the test works. Processing the second shows it working on some larger numbers. package main import ( "fmt" "math/big") var primes = []uint{3, 5, 7, 11, 13, 17, 19, 23, 29, 31, 37, 41, 43, 47, 53, 59, 61, 67, 71, 73, 79, 83, 89, 97, 101, 103, 107, 109, 113, 127} var mersennes = []uint{521, 607, 1279, 2203, 2281, 3217, 4253, 4423, 9689, 9941, 11213, 19937, 21701, 23209, 44497, 86243, 110503, 132049, 216091, 756839, 859433, 1257787, 1398269, 2976221, 3021377, 6972593, 13466917, 20996011, 24036583} func main() { llTest(primes) fmt.Println() llTest(mersennes)} func llTest(ps []uint) { var s, m big.Int one := big.NewInt(1) two := big.NewInt(2) for _, p := range ps { m.Sub(m.Lsh(one, p), one) s.SetInt64(4) for i := uint(2); i < p; i++ { s.Mod(s.Sub(s.Mul(&s, &s), two), &m) } if s.BitLen() == 0 { fmt.Printf("M%d ", p) } }} Output: M3 M5 M7 M13 M17 M19 M31 M61 M89 M107 M127 M521 M607 M1279 M2203 M2281 M3217 M4253 M4423 M9689 M9941 M11213 M19937...  ## Haskell Works with: GHCi version 6.8.2 Works with: GHC version 6.8.2 module Main where main = printMersennes$ take 45 $filter lucasLehmer$ sieve [2..] s mp 1 = 4 mod mps mp n = ((s mp $n-1)^2-2) mod mp lucasLehmer 2 = TruelucasLehmer p = s (2^p-1) (p-1) == 0 printMersennes = mapM_ (\x -> putStrLn$ "M" ++ show x)

It is pointed out on the Sieve of Eratosthenes page that the following "sieve" is inefficient. Nonetheless it takes very little time compared to the Lucas-Lehmer test itself.

sieve (p:xs) = p : sieve [x | x <- xs, x mod p > 0]

It takes about 30 minutes to get up to:

M2 M3 M5 M7 M13 M17 M19 M31 M61 M89 M107 M127 M521 M607 M1279 M2203 M2281 M3217 M4253 M4423 M9689 M9941 M11213


## HicEst

s = 0DO exponent = 2, 31  IF(exponent > 2) s = 4  n = 2^exponent - 1  DO i = 1, exponent-2    s = MOD(s*s - 2, n)  ENDDO  IF(s == 0) WRITE(Messagebox) 'M', exponent, ' is prime;', nENDDO END

## Java

We use arbitrary-precision integers in order to be able to test any arbitrary prime.

import java.math.BigInteger;public class Mersenne{     public static boolean isPrime(int p) {        if (p == 2)            return true;        else if (p <= 1 || p % 2 == 0)            return false;        else {            int to = (int)Math.sqrt(p);            for (int i = 3; i <= to; i += 2)                if (p % i == 0)                    return false;            return true;        }    }     public static boolean isMersennePrime(int p) {        if (p == 2)            return true;        else {            BigInteger m_p = BigInteger.ONE.shiftLeft(p).subtract(BigInteger.ONE);            BigInteger s = BigInteger.valueOf(4);            for (int i = 3; i <= p; i++)                s = s.multiply(s).subtract(BigInteger.valueOf(2)).mod(m_p);            return s.equals(BigInteger.ZERO);        }    }     // an arbitrary upper bound can be given as an argument    public static void main(String[] args) {        int upb;        if (args.length == 0)            upb = 500;        else            upb = Integer.parseInt(args[0]);         System.out.print(" Finding Mersenne primes in M[2.." + upb + "]:\nM2 ");        for (int p = 3; p <= upb; p += 2)            if (isPrime(p) && isMersennePrime(p))                System.out.print(" M" + p);        System.out.println();    }}

Finding Mersenne primes in M[2..2147483647]:
M2 M3 M5 M7 M13 M17 M19 M31 M61 M89 M107 M127 M521 M607 M1279 M2203 M2281 M3217 M4253 M4423 M9689 M9941 M11213


## Mathematica

This version is very speedy and is bounded.

Select[Table[M = 2^p - 1;   For[i = 1; s = 4, i <= p - 2, i++, s = Mod[s^2 - 2, M]];   If[s == 0, "M" <> ToString@p, p], {p,    Prime /@ Range[300]}], StringQ] => {M3, M5, M7, M13, M17, M19, M31, M61, M89, M107, M127, M521, M607, M1279}

This version is unbounded (and timed):

t = SessionTime[];For[p = 2, True, p = NextPrime[p], M = 2^p - 1;  For[i = 1; s = 4, i <= p - 2, i++, s = Mod[s^2 - 2, M]];  If[s == 0, Print["M" <> ToString@p]]](SessionTime[] - t) {Seconds, Minutes/60, Hours/3600, Days/86400}

I'll see what this gets.

## MATLAB

MATLAB suffers from a lack of an arbitrary precision math (bignums) library. It also doesn't have great support for 64-bit integer arithmetic...or at least MATLAB 2007 doesn't. So, the best precision we have is doubles; therefore, this script can only find up to M19 and no greater.

function [mNumber,mersennesPrime] = mersennePrimes()     function isPrime = lucasLehmerTest(thePrime)         llResidue = 4;        mersennesPrime = (2^thePrime)-1;         for i = ( 1:thePrime-2 )            llResidue = mod( ((llResidue^2) - 2),mersennesPrime );                end                  isPrime = (llResidue == 0);     end     %Because IEEE764 Double is the highest precision number we can    %represent in MATLAB, the highest Mersenne Number we can test is 2^52.    %In addition, because we have this cap, we can only test up to the    %number 30 for Mersenne Primeness. When we input 31 into the    %Lucas-Lehmer test, during the computation of the residue, the    %algorithm multiplies two numbers together the result of which is    %greater than 2^53. Because we require every digit to be significant,    %this leads to an error. The Lucas-Lehmer test should say that M31 is a    %Mersenne Prime, but because of the rounding error in calculating the    %residues caused by floating-point arithmetic, it does not. So M30 is    %the largest number we test.     mNumber = (3:30);     [isPrime] = arrayfun(@lucasLehmerTest,mNumber);     mNumber = [2 mNumber(isPrime)];    mersennesPrime = (2.^mNumber) - 1; end

Ouput:

[mNumber,mersennesPrime] = mersennePrimes mNumber =      2     3     5     7    13    17    19  mersennesPrime =            3           7          31         127        8191      131071      524287

## Perl 6

multi is_prime(2) { True }multi is_prime(Int $p) {$p %% none(2,3,5,7...^ * > sqrt($p)) } multi is_mersenne_prime(2) { True }multi is_mersenne_prime(Int$p) {    my $m_p = 2 **$p - 1;    my $s = 4;$s = ($s ** 2 - 2) %$m_p for 3 .. $p;$s == 0;} say "M$_" if is_prime($_) and is_mersenne_prime($_) for 2..*; Output: M2 M3 M5 M7 M13 M17 M19 M31 M61 M89 M107 M127 M521 M607 M1279 ^C  ## PicoLisp (de prime? (N) (or (= N 2) (and (> N 1) (bit? 1 N) (for (D 3 T (+ D 2)) (T (> D (sqrt N)) T) (T (=0 (% N D)) NIL) ) ) ) ) (de mersenne? (P) (or (= P 2) (let (MP (dec (>> (- P) 1)) S 4) (do (- P 2) (setq S (% (- (* S S) 2) MP)) ) (=0 S) ) ) ) Output: : (for N 10000 (and (prime? N) (mersenne? N) (println N)) ) 2 3 5 7 13 17 19 31 61 89 107 127 521 607 1279 2203 2281 3217 4253 4423 9689 9941 ## Pop11 Checking large numbers takes a lot of time so we limit p to be smaller than 1000. define Lucas_Lehmer_Test(p); lvars mp = 2**p - 1, sn = 4, i; for i from 2 to p - 1 do (sn*sn - 2) rem mp -> sn; endfor; sn = 0;enddefine; lvars p = 3;printf('M2', '%p\n');while p < 1000 do if Lucas_Lehmer_Test(p) then printf('M', '%p'); printf(p, '%p\n'); endif; p + 2 -> p;endwhile; The output (obtained in few seconds) is: M2M3M5M7M13M17M19M31M61M89M107M127M521M607 ## PureBasic PureBasic has no large integer support. Calculations are limited to the range of a signed quad integer type. Procedure Lucas_Lehmer_Test(p) Protected mp.q = (1 << p) - 1, sn.q = 4, i For i = 3 To p sn = (sn * sn - 2) % mp Next If sn = 0 ProcedureReturn #True EndIf ProcedureReturn #FalseEndProcedure #upperBound = SizeOf(Quad) * 8 - 1 ;equivalent to significant bits in a signed quad integerIf OpenConsole() Define p = 3 PrintN("M2") While p <= #upperBound If Lucas_Lehmer_Test(p) PrintN("M" + Str(p)) EndIf p + 2 Wend Print(#CRLF$ + #CRLF$+ "Press ENTER to exit"): Input() CloseConsole()EndIf Sample output: M2 M3 M5 M7 M13 M17 M19 M31 ## Python from sys import stdoutfrom math import sqrt, log def is_prime ( p ): if p == 2: return True # Lucas-Lehmer test only works on odd primes elif p <= 1 or p % 2 == 0: return False else: for i in range(3, int(sqrt(p))+1, 2 ): if p % i == 0: return False return True def is_mersenne_prime ( p ): if p == 2: return True else: m_p = ( 1 << p ) - 1 s = 4 for i in range(3, p+1): s = (s ** 2 - 2) % m_p return s == 0 precision = 20000 # maximum requested number of decimal places of 2 ** MP-1 #long_bits_width = precision * log(10, 2)upb_prime = int( long_bits_width - 1 ) / 2 # no unsigned #upb_count = 45 # find 45 mprimes if int was given enough bits # print (" Finding Mersenne primes in M[2..%d]:"%upb_prime) count=0for p in range(2, upb_prime+1): if is_prime(p) and is_mersenne_prime(p): print("M%d"%p), stdout.flush() count += 1 if count >= upb_count: breakprint Output:  Finding Mersenne primes in M[2..33218]: M2 M3 M5 M7 M13 M17 M19 M31 M61 M89 M107 M127 M521 M607 M1279 M2203 M2281 M3217 M4253 M4423 M9689 M9941 M11213 M19937 M21701 M23209 ###  Faster loop without division def isqrt(n): if n < 0: raise ValueError elif n < 2: return n else: a = 1 << ((1 + n.bit_length()) >> 1) while True: b = (a + n // a) >> 1 if b >= a: return a a = b def isprime(n): if n < 5: return n == 2 or n == 3 elif n%2 == 0: return False else: r = isqrt(n) k = 3 while k <= r: if n%k == 0: return False k += 2 return True def lucas_lehmer_fast(n): if n == 2: return True elif not isprime(n): return False else: m = 2**n - 1 s = 4 for i in range(2, n): sqr = s*s r = (sqr & m) + (sqr >> n) if r >= m: r -= m s = r - 2 return s == 0 ## Racket  #lang racket(require math) (define (mersenne-prime? p) (divides? (- (expt 2 p) 1) (S (- p 1)))) (define (S n) (if (= n 1) 4 (- (sqr (S (- n 1))) 2))) (define (loop p) (when (mersenne-prime? p) (displayln p)) (loop (next-prime p))) (loop 3)  ## REXX REXX won't have a problem with the large number of digits involved, but since it's an interpreted language, such massive number crunching is not conducive to find large primes. /*REXX program to use Lucas-Lehmer primality test for prime powers of 2.*/parse arg limit . /*get the argument (if any). */if limit=='' then limit=1000 /*if no argument, assume 1000. */list= /*placeholder for the results. */ do j=1 by 2 to limit /*only process so much... */ /*there're only so many hours in a day...*/ power=j if j==1 then power=2 /*special case for the even prime*/ if \isprime(power ) then iterate /*if not prime, then ignore it. */ list=list Lucas_Lehmer2(power) /*add to list (or maybe not). */ end /*j*/ list=space(list) /*remove extraneous blanks. */say; say center('list',60-3,"═") /*show a fancy-dancy header/title*/say do k=1 for words(list) /*show entries in list, 1/line. */ say right(word(list,k),30) /*and right-justify 'em. */ end /*k*/exit /*stick a fork in it, we're done.*//*──────────────────────────────────LUCAS_LEHMER2 subroutine────────────*/Lucas_Lehmer2: procedure; parse arg ? /*Lucas-Lehmer test on 2**? - 1 */numeric form /*ensure correct scientific form.*/ if ?==2 then s=0 else s=4 /* DIGITs of 1 million could be used, but that really */ /* slows up the whole works. So, we start with the */ /* default of 9 digits, find the 10's exponent in */ /* the product (2**?), double it, and then add 6. */ /* 2 is all that's needed, but 6 is a lot safer. */ /* The doubling is for the square of S (s*s, below).*/ q=2**? /*compute a power of 2, using only 9 decimal digits. */ if pos('E',q)\==0 then do /*is it in exponential notation? */ parse var q 'E' tenpow numeric digits tenpow*2+6 end else numeric digits digits()*2+6 /* 9*2 + 6 */ q=2**?-1 /*now, compute the real McCoy. */ do ?-2 /*apply, rinse, repeat ... */ s=(s*s-2) // q /*modulus in REXX is: // */ end if s\==0 then return '' /*return nuttin' if not prime. */say "Mersenne prime via Lucas Lehmer test:", right('M'?,10) right("["digits() "digits]",25)return 'M'? /*return modified (prime) number.*//*──────────────────────────────────ISPRIME subroutine──────────────────*/isprime: procedure; parse arg xif x<17 then return wordpos(x,'2 3 5 7 11 13')\==0 /*test special cases*/if x//2==0 then return 0 /*is it even? Not prime.*/if x//3==0 then return 0 /*divisible by three? Not prime.*/if right(x,1)==5 then return 0 /*divisible by three? Not prime.*/if x//7==0 then return 0 /*divisible by seven? Not prime.*/ do j=11 by 6 /*ensures J isn't divisible by 3.*/ if x// j ==0 then return 0 if x//(j+2)==0 then return 0 if j*j>x then return 1 /*past the sqrt of X? It's prime*/ end /*j*/ output when the following is used for input: 10000 ══════════════════════════list═══════════════════════════ M2 M3 M5 M7 M13 M17 M19 M31 M61 M89 M107 M127 M521 M607 M1279 M2203 M2281 M3217 M4253 M4423 M9689 M9941  On IBM VM and MVS Rexx programs can be compiled. The above program uses (elapsed tim) interpreted 100 seconds compiled 1.5 seconds (for limit=2000) 20 seconds --[[User:Walterpachl|Walterpachl]] 09:55, 12 July 2012 (UTC)  ## RPL  %%HP: T(3)A(R)F(.); ; ASCII transfer header \<< DUP LN DUP \pi * 4 SWAP / 1 + UNROT / * IP 2 { 2 } ROT 2 SWAP ; input n; n := Int(n/ln(n)*(1 + 4/(pi*ln(n)))), p:=2; (n ~ number of primes less then n, pi used here only as a convenience), 2 is assumed to be the 1st elemente in the list START SWAP NEXTPRIME DUP UNROT DUP 2 SWAP ^ 1 - 4 PICK3 2 - 1 SWAP ; for i := 2 to n, p := nextprime; s := 4; m := 2^p - 1; START SQ 2 - OVER MOD ; for j := 1 to p - 2; s := s^2 mod m; NEXT NIP NOT { + } { DROP } IFTE ; next j; if s = 0 then add p to the list else discard p; NEXT NIP ; next i; \>>  Output: Outputs for arguments 130, 607 and 2281, respectively: { 2 3 5 7 13 17 19 31 61 89 107 127 } { 2 3 5 7 13 17 19 31 61 89 107 127 521 607 } { 2 3 5 7 13 17 19 31 61 89 107 127 521 607 1279 2203 2281 } These take respectively 1m 22s on the real HP 50g, 4m 29s and 10h 29m 23s on the emulator (Debug4 running on PC under WinXP, Intel(R) Core(TM) Duo CPU T2350 @ 1.86GHz).  ## Ruby def is_prime ( p ) return true if p == 2 return false if p <= 1 || p.even? (3 .. Math.sqrt(p)).step(2) do |i| return false if p % i == 0 end trueend def is_mersenne_prime ( p ) return true if p == 2 m_p = ( 1 << p ) - 1 s = 4 (p-2).times { s = (s ** 2 - 2) % m_p } s == 0end precision = 20000 # maximum requested number of decimal places of 2 ** MP-1 #long_bits_width = precision / Math.log(2) * Math.log(10)upb_prime = (long_bits_width - 1).to_i / 2 # no unsigned #upb_count = 45 # find 45 mprimes if int was given enough bits # puts " Finding Mersenne primes in M[2..%d]:" % upb_prime count = 0for p in 2..upb_prime if is_prime(p) && is_mersenne_prime(p) print "M%d " % p count += 1 end break if count >= upb_countendputs Output:  Finding Mersenne primes in M[2..33218]: M2 M3 M5 M7 M13 M17 M19 M31 M61 M89 M107 M127 M521 M607 M1279 M2203 M2281 M3217 M4253 M4423 M9689 M9941 M11213 M19937 M21701 M23209 ## Scala Works with: Scala version 2.10.2 Library: Scala In accordance with definition of Mersenne primes it could only be a Mersenne number with prime exponent. object LLT extends App { import Stream._ def primeSieve(s: Stream[Int]): Stream[Int] = s.head #:: primeSieve(s.tail filter { _ % s.head != 0 }) val primes = primeSieve(from(2)) def mersenne(p: Int): BigInt = (BigInt(2) pow p) - 1 def s(mp: BigInt, p: Int): BigInt = { if (p == 1) 4 else ((s(mp, p - 1) pow 2) - 2) % mp } val upbPrime = 9999 println(s"Finding Mersenne primes in M[2..$upbPrime]")  ((primes takeWhile (_ <= upbPrime)).par map { p => (p, mersenne(p)) }    map { p => if (p._1 == 2) (p, 0) else (p, s(p._2, p._1 - 1)) } filter { _._2 == 0 })    .foreach { p =>      println(s"prime M${(p._1)._1}: " + { if ((p._1)._1 < 200) (p._1)._2 else s"(${(p._1)._2.toString.size} digits)" })    }  println("That's All Folks!")}
Output:
After approx 20 minutes (2.10 GHz dual core)
Finding Mersenne primes in M[2..9999]
prime M2: 3
prime M3: 7
prime M5: 31
prime M7: 127
prime M13: 8191
prime M17: 131071
prime M19: 524287
prime M31: 2147483647
prime M61: 2305843009213693951
prime M89: 618970019642690137449562111
prime M107: 162259276829213363391578010288127
prime M127: 170141183460469231731687303715884105727
prime M521: (157 digits)
prime M607: (183 digits)
prime M1279: (386 digits)
prime M2203: (664 digits)
prime M2281: (687 digits)
prime M3217: (969 digits)
prime M4253: (1281 digits)
prime M4423: (1332 digits)
prime M9689: (2917 digits)
prime M9941: (2993 digits)
That's All Folks!

## Scheme

;;;The heart of the algorithm(define (S n)  (let ((m (- (expt 2 n) 1)))    (let loop ((c (- n 2)) (a 4))      (if (zero? c)          a          (loop (- c 1) (remainder (- (* a a) 2) m)))))) (define (mersenne-prime? n)  (if (= n 2)    #t    (zero? (S n)))) ;;;Trivial unoptimized implementation for the base primes(define (next-prime x)  (if (prime? (+ x 1))      (+ x 1)      (next-prime (+ x 1)))) (define (prime? x)  (let loop ((c 2))    (cond ((>= c x) #t)          ((zero? (remainder x c)) #f)          (else (loop (+ c 1)))))) ;;Main loop(let loop ((i 45) (p 2))  (if (not (zero? i))      (if (mersenne-prime? p)          (begin            (display "M") (display p) (display " ")            (loop (- i 1) (next-prime p)))          (loop i (next-prime p)))))
M2 M3 M5 M7 M13...


## Seed7

To get maximum speed the program should be compiled with -O2.

$include "seed7_05.s7i"; include "bigint.s7i"; const func boolean: isPrime (in integer: number) is func result var boolean: prime is FALSE; local var integer: upTo is 0; var integer: testNum is 3; begin if number = 2 then prime := TRUE; elsif number rem 2 = 0 or number <= 1 then prime := FALSE; else upTo := sqrt(number); while number rem testNum <> 0 and testNum <= upTo do testNum +:= 2; end while; prime := testNum > upTo; end if; end func; const func boolean: lucasLehmerTest (in integer: p) is func result var boolean: prime is TRUE; local var bigInteger: m_p is 0_; var bigInteger: s is 4_; var integer: i is 0; begin if p <> 2 then m_p := 2_ ** p - 1_; for i range 2 to pred(p) do s := (s ** 2 - 2_) rem m_p; end for; prime := s = 0_; end if; end func; const proc: main is func local var integer: p is 2; begin writeln(" Mersenne primes:"); while p <= 3217 do if isPrime(p) and lucasLehmerTest(p) then write(" M" <& p); flush(OUT); end if; incr(p); end while; writeln; end func; Original source: lucasLehmerTest, isPrime Output:  Mersenne primes: M2 M5 M7 M13 M17 M19 M31 M61 M89 M107 M127 M521 M607 M1279 M2203 M2281 M3217  ## Tcl Translation of: Pop11 proc main argv { set n 0 set t [clock seconds] show_mersenne 2 [incr n] t for {set p 3} {$p <= [lindex $argv 0]} {incr p 2} { if {![prime$p]} continue        if {[LucasLehmer $p]} { show_mersenne$p [incr n] t        }    }}proc show_mersenne {p n timevar} {    upvar 1 $timevar time set now [clock seconds] puts [format "%2d: %5ds M%s"$n [expr {$now -$time}] $p] set time$now}proc prime i {   if {$i in {2 3}} {return 1} prime0$i [expr {int(sqrt($i))}]}proc prime0 {i div} { expr {!($i % $div)? 0:$div <= 2? 1: [prime0 $i [incr div -1]]}}proc LucasLehmer p { set mp [expr {2**$p-1}]    set s  4    for {set i 2} {$i <$p} {incr i} {        set s [expr {($s**2 - 2) %$mp}]    }    expr {\$s == 0}} main 33218

Sample output. The program was still running, but as the next Mersenne prime is 19937 there will be a long wait until the program finds it.

 1:     0s  M2
2:     0s  M3
3:     0s  M5
4:     0s  M7
5:     0s  M13
6:     0s  M17
7:     0s  M19
8:     0s  M31
9:     0s  M61
10:     0s  M89
11:     0s  M107
12:     0s  M127
13:     1s  M521
14:     0s  M607
15:     4s  M1279
16:    21s  M2203
17:     4s  M2281
18:    69s  M3217
19:   180s  M4253
20:    39s  M4423
21:  5543s  M9689
22:   655s  M9941
23:  3546s  M11213