Lucas-Lehmer test: Difference between revisions

Content added Content deleted

Inline

Revision as of 14:35, 21 February 2008

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

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

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

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

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);
		
		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);
	
	for(mpz_class i(1); i < exp; ++i)
		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

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 = 2; p <= upb; p++)
            if (isPrime(p) && isMersennePrime(p))
                System.out.print(" M" + p);
        System.out.println();
    }
}

Output:

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

Python

from sys import stdout
from math import sqrt, log
fi = od = None

def is_prime ( p ):
  if p == 2: return True
  elif p <= 1 or p % 2 == 0: return False
  else:
    prime = True;
    for i in range(3, int(sqrt(p))+1, 2 ): # do
      if p % i == 0: return False
    else:
      return True
    fi
  fi;

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): # do
      s = (s ** 2 - 2) % m_p
    od;
    return s == 0
  fi;

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): # do
  if is_prime(p) and is_mersenne_prime(p):
    print("M%d"%p),
    stdout.flush();
    count += 1
  fi;
  if count >= upb_count: break
od
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

@@ Line 40: / Line 40: @@
 =={{header|ALGOL 68}}==
-Interpretor: algol68g-mk11
+{{works with|algol68g-mk11}}
  main:(
    PRAGMAT stack=1M precision=20000 PRAGMAT
@@ Line 84: / Line 84: @@
  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
 =={{header|C}}==
 {{works with|gcc|4.1.2 20070925 (Red Hat 4.1.2-27)}}