Lucas-Lehmer test: Difference between revisions
m →{{header|C++}}: Adjust commentary text |
m C/GMP, add libheader for GMP. |
||
Line 343:
Takes an optional argument to test up to the given value.
{{libheader|GMP}}
<lang C>#include <stdio.h>
#include <stdlib.h>
|
Revision as of 20:27, 21 November 2014
You are encouraged to solve this task according to the task description, using any language you may know.
Lucas-Lehmer Test: for an odd prime, the Mersenne number is prime if and only if divides where , and .
The following programs calculate all Mersenne primes up to the implementation's maximum precision, or the 47th Mersenne prime. (Which ever comes first).
360 Assembly
For maximum compatibility, this program uses only the basic instruction set. <lang 360asm>* Lucas-Lehmer test LUCASLEH CSECT
USING LUCASLEH,R12
SAVEARA B STM-SAVEARA(R15)
DC 17F'0' DC CL8'LUCASLEH'
STM STM R14,R12,12(R13) save calling context
ST R13,4(R15) ST R15,8(R13) LR R12,R15 set addessability
- ---- CODE
LA R2,2 R2=2 LA R11,0 R11:N' BCTR R11,0 N':=X'FFFFFFFF' LA R10,1 R10:N N=1 LA R4,1 R4:IEXP LA R6,1 step LH R7,IEXPMAX R7:IEXPMAX limit
LOOPE BXH R4,R6,ENDLOOPE do iexp=2 to iexpmax
SR R3,R3 R3:S S=0 CR R4,R2 if iexp=2 then S=0 BE OKS LA R3,4 else S=4
OKS EQU *
SLDA R10,1 n=(n+1)*2-1 LA R5,0 I LA R8,1 step LR R9,R4 IEXP SR R9,R2 IEXP-2 limit
LOOPI BXH R5,R8,ENDLOOPI do i=1 to iexp-2
- ---- compute s=(s*s-2) MOD n
SR R14,R14 R14=0 LR R15,R3 R15=S MR R14,R3 R{14-15}=S*S SLR R15,R2 R15=R15-2=S*S-2 BNM *+6 skip next if no borrow BCTR R14,0 perform borrow DR R14,R10 R10=N LR R3,R14 R14=MOD B LOOPI
ENDLOOPI EQU *
LTR R3,R3 BNZ NOPRT if s<>0 then no print CVD R4,P store to packed P UNPK Z,P Z=P MVC C,Z C=Z OI C+L'C-1,X'F0' zap sign MVC WTOBUF(4),C+12 MVI WTOBUF,C'M' WTO MF=(E,WTOMSG)
NOPRT EQU *
B LOOPE
ENDLOOPE EQU *
- ---- END CODE
RETURN EQU *
LM R14,R12,12(R13) XR R15,R15 BR R14
- ---- DATA
IEXPMAX DC H'31' I DS H IEXP DS H S DS F N DS F P DS PL8 packed Z DS ZL16 zoned C DS CL16 character WTOMSG DS 0F
DC H'80',XL2'0000'
WTOBUF DC 80C' '
LTORG YREGS END LUCASLEH</lang>
- Output:
M002 M003 M005 M007 M013 M017 M019 M031
Ada
<lang ada>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;</lang>
- 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. <lang agena>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 fi
end;
for i from 3 to 64 do
if mersenne(i) then write('M' & i & ' ') fi
od;</lang> produces: <lang agena>M3 M5 M7 M13 M17 M19 M31</lang>
ALGOL 68
<lang algol68>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
)</lang>
- 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
See also: http://www.xs4all.nl/~jmvdveer/mersenne.a68.html
AWK
<lang AWK>
- syntax: GAWK -f LUCAS-LEHMER_TEST.AWK
- converted from Pascal
BEGIN {
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)
} </lang>
- Output:
Mersenne primes: M2 M3 M5 M7 M13 M17 M19
BBC BASIC
Using its native arithmetic BBC BASIC can only test up to M23. <lang bbcbasic> *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)</lang>
- Output:
Mersenne Primes: M2 M3 M5 M7 M13 M17 M19
Bracmat
Only exponent
s 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^%
.
<lang bracmat> ( 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 );</lang>
- 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
GMP
This uses some pre-tests to show how we can skip some numbers with relatively inexpensive methods. This also does a simple optimization of the modulus. It takes about 30 seconds to get to M11213. This is substantially faster than many of the other solutions, though certainly not comparable to dedicated programs such as Prime95.
Takes an optional argument to test up to the given value.
<lang C>#include <stdio.h>
- include <stdlib.h>
- include <limits.h>
- include <gmp.h>
int lucas_lehmer(unsigned long p) {
mpz_t V, mp, t; unsigned long k, tlim; int res;
if (p == 2) return 1; if (!(p&1)) return 0;
mpz_init_set_ui(t, p); if (!mpz_probab_prime_p(t, 25)) /* if p is composite, 2^p-1 is not prime */ { mpz_clear(t); return 0; }
if (p < 23) /* trust the PRP test for these values */ { mpz_clear(t); return (p != 11); }
mpz_init(mp); mpz_setbit(mp, p); mpz_sub_ui(mp, mp, 1);
/* If p=3 mod 4 and p,2p+1 both prime, then 2p+1 | 2^p-1. Cheap test. */ if (p > 3 && p % 4 == 3) { mpz_mul_ui(t, t, 2); mpz_add_ui(t, t, 1); if (mpz_probab_prime_p(t,25) && mpz_divisible_p(mp, t)) { mpz_clear(mp); mpz_clear(t); return 0; } }
/* Do a little trial division first. Saves quite a bit of time. */ tlim = p/2; if (tlim > (ULONG_MAX/(2*p))) tlim = ULONG_MAX/(2*p); for (k = 1; k < tlim; k++) { unsigned long q = 2*p*k+1; /* factor must be 1 or 7 mod 8 and a prime */ if ( (q%8==1 || q%8==7) && q % 3 && q % 5 && q % 7 && mpz_divisible_ui_p(mp, q) ) { mpz_clear(mp); mpz_clear(t); return 0; } }
mpz_init_set_ui(V, 4); for (k = 3; k <= p; k++) { mpz_mul(V, V, V); mpz_sub_ui(V, V, 2); /* mpz_mod(V, V, mp) but more efficiently done given mod 2^p-1 */ if (mpz_sgn(V) < 0) mpz_add(V, V, mp); /* while (n > mp) { n = (n >> p) + (n & mp) } if (n==mp) n=0 */ /* but in this case we can have at most one loop plus a carry */ mpz_tdiv_r_2exp(t, V, p); mpz_tdiv_q_2exp(V, V, p); mpz_add(V, V, t); while (mpz_cmp(V, mp) >= 0) mpz_sub(V, V, mp); } res = !mpz_sgn(V); mpz_clear(t); mpz_clear(mp); mpz_clear(V); return res;
}
int main(int argc, char* argv[]) {
unsigned long i, n = 43112609; if (argc >= 2) n = strtoul(argv[1], 0, 10); for (i = 1; i <= n; i++) { if (lucas_lehmer(i)) { printf("M%lu ", i); fflush(stdout); } } printf("\n"); return 0;
}</lang>
- Output:
(partial output after 50 minutes)
M2 M3 M5 M7 M13 M17 M19 M31 M61 M89 M107 M127 M521 M607 M1279 M2203 M2281 M3217 M4253 M4423 M9689 M9941 M11213 M11213 M19937 M21701 M23209 M44497
Small inputs with native types
Compiler options: gcc -std=c99 -lm Lucas-Lehmer_test.c -o Lucas-Lehmer_test
<lang c>#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");
}</lang>
- Output:
Mersenne primes: M2 M3 M5 M7 M13 M17 M19 M31
C++
Straightforward method.
<lang cpp>#include <iostream>
- include <gmpxx.h>
static bool is_mersenne_prime(mpz_class p) {
if( 2 == p ) return true; else { mpz_class s(4); mpz_class div( (mpz_class(1) << p.get_ui()) - 1 ); for( mpz_class i(3); i <= p; ++i ) { s = (s * s - mpz_class(2)) % div ; }
return ( s == mpz_class(0) ); }
}
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; } }
}</lang>
- 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#
<lang csharp>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(); } }
}</lang>
- 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
<lang lisp>(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)))</lang>
- 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
<lang d>import std.stdio, std.math, std.bigint;
bool isPrime(in uint p) pure nothrow @safe @nogc {
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 uint p) pure nothrow /*@safe*/ {
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; }
}</lang>
- 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. <lang delphi>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();</lang>
- Output:
M2 M3 M5 M7 M13 M17 M19 M31
Erlang
<lang 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.</lang>
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: <lang fsharp>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 ]</lang>
Tail-recursive version: <lang fsharp>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
</lang>
Version using library folding function (way shorter and faster than the above): <lang fsharp>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
</lang>
Fortran
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. <lang fortran>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</lang>
FunL
<lang funl>def mersenne( p ) =
if p == 2 then return true var s = 4 var M = 2^p - 1 repeat p - 2 s = (s*s - 2) mod M
s == 0
import integers.primes
for p <- primes().filter( mersenne ).take( 20 )
println( 'M' + p )</lang>
- Output:
M2 M3 M5 M7 M13 M17 M19 M31 M61 M89 M107 M127 M521 M607 M1279 M2203 M2281 M3217 M4253 M4423
GAP
<lang 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]</lang>
Go
Processing the first list indicates that the test works. Processing the second shows it working on some larger numbers. <lang go>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) } }
}</lang>
- 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
<lang haskell>module Main
where
main = printMersennes $ take 45 $ filter lucasLehmer $ sieve [2..]
s mp 1 = 4 `mod` mp s mp n = ((s mp $ n-1)^2-2) `mod` mp
lucasLehmer 2 = True lucasLehmer p = s (2^p-1) (p-1) == 0
printMersennes = mapM_ (\x -> putStrLn $ "M" ++ show x)</lang> 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. <lang haskell>sieve (p:xs) = p : sieve [x | x <- xs, x `mod` p > 0]</lang> 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
<lang HicEst>s = 0 DO 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;', n
ENDDO
END</lang>
J
See Primality Tests essay on the J wiki.
Java
We use arbitrary-precision integers in order to be able to test any arbitrary prime.
<lang java>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(); }
}</lang>
- Output:
(after about eight hours)
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. <lang Mathematica>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}</lang>
This version is unbounded (and timed): <lang Mathematica>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}</lang>
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. <lang MATLAB>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</lang>
- Output:
<lang MATLAB>[mNumber,mersennesPrime] = mersennePrimes
mNumber =
2 3 5 7 13 17 19
mersennesPrime =
3 7 31 127 8191 131071 524287</lang>
Maxima
<lang maxima>lucas_lehmer(p) := block([s, n, i],
if not primep(p) then false elseif p = 2 then true else (s: 4, n: 2^p - 1, for i: 2 thru p - 1 do s: mod(s*s - 2, n), is(s = 0))
)$
sublist(makelist(i, i, 1, 200), lucas_lehmer); /* [2, 3, 5, 7, 13, 17, 19, 31, 61, 89, 107, 127] */</lang>
Modula-3
Modula-3 uses L as the literal for LONGINT. <lang modula3>MODULE LucasLehmer EXPORTS Main;
IMPORT IO, Fmt, Long;
PROCEDURE Mersenne(p: CARDINAL): BOOLEAN =
VAR s := 4L; m := Long.Shift(1L, p) - 1L; (* 2^p - 1 *) BEGIN IF p = 2 THEN RETURN TRUE; ELSE FOR i := 3 TO p DO s := (s * s - 2L) MOD m; END; RETURN s = 0L; END; END Mersenne;
BEGIN
FOR i := 2 TO 63 DO IF Mersenne(i) THEN IO.Put("M" & Fmt.Int(i) & " "); END; END; IO.Put("\n");
END LucasLehmer.</lang>
- Output:
M2 M3 M5 M7 M13 M17 M19 M31
Nimrod
<lang nimrod>import math
proc isPrime(a: int): bool =
if a == 2: return true if a < 2 or a mod 2 == 0: return false for i in countup(3, int sqrt(float a), 2): if a mod i == 0: return false return true
proc isMersennePrime(p: int): bool =
if p == 2: return true let mp = (1'i64 shl p) - 1 var s = 4'i64 for i in 3 .. p: s = (s * s - 2) mod mp result = s == 0
let upb = int((log2 float int64.high) / 2) echo " Mersenne primes:" for p in 2 .. upb:
if isPrime(p) and isMersennePrime(p): stdout.write " M",p
echo ""</lang>
- Output:
Mersenne primes: M2 M3 M5 M7 M13 M17 M19 M31
Oz
Oz's multiple precision number system use GMP core. <lang oz>%% compile : ozc -x <file.oz> functor import
Application System
define
fun {Arg Idx Default} Cmd = {Application.getArgs plain} Len = {Length Cmd} in if Len < Idx then Default else {StringToInt {Nth Cmd Idx}} end end fun {LLtest N} Mp = {Pow 2 N} - 1 fun {S K} X T in if K == 1 then 4 else T = {S K-1} X = T * T - 2 X mod Mp end end in if N == 2 then true else {S N-1} == 0 end end proc {FindLL X} fun {Sieve Ls} case Ls of nil then nil [] X|Xs then fun {DIV M} M mod X \= 0 end in X|{Sieve {Filter Xs DIV}} end end in if {IsList X} then case X of nil then skip [] M|Ms then {System.printInfo "M"#M#" "} {FindLL Ms} end else {FindLL {Filter {Sieve 2|{List.number 3 X 2}} LLtest}} end end Num = {Arg 1 607}
{FindLL Num} {Application.exit 0}
end</lang>
PARI/GP
<lang parigp>LL(p)={
my(m=Mod(4,1<<p-1)); for(i=3,p,m=m^2-2); m==0
};
search()={
my(t=1); print("2^2-1"); forprime(p=3,43112609, if(LL(p), print("2^"p"-1"); if(t++==47,return()) ) )
};</lang>
Pascal
int64 is good enough up to M31: <lang pascal>Program LucasLehmer(output); var
s, n: int64; i, exponent: integer;
begin
n := 1; for exponent := 2 to 31 do begin if exponent = 2 then s := 0 else s := 4; n := (n + 1)*2 - 1; // This saves from needing the math unit for exponentiation for i := 1 to exponent-2 do s := (s*s - 2) mod n; if s = 0 then writeln('M', exponent, ' is PRIME!'); end;
end.</lang>
- Output:
:> ./LucasLehmer M2 is PRIME! M3 is PRIME! M5 is PRIME! M7 is PRIME! M13 is PRIME! M17 is PRIME! M19 is PRIME! M31 is PRIME!
Perl
Using Math::GMP: <lang perl>use Math::GMP qw/:constant/;
sub is_prime { Math::GMP->new(shift)->probab_prime(12); }
sub is_mersenne_prime {
my $p = shift; return 1 if $p == 2; my $mp = 2 ** $p - 1; my $s = 4; $s = ($s * $s - 2) % $mp for 3..$p; $s == 0;
}
foreach my $p (2 .. 38_000_000) {
print "M$p\n" if is_prime($p) && is_mersenne_prime($p);
}</lang>
Another possibility is to use the modular lucas sequence from Math::Prime::Util, though it will be a bit slower as it calculates both and : <lang perl>use Math::Prime::Util qw/:all/; use bigint try=>"GMP,Pari"; forprimes {
my $p = $_; my $mp1 = 2**$p; print "M$p\n" if $p == 2 || 0 == (lucas_sequence($mp1-1, 4, 1, $mp1))[0];
} 38_000_000;</lang>
Perl's Math::BigInt
core module is a class that implements arbitrary-precision integers.
Here we use the bigint
pragma to transparently have number constants and operations be bigints:
<lang perl>sub is_prime {
my $p = shift; if ($p == 2) { return 1; } elsif ($p <= 1 || $p % 2 == 0) { return 0; } else { my $limit = sqrt($p); for (my $i = 3; $i <= $limit; $i += 2) { return 0 if $p % $i == 0; } return 1; }
}
sub is_mersenne_prime {
use bigint; my $p = shift; if ($p == 2) { return 1; } else { my $m_p = 2 ** $p - 1; my $s = 4; foreach my $i (3 .. $p) { $s = ($s ** 2 - 2) % $m_p; } return $s == 0; }
}
my $precision = 20000; # maximum requested number of decimal places of 2 ** MP-1 # my $long_bits_width = $precision / log(2) * log(10); my $upb_prime = int(($long_bits_width - 1)/2); # no unsigned # my $upb_count = 45; # find 45 mprimes if int was given enough bits #
print " Finding Mersenne primes in M[2..$upb_prime]:\n";
my $count = 0; foreach my $p (2 .. $upb_prime) {
if (is_prime($p) && is_mersenne_prime($p)) { print "M$p\n"; $count++; } last if $count >= $upb_count;
}</lang>
Perl 6
<lang perl6>multi is_mersenne_prime(2) { True } multi is_mersenne_prime(Int $p) {
my $m_p = 2 ** $p - 1; my $s = 4; # Alternate but slightly slower: $s = ($s * $s - 2) % $m_p for 3..$p; for (3 .. $p) { $s = $s.expmod(2, $m_p) - 2; $s += $m_p if $s < 0; } $s == 0;
}
say "M$_" if is-prime($_) and is_mersenne_prime($_) for 2..*;</lang>
- Output:
M2 M3 M5 M7 M13 M17 M19 M31 M61 M89 M107 M127 M521 M607 M1279 M2203 M2281 M3217 M4253 M4423 ^C
PicoLisp
<lang 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) ) ) )</lang>
- 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.
<lang pop11>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;</lang>
- Output:
(obtained in few seconds)
<lang pop11>M2 M3 M5 M7 M13 M17 M19 M31 M61 M89 M107 M127 M521 M607</lang>
PureBasic
PureBasic has no large integer support. Calculations are limited to the range of a signed quad integer type. <lang PureBasic>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 #False
EndProcedure
- upperBound = SizeOf(Quad) * 8 - 1 ;equivalent to significant bits in a signed quad integer
If 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</lang>
- Output:
M2 M3 M5 M7 M13 M17 M19 M31
Python
<lang python>from sys import stdout from 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=0 for 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: break
print</lang>
- 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
<lang python>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</lang>
Racket
<lang 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) </lang>
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 isn't conducive in searching for large primes.
<lang rexx>/*REXX program to use Lucas-Lehmer primality test for prime powers of 2.*/
parse arg limit . /*get the optional arg from C.L. */
if limit== then limit=1000 /*No argument? Then assume 1000.*/
list= /*placeholder for the results. */
/* [↓] only process up to LIMIT,*/ do j=1 by 2 to limit /*···only so many hours in a day.*/ power=j + (j==1) /*POWER ≡ J except for when J=1.*/ 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 all 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) /*right-justify 'em to look nice.*/ end /*k*/
exit /*stick a fork in it, we're done.*/ /*──────────────────────────────────ISPRIME subroutine──────────────────*/ isPrime: procedure; parse arg x /*get # to be tested*/ if 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 /*right-most dig ≡ 5? Not prime.*/ if x//7==0 then return 0 /*divisible by seven? Not prime.*/
do j=11 by 6 until j*j>x /*ensures J isn't divisible by 3.*/ if x// j ==0 then return 0 /*is divisible by J ? */ if x//(j+2)==0 then return 0 /* " " " J+2 ? ___*/ end /*j*/ /* [↑] perform loop through √ x */
return 1 /*indicate the number X is prime.*/ /*──────────────────────────────────LUCAS_LEHMER2 subroutine────────────*/ Lucas_Lehmer2: procedure; parse arg ? /*Lucas-Lehmer test on 2**? - 1 */ if form()\=='SCIENTIFIC' then numeric form /*ensure correct # form.*/ if ?==2 then s=0 /*handle the special even case.*/
else s=4
q=2**? /*╔═════════════════════════════════════════════════════════════╗
║Compute a power of 2, using only 9 decimal digits. DIGITs ║ ║of 1 million could be used, but that really gums up the whole║ ║works. So, we start with the default of 9 digits, find the ║ ║ten'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 squaring of S (below, for s*s).║ ╚═════════════════════════════════════════════════════════════╝*/
if pos('E',q)\==0 then do /*is # 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
do ?-2 /*apply, rinse, repeat ··· */ s=(s*s-2) // q /*remainder in REXX is: // */ end /* [↑] compute the real McCoy. */
if s\==0 then return /*return nuttin' if not prime. */
return 'M'? /*return modified (prime) number.*/</lang>
- Output:
when the following is used for input
══════════════════════════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 time) interpreted 100 seconds compiled 1.5 seconds (for limit=2000) 20 seconds --[[User:Walterpachl|Walterpachl]] 09:55, 12 July 2012 (UTC)
RPL
<lang 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;
\>> </lang>
- 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
<lang 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 true
end
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 == 0
end
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 = 0 for p in 2..upb_prime
if is_prime(p) && is_mersenne_prime(p) print "M%d " % p count += 1 end break if count >= upb_count
end puts</lang>
- 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
In accordance with definition of Mersenne primes it could only be a Mersenne number with prime exponent. <lang Scala>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 = 9941 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!")
}</lang>
- 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
<lang 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)))))</lang>
M2 M3 M5 M7 M13...
Scilab
<lang> iexpmax=30
n=1 for iexp=2:iexpmax if iexp==2 then s=0; else s=4; end n=(n+1)*2-1 for i=1:iexp-2 s=modulo((s*s-2),n) end if s==0 then printf("M%d ",iexp); end end</lang>
- Output:
M2 M3 M5 M7 M13 M17 M19
Seed7
To get maximum speed the program should be compiled with -O2.
<lang seed7>$ 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;</lang>
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
<lang Tcl>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</lang>
- 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
TI-83 BASIC
<lang ti83b>19→M 1→N For(E,2,M) If E=2 Then:0→S Else:4→S End (N+1)*2-1→N For(I,1,E-2) Reminder(S*S-2,N)→S End If S=0 Then:Disp E End End</lang>
- Output:
2 3 5 7 13 17 19
VBScript
<lang vb>iexpmax = 15 n=1 out="" For iexp = 2 To iexpmax If iexp = 2 Then s = 0 Else s = 4 End If n = (n + 1) * 2 - 1 For i = 1 To iexp - 2 s = (s * s - 2) Mod n Next If s = 0 Then out=out & "M" & iexp & " " End If Next Wscript.echo out</lang>
- Output:
M2 M3 M5 M7 M13
Visual Basic .NET
<lang vbnet>Public Class LucasLehmer
Private Sub btnGo_Click(sender As Object, e As EventArgs) Handles btnGo.Click Const iexpmax = 31 Dim s, n As Long Dim i, iexp As Integer n = 1 txtOut.Text = "" For iexp = 2 To iexpmax If iexp = 2 Then s = 0 Else s = 4 End If n = (n + 1) * 2 - 1 For i = 1 To iexp - 2 s = (s * s - 2) Mod n Next i If s = 0 Then txtOut.Text = txtOut.Text & "M" & iexp & " " End If Next iexp End Sub
End Class</lang>
- Output:
M2 M3 M5 M7 M13 M17 M19 M31
zkl
Using Extensible prime generator#zkl and the GMP library. <lang zkl>var BN=Import("zklBigNum"); // lib GMP primes:=Utils.Generator(Import("sieve").postponed_sieve); fcn isMersennePrime(p){
if(p==2) return(True); mp:=BN(1).shiftLeft(p) - 1; // 2^p - 1, a BIG number, like 1000s of digits s:=BN(4); do(p-2){ s.mul(s).sub(2).mod(mp) } // the % REALLY cuts down on mem usage return(s==0);
}</lang> Calculating S(n) is done in place (overwriting the value in the BigInt with the result); this really cuts down on memory usage. <lang>mersennePrimes:=primes.tweak(fcn(p){ isMersennePrime(p) and p or Void.Skip }); println("Mersenne primes:"); foreach mp in (mersennePrimes) { print(" M",mp); }</lang> This will "continuously" spew out Mersenne Primes.
Tweaking a Walker (aka iterator, Generators are a class of Walker) basically puts a filter on the underlying iterator, in this case, ignoring prime numbers that are not Mersenne primes and passing those that are.
- Output:
Mersenne primes: 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
Additionally, this problem is readily threaded and has a linear speedup. Since there are lots of calculations between results, the [bigger] results are basically time sorted. However, N times faster doesn't mean much given the huge calculations used by the LL test (math with thousands of digits ain't quick).
Using five threads: <lang zkl>ps:=Thread.Pipe(); mpOut:=Thread.Pipe(); // how the threads will communicate fcn(ps){ // a thread to generate primes, sleeps most of the time
Utils.Generator(Import("sieve").postponed_sieve).pump(ps)
}.launch(ps);
do(4){ // four threads to perform the Lucas-Lehmer test
fcn(ps,out){ ps.filter2s(out,isMersennePrime) }.launch(ps,mpOut)
} println("Mersenne primes:"); foreach mp in (mpOut) { print(" M",mp); }</lang>
- Programming Tasks
- Prime Numbers
- Arbitrary precision
- Arithmetic operations
- 360 Assembly
- Ada
- Agena
- ALGOL 68
- AWK
- BBC BASIC
- Bracmat
- C
- GMP
- C++
- C sharp
- Clojure
- D
- DWScript
- Erlang
- F Sharp
- Fortran
- FunL
- GAP
- Go
- Haskell
- HicEst
- J
- Java
- Mathematica
- MATLAB
- Maxima
- Modula-3
- Nimrod
- Oz
- PARI/GP
- Pascal
- Perl
- Perl 6
- PicoLisp
- Pop11
- PureBasic
- Python
- Racket
- REXX
- RPL
- Ruby
- Scala
- Scheme
- Scilab
- Seed7
- Tcl
- TI-83 BASIC
- VBScript
- Visual Basic .NET
- Zkl
- GUISS/Omit