Lucas-Lehmer test: Difference between revisions
Content added Content deleted
(added Fortran) |
(added ruby) |
||
Line 7: | Line 7: | ||
=={{header|Ada}}== |
=={{header|Ada}}== |
||
<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;</ada> |
|||
Output: |
Output: |
||
Mersenne primes: |
Mersenne primes: |
||
Line 89: | Line 89: | ||
{{works with|C99}} |
{{works with|C99}} |
||
Compiler options: gcc -std=c99 -lm Lucas-Lehmer_test.c -o Lucas-Lehmer_test<br> |
Compiler options: gcc -std=c99 -lm Lucas-Lehmer_test.c -o Lucas-Lehmer_test<br> |
||
<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) % 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"); |
|||
}</java> |
|||
} |
|||
Output: |
Output: |
||
Line 142: | Line 142: | ||
{{libheader|GMP}} |
{{libheader|GMP}} |
||
< |
<cpp>#include <iostream> |
||
#include <gmpxx.h> |
#include <gmpxx.h> |
||
Line 198: | Line 198: | ||
//std::cout << "pow2( " << exp << " ) = " << ret << std::endl; |
//std::cout << "pow2( " << exp << " ) = " << ret << std::endl; |
||
return ret; |
return ret; |
||
⚫ | |||
} |
|||
⚫ | |||
Output: (Incomplete; It takes a long time.) |
Output: (Incomplete; It takes a long time.) |
||
Line 255: | Line 254: | ||
We use arbitrary-precision integers in order to be able to test any arbitrary prime. |
We use arbitrary-precision integers in order to be able to test any arbitrary prime. |
||
<java>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(); |
|||
} |
|||
}</java> |
|||
} |
|||
Output (after about eight hours): |
Output (after about eight hours): |
||
Finding Mersenne primes in M[2..2147483647]: |
Finding Mersenne primes in M[2..2147483647]: |
||
Line 416: | Line 415: | ||
=={{header|Python}}== |
=={{header|Python}}== |
||
<python> |
<python>from sys import stdout |
||
from math import sqrt, log |
|||
def is_prime ( p ): |
|||
if p == 2: return True |
|||
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</python> |
|||
Output: |
|||
<pre> 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</pre> |
|||
=={{header|Ruby}}== |
|||
<python>def is_prime ( p ) |
|||
if p == 2 |
|||
return true |
|||
elsif p <= 1 or 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 |
|||
for i in 3..p |
|||
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</python> |
|||
Output: |
Output: |
||
Line 458: | Line 509: | ||
=={{header|Scheme}}== |
=={{header|Scheme}}== |
||
<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)))))</scheme> |
|||
M2 M3 M5 M7 M13... |
M2 M3 M5 M7 M13... |