Lucas-Lehmer test

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

Lucas-Lehmer Test: for p 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...