Lucas-Lehmer test: Difference between revisions
(→{{header|Python}}: python colorization, remove nonpythonic syntax) |
(added Fortran) |
||
Line 4:
The following programs calculate all Mersenne primes up to the implementation's
maximium precision, or the
=={{header|Ada}}==
Line 203:
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
=={{header|Fortran}}==
{{works with|Fortran|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
=={{header|Haskell}}==
|
Revision as of 17:18, 2 December 2008
You are encouraged to solve this task according to the task description, using any language you may know.
Lucas-Lehmer Test: for p a prime, the Mersenne number 2**p-1 is prime if and only if 2**p-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 maximium precision, or the 47th Mersenne prime. (Which ever comes first).
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;
Output:
Mersenne primes: M2 M3 M5 M7 M13 M17 M19 M31
ALGOL 68
main:( 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; 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
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
#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++
#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
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.
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
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 [] = return () printMersennes (x:xs) = do putStrLn $ "M"++(show x) printMersennes xs
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
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.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(); } }
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
Oz
Oz's multiple precision number system use GMP core. <ocaml>%% 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</ocaml>
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:
M2 M3 M5 M7 M13 M17 M19 M31 M61 M89 M107 M127 M521 M607
Python
<python> from sys import stdout
from math import sqrt, log def is_prime ( p ): if p == 2: return True 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</python>
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
;;;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...