Lucas-Lehmer test

From Rosetta Code

Jump to: navigation, search
Lucas-Lehmer test is a programming task. Visitors like you are encouraged to solve it according to the task description, using any language they may happen to know.
Add to BlogMarksAdd to del.icio.usAdd to diggAdd to NewsvineAdd to redditAdd to Slashdot

Lucas-Lehmer Test: for p a prime, the Mersenne number 2p − 1 is prime if and only if 2p − 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 maximum precision, or the 47th 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] Agena

Because of the very large numbers computed, the mapm binding is used to calculate with arbitrary precision. 

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;

produces: 

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] Clojure

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

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


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

[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 = mapM_ (\x -> putStrLn $ "M" ++ show x)

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] J

See Primality Tests essay on the J wiki.

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

This version is very speedy and is bounded. 

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}

This version is unbounded (and timed): 

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}

I'll see what this gets.

[edit] Modula-3

Modula-3 uses L as the literal for LONGINT. 

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.

Output: 

M2 M3 M5 M7 M13 M17 M19 M31

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

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

[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 # 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

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

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

[edit] Tcl

Modeled after the Pop11 solution. 

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

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