Lucas-Lehmer test
From Rosetta Code
You are encouraged to solve this task according to the task description, using any language you may know.
Lucas-Lehmer Test: for p an odd 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).
[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
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;
test:(
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
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++
#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 M44497
[edit] C#
using System;
using System.Collections.Generic;
using System.Numerics;
using System.Threading.Tasks;
namespace LucasLehmerTestForRosettaCode
{
public class LucasLehmerTest
{
static BigInteger ZERO = new BigInteger(0);
static BigInteger ONE = new BigInteger(1);
static BigInteger TWO = new BigInteger(2);
static BigInteger FOUR = new BigInteger(4);
private static bool isMersennePrime(int p)
{
if (p % 2 == 0) return (p == 2);
else {
for (int i = 3; i <= (int)Math.Sqrt(p); i += 2)
if (p % i == 0) return false; //not prime
BigInteger m_p = BigInteger.Pow(TWO, p) - ONE;
BigInteger s = FOUR;
for (int i = 3; i <= p; i++)
s = (s * s - TWO) % m_p;
return s == ZERO;
}
}
public static int[] GetMersennePrimeNumbers(int upTo)
{
List<int> response = new List<int>();
Parallel.For(2, upTo, i => {
if (isMersennePrime(i)) response.Add(i);
});
response.Sort();
return response.ToArray();
}
static void Main(string[] args)
{
int[] mersennePrimes = LucasLehmerTest.GetMersennePrimeNumbers(11213);
foreach (int mp in mersennePrimes)
Console.Write("M" + mp+" ");
Console.ReadLine();
}
}
}
Output: (Run only to 11213)
M2 M3 M5 M7 M13 M17 M19 M31 M61 M89 M107 M127 M521 M607 M1279 M2203 M2281 M3217 M4253 M4423 M9689 M9941 M11213
[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] D
import std.stdio, std.math, std.bigint;
bool isPrime(int p) pure nothrow {
if (p < 2 || p % 2 == 0)
return p == 2;
foreach (i; 3 .. cast(uint)(sqrt(cast(real)p)) + 1)
if (p % i == 0)
return false;
return true;
}
bool isMersennePrime(int p) /*pure nothrow*/ {
if (!isPrime(p))
return false;
if (p == 2)
return true;
auto mp = (BigInt(1) << p) - 1;
auto s = BigInt(4);
foreach (_; 3 .. p + 1)
s = (s ^^ 2 - 2) % mp;
return s == 0;
}
void main() {
foreach (p; 2 .. 2_300)
if (isMersennePrime(p)) {
write(" M", p);
stdout.flush();
}
}
- Output:
M2 M3 M5 M7 M13 M17 M19 M31 M61 M89 M107 M127 M521 M607 M1279 M2203 M2281
With p up to 10_000 it prints:
M2 M3 M5 M7 M13 M17 M19 M31 M61 M89 M107 M127 M521 M607 M1279 M2203 M2281 M3217 M4253 M4423 M9941
[edit] DWScript
Using Integer type, which is 64bit, limits the search to M31.
function IsMersennePrime(p : Integer) : Boolean;
var
i, s, m_p : Integer;
begin
if p=2 then
Result:=True
else begin
m_p := (1 shl p)-1;
s := 4;
for i:=3 to p do
s:=(s*s-2) mod m_p;
Result:=(s=0);
end;
end;
const upperBound = Round(Log2(High(Integer))/2);
PrintLn('Finding Mersenne primes in M[2..' + IntToStr(upperBound) + ']: ');
Print('M2');
var p : Integer;
for p:=3 to upperBound step 2 do begin
if IsMersennePrime(p) then
Print(' M'+IntToStr(p));
end;
PrintLn('');
Output:
M2 M3 M5 M7 M13 M17 M19 M31
[edit] Erlang
-module(mp).
-export([main/0]).
main() -> [ io:format("M~p ", [P]) || P <- lists:seq(2,700), (P == 2) orelse (s((1 bsl P) - 1, P-1) == 0) ].
s(MP,1) -> 4 rem MP;
s(MP,N) -> X=s(MP,N-1), (X*X - 2) rem MP.
In 3 seconds will print
M2 M3 M5 M7 M13 M17 M19 M31 M61 M89 M107 M127 M521 M607
Testing larger numbers (i.e. 5000) is possible but will take few minutes.
[edit] F#
Simple arbitrary-precision version:
let rec s mp n =
if n = 1 then 4I % mp else ((s mp (n - 1)) ** 2 - 2I) % mp
[ for p in 2..47 do
if p = 2 || s ((1I <<< p) - 1I) (p - 1) = 0I then
yield p ]
[edit] Fortran
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] Go
Processing the first list indicates that the test works. Processing the second shows it working on some larger numbers.
package main
import (
"fmt"
"math/big"
)
var primes = []uint{3, 5, 7, 11, 13, 17, 19, 23, 29, 31, 37, 41, 43, 47,
53, 59, 61, 67, 71, 73, 79, 83, 89, 97, 101, 103, 107, 109, 113, 127}
var mersennes = []uint{521, 607, 1279, 2203, 2281, 3217, 4253, 4423, 9689,
9941, 11213, 19937, 21701, 23209, 44497, 86243, 110503, 132049, 216091,
756839, 859433, 1257787, 1398269, 2976221, 3021377, 6972593, 13466917,
20996011, 24036583}
func main() {
llTest(primes)
fmt.Println()
llTest(mersennes)
}
func llTest(ps []uint) {
var s, m big.Int
one := big.NewInt(1)
two := big.NewInt(2)
for _, p := range ps {
m.Sub(m.Lsh(one, p), one)
s.SetInt64(4)
for i := uint(2); i < p; i++ {
s.Mod(s.Sub(s.Mul(&s, &s), two), &m)
}
if s.BitLen() == 0 {
fmt.Printf("M%d ", p)
}
}
}
- Output:
M3 M5 M7 M13 M17 M19 M31 M61 M89 M107 M127 M521 M607 M1279 M2203 M2281 M3217 M4253 M4423 M9689 M9941 M11213 M19937...
[edit] Haskell
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] HicEst
s = 0
DO exponent = 2, 31
IF(exponent > 2) s = 4
n = 2^exponent - 1
DO i = 1, exponent-2
s = MOD(s*s - 2, n)
ENDDO
IF(s == 0) WRITE(Messagebox) 'M', exponent, ' is prime;', n
ENDDO
END
[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] MATLAB
MATLAB suffers from a lack of an arbitrary precision math (bignums) library. It also doesn't have great support for 64-bit integer arithmetic...or at least MATLAB 2007 doesn't. So, the best precision we have is doubles; therefore, this script can only find up to M19 and no greater.
function [mNumber,mersennesPrime] = mersennePrimes()
function isPrime = lucasLehmerTest(thePrime)
llResidue = 4;
mersennesPrime = (2^thePrime)-1;
for i = ( 1:thePrime-2 )
llResidue = mod( ((llResidue^2) - 2),mersennesPrime );
end
isPrime = (llResidue == 0);
end
%Because IEEE764 Double is the highest precision number we can
%represent in MATLAB, the highest Mersenne Number we can test is 2^52.
%In addition, because we have this cap, we can only test up to the
%number 30 for Mersenne Primeness. When we input 31 into the
%Lucas-Lehmer test, during the computation of the residue, the
%algorithm multiplies two numbers together the result of which is
%greater than 2^53. Because we require every digit to be significant,
%this leads to an error. The Lucas-Lehmer test should say that M31 is a
%Mersenne Prime, but because of the rounding error in calculating the
%residues caused by floating-point arithmetic, it does not. So M30 is
%the largest number we test.
mNumber = (3:30);
[isPrime] = arrayfun(@lucasLehmerTest,mNumber);
mNumber = [2 mNumber(isPrime)];
mersennesPrime = (2.^mNumber) - 1;
end
Ouput:
[mNumber,mersennesPrime] = mersennePrimes
mNumber =
2 3 5 7 13 17 19
mersennesPrime =
3 7 31 127 8191 131071 524287
[edit] Maxima
lucas_lehmer(p) := block([s, n, i],
if not primep(p) then false elseif p = 2 then true else
(s: 4,
n: 2^p - 1,
for i: 2 thru p - 1 do s: mod(s*s - 2, n),
is(s = 0))
)$
sublist(makelist(i, i, 1, 200), lucas_lehmer);
/* [2, 3, 5, 7, 13, 17, 19, 31, 61, 89, 107, 127] */
[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] PARI/GP
LL(p)={
my(m=Mod(4,1<<p-1));
for(i=3,p,m=m^2-2);
m==0
};
search()={
my(t=1);
print("2^2-1");
forprime(p=3,43112609,
if(LL(p),
print("2^"p"-1");
if(t++==47,return())
)
)
};
[edit] Pascal
int64 is good enough up to M31:
Program LucasLehmer(output);
var
s, n: int64;
i, exponent: integer;
begin
n := 1;
for exponent := 2 to 31 do
begin
if exponent = 2 then
s := 0
else
s := 4;
n := (n + 1)*2 - 1; // This saves from needing the math unit for exponentiation
for i := 1 to exponent-2 do
s := (s*s - 2) mod n;
if s = 0 then
writeln('M', exponent, ' is PRIME!');
end;
end.
Output:
:> ./LucasLehmer M2 is PRIME! M3 is PRIME! M5 is PRIME! M7 is PRIME! M13 is PRIME! M17 is PRIME! M19 is PRIME! M31 is PRIME!
[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] Perl 6
Precision limited to 18 because rakudo does not yet implement arbitrary precision Int as specced.
multi is_prime(Int $p where ($p <= 1)) { False }
multi is_prime(Int $p) { $p %% none(2,3,5,7...^ * > sqrt($p)) }
multi is_mersenne_prime(2) { True }
multi is_mersenne_prime(Int $p) {
my $m_p = 2 ** $p - 1;
my $s = 4;
$s = ($s ** 2 - 2) % $m_p for 3 .. $p;
$s == 0;
}
my $precision = 18; # maximum requested number of decimal places of 2 ** MP-1
my $long_bits_width = $precision / log(2) * log(10);
my $max_prime = floor(($long_bits_width - 1)/2);
my $max_count = 45;
say " Finding Mersenne primes in M[2..$max_prime]:";
my $count = 0;
for 2 .. $max_prime -> $p {
if is_prime($p) && is_mersenne_prime($p) {
say "M$p";
last if ++$count >= $max_count;
}
}
[edit] PicoLisp
(de prime? (N)
(or
(= N 2)
(and
(> N 1)
(bit? 1 N)
(for (D 3 T (+ D 2))
(T (> D (sqrt N)) T)
(T (=0 (% N D)) NIL) ) ) ) )
(de mersenne? (P)
(or
(= P 2)
(let (MP (dec (>> (- P) 1)) S 4)
(do (- P 2)
(setq S (% (- (* S S) 2) MP)) )
(=0 S) ) ) )
Output:
: (for N 10000 (and (prime? N) (mersenne? N) (println N)) ) 2 3 5 7 13 17 19 31 61 89 107 127 521 607 1279 2203 2281 3217 4253 4423 9689 9941
[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] PureBasic
PureBasic has no large integer support. Calculations are limited to the range of a signed quad integer type.
Procedure Lucas_Lehmer_Test(p)
Protected mp.q = (1 << p) - 1, sn.q = 4, i
For i = 3 To p
sn = (sn * sn - 2) % mp
Next
If sn = 0
ProcedureReturn #True
EndIf
ProcedureReturn #False
EndProcedure
#upperBound = SizeOf(Quad) * 8 - 1 ;equivalent to significant bits in a signed quad integer
If OpenConsole()
Define p = 3
PrintN("M2")
While p <= #upperBound
If Lucas_Lehmer_Test(p)
PrintN("M" + Str(p))
EndIf
p + 2
Wend
Print(#CRLF$ + #CRLF$ + "Press ENTER to exit"): Input()
CloseConsole()
EndIf
Sample output:
M2 M3 M5 M7 M13 M17 M19 M31
[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(10, 2)
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
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] REXX
REXX won't have a problem with the large number of digits involved, but since it's an interpreted language, such massive number crunching is not conducive to find large primes.
/*REXX program to use Lucas-Lehmer primality test for prime powers of 2.*/
parse arg limit . /*get the argument (if any). */
if limit=='' then limit=1000 /*if no argument, assume 1000. */
list='' /*placeholder for the results. */
do j=1 by 2 to limit /*only process so much... */
/*there's only so many hours in a day...*/
power=j
if j==1 then power=2 /*special case for the even prime*/
if \isprime(power) then iterate /*if not prime, then ignore it. */
list=list Lucas_Lehmer2(power) /*add to list (or maybe not). */
end /*j*/
list=space(list) /*remove extraneous blanks. */
say
say center('list',60-3,"-") /*show a fancy-dancy header/title*/
say
do j=1 for words(list) /*show entries in list, 1/line. */
say right(word(list,j),30) /*and right-justify 'em. */
end
say
exit
/*-------------------------------------LUCAS_LEHMER2 subroutine---------*/
Lucas_Lehmer2: procedure; parse arg ? /*Lucas-Lehmer test on 2**? - 1 */
if ?==2 then s=0
else s=4
/* DIGITs of 1 million could be used, but that really */
/* slows up the whole works. So, we start with the */
/* default of 9 digits, find the 10's exponent in */
/* the product (2**?), double it, and then add 6. */
/* 2 is all that's needed, but 6 is a lot safer. */
/* The doubling is for the square of S (s*s, below).*/
q=2**? /*compute a power of 2, using only 9 decimal digits. */
if pos('E',q)\==0 then do /*is it in exponential notation? */
parse var q 'E' tenpow
numeric digits tenpow*2+6
end
else numeric digits digits()*2+6 /* 9*2 + 6 */
q=2**?-1 /*now, compute the real McCoy. */
do ?-2 /*apply, rinse, repeat ... */
s=(s*s-2)//q /*modulus in REXX is: // */
end
if s\==0 then return '' /*return nuttin' if not prime. */
return 'M'? /*return modified (prime) number.*/
/*-------------------------------------ISPRIME subroutine---------------*/
isprime: procedure; parse arg x; if x<2 then return 0 /*idiot test.*/
if wordpos(x,'2 3 5 7')\==0 then return 1 /*test for special cases.*/
if x//3==0 then return 0 /*divisible by three? Not prime.*/
if x//2==0 then return 0 /*is it even? Not prime.*/
do j=5 by 6 /*ensures J isn't divisible by 3.*/
if x//j==0 | x//(j+2)==0 then return 0
if j*j>x then return 1 /*past the sqrt of X? It's prime*/
end
Output from the program when the following is specified:
10000
--------------------------list---------------------------
M2
M3
M5
M7
M13
M17
M19
M31
M61
M89
M107
M127
M521
M607
M1279
M2203
M2281
M3217
M4253
M4423
M9689
M9941
[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] Seed7
To get maximum speed the program should be compiled with -O2.
$ include "seed7_05.s7i";
include "bigint.s7i";
const func boolean: is_prime (in integer: number) is func
result
var boolean: prime is FALSE;
local
var integer: upTo is 0;
var integer: testNum is 3;
begin
if number = 2 then
prime := TRUE;
elsif number rem 2 = 0 or number <= 1 then
prime := FALSE;
else
upTo := sqrt(number);
while number rem testNum <> 0 and testNum <= upTo do
testNum +:= 2;
end while;
prime := testNum > upTo;
end if;
end func;
const func boolean: lucas_lehmer_test (in integer: p) is func
result
var boolean: prime is TRUE;
local
var bigInteger: m_p is 0_;
var bigInteger: s is 4_;
var integer: i is 0;
begin
if p <> 2 then
m_p := 2_ ** p - 1_;
for i range 2 to pred(p) do
s := (s ** 2 - 2_) rem m_p;
end for;
prime := s = 0_;
end if;
end func;
const proc: main is func
local
var integer: p is 2;
begin
writeln(" Mersenne primes:");
while p <= 3217 do
if is_prime(p) and lucas_lehmer_test(p) then
write(" M" <& p);
flush(OUT);
end if;
incr(p);
end while;
writeln;
end func;
Original source: [1]
Output:
Mersenne primes: M2 M5 M7 M13 M17 M19 M31 M61 M89 M107 M127 M521 M607 M1279 M2203 M2281 M3217
[edit] Tcl
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
