Lucas-Lehmer test: Difference between revisions
m (→[[Lucas-Lehmer_test#ALGOL 68]]: fix pre tag) |
m (→[[Lucas-Lehmer_test#ALGOL 68]]: rm new line) |
||
Line 71:
{{works with|ALGOL 68G|Any - tested with release [http://sourceforge.net/projects/algol68/files/algol68g/algol68g-1.18.0/algol68g-1.18.0-9h.tiny.el5.centos.fc11.i386.rpm/download 1.18.0-9h.tiny]}}
{{wont work with|ELLA ALGOL 68|Any (with appropriate job cards) - tested with release [http://sourceforge.net/projects/algol68/files/algol68toc/algol68toc-1.8.8d/algol68toc-1.8-8d.fc9.i386.rpm/download 1.8-8d] - due to extensive use of FORMATted transput}}
<lang algol68>PRAGMAT stack=1M precision=20000 PRAGMAT
PROC is prime = ( INT p )BOOL:
|
Revision as of 07:23, 22 July 2010
You are encouraged to solve this task according to the task description, using any language you may know.
Lucas-Lehmer Test: for a 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).
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
C
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++
<lang cpp>#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;
}</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
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...
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>
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>
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.println(" Finding Mersenne primes in M[2.." + upb + "]: "); 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>
Ouput: <lang MATLAB>[mNumber,mersennesPrime] = mersennePrimes
mNumber =
2 3 5 7 13 17 19
mersennesPrime =
3 7 31 127 8191 131071 524287</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
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>
Perl
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>
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>
The output (obtained in few seconds) is: <lang pop11>M2 M3 M5 M7 M13 M17 M19 M31 M61 M89 M107 M127 M521 M607</lang>
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(2) * log(10) 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
Ruby
<lang ruby>def is_prime ( p )
if p == 2 return true elsif p <= 1 || p % 2 == 0 return false else (3 .. Math.sqrt(p)).step(2) do |i| if p % i == 0 return false end end return true end
end
def is_mersenne_prime ( p )
if p == 2 return true else m_p = ( 1 << p ) - 1 s = 4 (p-2).times do s = (s ** 2 - 2) % m_p end return s == 0 end
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 if count >= upb_count break end
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
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...
Tcl
Modeled after the Pop11 solution.
<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> 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