Lucas-Lehmer test: Difference between revisions
Content added Content deleted
(→{{header|ALGOL 68}}: calculate all Mersenne primes up to the implementation's max precision, or the 45th Mersenne prime. (Which ever comes first).) |
m (→{{header|Python}}: straight from A68 ☺ ... python does rather well with arbitary precision, fast too.) |
||
Line 188: | Line 188: | ||
} |
} |
||
} |
} |
||
Output: |
Output: |
||
Finding Mersenne primes in M[2..500]: |
Finding Mersenne primes in M[2..500]: |
||
M2 M3 M5 M7 M13 M17 M19 M31 M61 M89 M107 M127 |
M2 M3 M5 M7 M13 M17 M19 M31 M61 M89 M107 M127 |
||
=={{header|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 ): |
|||
# print "p:",p |
|||
if p % i == 0 : return False |
|||
else: |
|||
return True |
|||
fi |
|||
fi; |
|||
def is_mersenne_prime ( p ): |
|||
if p == 2 : return True |
|||
else: |
|||
m_p = 2 ** p - 1; s = 4; |
|||
for i in range(3, p+1): |
|||
s = (s * s - 2) % m_p |
|||
od; |
|||
return s == 0 |
|||
fi; |
|||
precision = 20000; |
|||
long_long_bits_width = precision / log(2) * log(10); |
|||
upb_prime = int( long_long_bits_width - 1 ) / 2; # no unsigned # |
|||
upb_count = 45; # find 45 mprimes if int has 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) : |
|||
if is_mersenne_prime(p) : |
|||
stdout.write(" M%d"%p); |
|||
stdout.flush(); |
|||
count += 1 |
|||
fi |
|||
fi; |
|||
if count >= upb_count : break; fi |
|||
od |
|||
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 |
Revision as of 11:19, 16 February 2008
![Task](http://static.miraheze.org/rosettacodewiki/thumb/b/ba/Rcode-button-task-crushed.png/64px-Rcode-button-task-crushed.png)
You are encouraged to solve this task according to the task description, using any language you may know.
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; 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=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 * s - 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
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
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 ): # print "p:",p if p % i == 0 : return False else: return True fi fi; def is_mersenne_prime ( p ): if p == 2 : return True else: m_p = 2 ** p - 1; s = 4; for i in range(3, p+1): s = (s * s - 2) % m_p od; return s == 0 fi; precision = 20000; long_long_bits_width = precision / log(2) * log(10); upb_prime = int( long_long_bits_width - 1 ) / 2; # no unsigned # upb_count = 45; # find 45 mprimes if int has 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) : if is_mersenne_prime(p) : stdout.write(" M%d"%p); stdout.flush(); count += 1 fi fi; if count >= upb_count : break; fi od
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