Lucas-Lehmer test: Difference between revisions

From Rosetta Code
Content added Content deleted
(added Fortran)
(added ruby)
Line 7:
 
=={{header|Ada}}==
<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;</ada>
Output:
Mersenne primes:
Line 89:
{{works with|C99}}
Compiler options: gcc -std=c99 -lm Lucas-Lehmer_test.c -o Lucas-Lehmer_test<br>
<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");
}</java>
}
 
Output:
Line 142:
{{libheader|GMP}}
 
<precpp>#include <iostream>
#include <gmpxx.h>
 
Line 198:
//std::cout << "pow2( " << exp << " ) = " << ret << std::endl;
return ret;
}</precpp>▼
}
▲</pre>
 
Output: (Incomplete; It takes a long time.)
Line 255 ⟶ 254:
We use arbitrary-precision integers in order to be able to test any arbitrary prime.
 
<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();
}
}</java>
}
Output (after about eight hours):
Finding Mersenne primes in M[2..2147483647]:
Line 416 ⟶ 415:
 
=={{header|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:
<pre> 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</pre>
 
=={{header|Ruby}}==
<python>def is_prime ( p )
if p == 2
return true
elsif p <= 1 or 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
for i in 3..p
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</python>
 
Output:
Line 458 ⟶ 509:
 
=={{header|Scheme}}==
<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)))))</scheme>
 
M2 M3 M5 M7 M13...

Revision as of 23:20, 6 December 2008

Task
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 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

<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;</ada> Output: 

Mersenne primes:
M2 M3 M5 M7 M13 M17 M19 M31

ALGOL 68

Works with: algol68g-mk11
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

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
<c>#include <math.h>

  1. include <stdio.h>
  2. include <limits.h>
  3. 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");

}</java>

Output: 

Mersenne primes:
M2 M3 M5 M7 M13 M17 M19 M31

C++

Library: GMP

<cpp>#include <iostream>

  1. 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;

}</cpp>

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

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

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` 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.

<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();
   }

}</java> 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

Ruby

<python>def is_prime ( p ) 

 if p == 2
   return true
 elsif p <= 1 or 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
   for i in 3..p
     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</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

<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)))))</scheme>

M2 M3 M5 M7 M13...
Retrieved from "https://rosettacode.org/wiki/Lucas-Lehmer_test?oldid=57858"