Lucas-Lehmer test

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

Ada

with Ada.Text_Io; use Ada.Text_Io;
with Ada.Integer_Text_Io; use Ada.Integer_Text_Io;
with Ada.Numerics.Elementary_Functions; use Ada.Numerics.Elementary_Functions;

procedure Lucas_Lehmer_Test is
   function Mersenne(Item : Integer) return Boolean is
      S : Long_Long_Integer := 4;
      MP : Long_Long_Integer := 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 := Integer(Float'Rounding(Log(10.0) / Log(2.0) * 10_000_000.0));
   M_Count : Natural := 0;
begin
   Put_Line(" Mersenne primes:");
   for P in 2..Upper_Bound loop
      if Mersenne(P) then
         Put(" M");
         Put(Item => P, Width => 1);
         M_Count := M_Count + 1;
         exit when M_Count = 45;
      end if;
   end loop;
end Lucas_Lehmer_Test;

Output: (Output was truncated by an arithmetic overflow exception.)

Mersenne primes:
M2 M3 M5 M7 M13 M17 M19 M31

ALGOL 68

Interpretor: algol68g-mk11

main:(
  PRAGMAT stack=1M precision=20001 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 * s - 2) MOD m p
      OD;
      s = 0
    FI;

  INT upb = ENTIER(SHORTEN long log(SHORTEN LONG LONG REAL(long long max int))/log(2)/2);

  printf(($" Finding Mersenne primes in M[2.."g(0)"]: "l$,upb));

  FOR p FROM 2 TO upb DO
    IF is prime(p) THEN
      IF is mersenne prime(p) THEN
        printf (($" M"g(0)$,p))
      FI
    FI
  OD
)

Output:

Finding Mersenne primes in M[2..33253]: 
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

#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);
      s = s % 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");
}

Compiler version: gcc version 4.1.2 20070925 (Red Hat 4.1.2-27)
Compiler options: gcc -std=c99 -lm Lucas-Lehmer_test.c -o Lucas-Lehmer_test
Output:

Mersenne primes:
M2 M3 M5 M7 M13 M17 M19 M31

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..500]:
M2 M3 M5 M7 M13 M17 M19 M31 M61 M89 M107 M127