Lucas-Lehmer test

From Rosetta Code

Jump to: navigation, search

Programming Task
This is a programming task. It lays out a problem which Rosetta Code users are encouraged to solve, using languages they know.

Code examples should be formatted along the lines of one of the existing prototypes.

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 45th Mersenne prime. (Which ever comes first).

Contents

[edit] 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

[edit] 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

[edit] 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

#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

[edit] C++

Library: GMP 


#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

[edit] 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

[edit] 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

[edit] Oz

Oz's multiple precision number system use GMP core. 

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

[edit] 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

[edit] 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

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

[edit] 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...
Retrieved from "http://rosettacode.org/wiki/Lucas-Lehmer_test"
Personal tools