Mersenne primes: Difference between revisions
m (Fixed minor typo) |
|||
(46 intermediate revisions by 24 users not shown) | |||
Line 1: | Line 1: | ||
{{draft task}} |
{{draft task|Prime Numbers}} |
||
Mersenne primes: |
|||
;Task: |
|||
Challenge: |
|||
Create code that will list (preferably calculate) all of the ''Mersenne primes'' until some limitation is reached. |
|||
Create code that will list (preferably calculate) all of the Mersenne primes until some limitation is reached. For information on what a Mersenne prime is, go to this link: [[https://en.wikipedia.org/wiki/Mersenne_prime]] |
|||
The number of ''known'' Mersenne primes is '''51''' (as of June, 2020), and the largest known Mersenne prime contains contains '''24,862,048''' decimal digits. |
|||
;Also see: |
|||
* the Wikipedia entry: [https://en.wikipedia.org/wiki/Mersenne_prime Mersenne prime]. |
|||
* the MathWorld entry; [https://mathworld.wolfram.com/MersennePrime.html Mersenne prime]. |
|||
* For a list of all the know Mersenne primes: [https://primes.utm.edu/mersenne/index.html#known list of Mersenne |
|||
* For a general website about primes: [https://primes.utm.edu/ prime pages]. |
|||
* the OEIS entry: [https://oeis.org/wiki/Mersenne_primes Mersenne primes]. |
|||
* the OEIS entry: [https://oeis.org/A000043 A000043 Mersenne exponents: primes p such that 2^p - 1 is prime. Then 2^p - 1 is called a Mersenne prime]. |
|||
<br><br> |
|||
=={{header|11l}}== |
|||
{{trans|D}} |
|||
<syntaxhighlight lang="11l">F is_prime(BigInt bi) |
|||
I bi < 2 {R 0B} |
|||
I bi % 2 == 0 {R bi == 2} |
|||
I bi % 3 == 0 {R bi == 3} |
|||
V test = BigInt(5) |
|||
L test * test < bi |
|||
I bi % test == 0 |
|||
R 0B |
|||
test += 2 |
|||
I bi % test == 0 |
|||
R 0B |
|||
test += 4 |
|||
R 1B |
|||
V base = BigInt(2) |
|||
L(p) 1..31 |
|||
I is_prime(base - 1) |
|||
print(‘2 ^ ’p‘ - 1’) |
|||
base *= 2</syntaxhighlight> |
|||
{{out}} |
|||
<pre> |
|||
2 ^ 2 - 1 |
|||
2 ^ 3 - 1 |
|||
2 ^ 5 - 1 |
|||
2 ^ 7 - 1 |
|||
2 ^ 13 - 1 |
|||
2 ^ 17 - 1 |
|||
2 ^ 19 - 1 |
|||
2 ^ 31 - 1 |
|||
</pre> |
|||
=={{header|ALGOL 60}}== |
|||
{{works with|A60}} |
|||
<syntaxhighlight lang="ALGOL"> |
|||
begin |
|||
integer procedure mersenne(n); |
|||
value n; integer n; |
|||
begin |
|||
integer i, m; |
|||
m := 1; |
|||
for i := 1 step 1 until n do |
|||
m := m * 2; |
|||
mersenne := m - 1; |
|||
end; |
|||
boolean procedure isprime(n); |
|||
value n; integer n; |
|||
begin |
|||
if n < 2 then |
|||
isprime := false |
|||
else if entier(n / 2) * 2 = n then |
|||
isprime := (n = 2) |
|||
else |
|||
begin |
|||
comment - check odd divisors up to sqrt(n); |
|||
integer i, limit; |
|||
boolean divisible; |
|||
i := 3; |
|||
limit := entier(sqrt(n)); |
|||
divisible := false; |
|||
for i := i while i <= limit and not divisible do |
|||
begin |
|||
if entier(n / i) * i = n then |
|||
divisible := true; |
|||
i := i + 2 |
|||
end; |
|||
isprime := if divisible then false else true; |
|||
end; |
|||
end; |
|||
comment - main code begins here; |
|||
integer i, m; |
|||
outstring(1,"Searching to M(31) for Mersenne primes\n"); |
|||
for i := 1 step 1 until 31 do |
|||
begin |
|||
m := mersenne(i); |
|||
if isprime(m) then |
|||
begin |
|||
outstring(1,"M("); |
|||
outinteger(1,i); |
|||
outstring(1,") : "); |
|||
outinteger(1,m); |
|||
outstring(1,"\n"); |
|||
end; |
|||
end; |
|||
end |
|||
</syntaxhighlight> |
|||
{{out}} |
|||
<pre> |
|||
Searching to M(31) for Mersenne primes |
|||
M( 2 ) : 3 |
|||
M( 3 ) : 7 |
|||
M( 5 ) : 31 |
|||
M( 7 ) : 127 |
|||
M( 13 ) : 8191 |
|||
M( 17 ) : 131071 |
|||
M( 19 ) : 524287 |
|||
M( 31 ) : 2147483647 |
|||
</pre> |
|||
=={{header|ALGOL 68}}== |
|||
{{works with|ALGOL 68G|Any - tested with release 2.8.3.win32}} |
|||
<syntaxhighlight lang="algol68"> |
|||
BEGIN # find some Mersenne Primrs - primes of the form 2^n - 1, n must be prime # |
|||
# This assumes LONG INT is at least 64 bits (as in e.g. Algol 68G) # |
|||
# we handle 2 as a special case and then the odd numbers starting at 3 # |
|||
PR read "primes.incl.a68" PR # include prime utilities # |
|||
# 2^0 - 1 = 0 and 2^1 - 1 = 1, neither of which is prime # |
|||
# start from 2^2 # |
|||
INT n := 2; |
|||
LONG INT p2 := 4; |
|||
IF is probably prime( p2 - 1 ) THEN print( ( " ", whole( n, 0 ) ) ) FI; |
|||
# 2^3, 2^5, etc. # |
|||
n +:= 1; |
|||
p2 *:= 2; |
|||
WHILE n < 61 DO |
|||
IF is probably prime( p2 - 1 ) THEN print( ( " ", whole( n, 0 ) ) ) FI; |
|||
n +:= 2; |
|||
p2 *:= 4 |
|||
OD; |
|||
IF is probably prime( p2 - 1 ) THEN print( ( " ", whole( n, 0 ) ) ) FI |
|||
END |
|||
</syntaxhighlight> |
|||
{{out}} |
|||
<pre> |
|||
2 3 5 7 13 17 19 31 61 |
|||
</pre> |
|||
=={{header|ALGOL W}}== |
|||
<syntaxhighlight lang="algolw"> |
|||
begin % find some Mersenne Primrs - primes of the form 2^n - 1, n must be prime % |
|||
% integers are 32 bit in Algol W, so there won't be many numbers to check % |
|||
% we handle 2 as a special case and then the odd numbers starting at 3 % |
|||
% simplified primality by trial division test - no need to check for even % |
|||
% numbers as 2^n - 1 is always odd for n >= 2 % |
|||
logical procedure oddOnlyPrimalityTest ( integer value n ) ; |
|||
begin |
|||
logical isPrime; |
|||
isPrime := true; |
|||
for i := 3 step 2 until entier( sqrt( n ) ) do begin; |
|||
isPrime := n rem i not = 0; |
|||
if not isPrime then goto endPrimalityTest |
|||
end for_i; |
|||
endPrimalityTest: isPrime |
|||
end oddOnlyPrimalityTest ; |
|||
integer p2, n; |
|||
n := 1; |
|||
p2 := 2; |
|||
while |
|||
begin |
|||
if n < 3 then begin |
|||
n := n + 1; |
|||
p2 := p2 * 2 |
|||
end |
|||
else begin |
|||
n := n + 2; |
|||
p2 := p2 * 4 |
|||
end if_n_le_3__ ;; |
|||
if oddOnlyPrimalityTest( p2 - 1 ) then writeon( i_w := 1, s_w := 0, " ", n ); |
|||
n < 29 |
|||
end |
|||
do begin end ; |
|||
% MAXINTEGER is 2**31 - 1 % |
|||
if oddOnlyPrimalityTest( MAXINTEGER ) then writeon( i_w := 1, s_w := 0, " ", 31 ); |
|||
end. |
|||
</syntaxhighlight> |
|||
{{out}} |
|||
<pre> |
|||
2 3 5 7 13 17 19 31 |
|||
</pre> |
|||
=={{header|AppleScript}}== |
=={{header|AppleScript}}== |
||
<syntaxhighlight lang="applescript"> |
|||
<lang> |
|||
on isPrime(integ) |
on isPrime(integ) |
||
set isComposite to "" |
set isComposite to "" |
||
Line 43: | Line 238: | ||
set x to x + 1 |
set x to x + 1 |
||
end repeat |
end repeat |
||
</syntaxhighlight> |
|||
</lang> |
|||
---- |
|||
Or more efficiently: |
|||
<syntaxhighlight lang="applescript">on isPrime(n) |
|||
if (n < 4) then return (n > 1) |
|||
if ((n mod 2 is 0) or (n mod 3 is 0)) then return false |
|||
set {i, max} to {5, (n ^ 0.5) div 1} |
|||
repeat until (i > max) |
|||
if ((n mod i is 0) or (n mod (i + 2) is 0)) then return false |
|||
set i to i + 6 |
|||
end repeat |
|||
return true |
|||
end isPrime |
|||
on join(lst, delim) |
|||
set astid to AppleScript's text item delimiters |
|||
set AppleScript's text item delimiters to delim |
|||
set txt to lst as text |
|||
set AppleScript's text item delimiters to astid |
|||
return txt |
|||
end join |
|||
on mersennePrimes() |
|||
set output to {"Mersenne primes within AppleScript's number precision:"} |
|||
-- Special-case 2 ^ 2 - 1. |
|||
set end of output to "2 ^ 2 - 1 = 3" |
|||
set p to 1 -- Otherwise test odd-numbered powers of 2. |
|||
try -- Survive the "numeric operation too large" error when it occurs. |
|||
repeat |
|||
set p to p + 2 |
|||
if ((isPrime(p)) and (isPrime(2 ^ p - 1))) then ¬ |
|||
set end of output to "2 ^ " & p & " - 1 = " & (2 ^ p div 1 - 1) |
|||
end repeat |
|||
end try |
|||
return join(output, linefeed) |
|||
end mersennePrimes |
|||
mersennePrimes()</syntaxhighlight> |
|||
{{output}} |
|||
<syntaxhighlight lang="applescript">"Mersenne primes within AppleScript's number precision: |
|||
2 ^ 2 - 1 = 3 |
|||
2 ^ 3 - 1 = 7 |
|||
2 ^ 5 - 1 = 31 |
|||
2 ^ 7 - 1 = 127 |
|||
2 ^ 13 - 1 = 8191 |
|||
2 ^ 17 - 1 = 131071 |
|||
2 ^ 19 - 1 = 524287 |
|||
2 ^ 31 - 1 = 2.147483647E+9"</syntaxhighlight> |
|||
=={{header|Arturo}}== |
|||
<syntaxhighlight lang="arturo">mersenne?: function [n][ |
|||
prime? dec 2 ^ n |
|||
] |
|||
1..31 | select => mersenne? |
|||
| loop 'x -> print ["M (" x ") = 2 ^" x " - 1 =" (2^x)-1 "is a prime"]</syntaxhighlight> |
|||
{{out}} |
|||
<pre>M ( 2 ) = 2 ^ 2 - 1 = 3 is a prime |
|||
M ( 3 ) = 2 ^ 3 - 1 = 7 is a prime |
|||
M ( 5 ) = 2 ^ 5 - 1 = 31 is a prime |
|||
M ( 7 ) = 2 ^ 7 - 1 = 127 is a prime |
|||
M ( 13 ) = 2 ^ 13 - 1 = 8191 is a prime |
|||
M ( 17 ) = 2 ^ 17 - 1 = 131071 is a prime |
|||
M ( 19 ) = 2 ^ 19 - 1 = 524287 is a prime |
|||
M ( 31 ) = 2 ^ 31 - 1 = 2147483647 is a prime</pre> |
|||
=={{header|AWK}}== |
=={{header|AWK}}== |
||
<syntaxhighlight lang="awk"> |
|||
<lang AWK> |
|||
# syntax: GAWK --bignum -f MERSENNE_PRIMES.AWK |
# syntax: GAWK --bignum -f MERSENNE_PRIMES.AWK |
||
BEGIN { |
BEGIN { |
||
Line 71: | Line 334: | ||
return(1) |
return(1) |
||
} |
} |
||
</syntaxhighlight> |
|||
</lang> |
|||
{{out}} |
{{out}} |
||
<pre> |
<pre> |
||
Line 84: | Line 347: | ||
2 ^ 61 - 1 |
2 ^ 61 - 1 |
||
</pre> |
</pre> |
||
=={{header|BASIC}}== |
|||
==={{header|BASIC256}}=== |
|||
{{trans|FreeBASIC}} |
|||
<syntaxhighlight lang="basic256">m = 0 |
|||
while True |
|||
m += 1 |
|||
if isPrime((2^m)-1) = True then print m |
|||
end while |
|||
end |
|||
function isPrime(v) |
|||
if v <= 1 then return False |
|||
for i = 2 To int(sqr(v)) |
|||
if v % i = 0 then return False |
|||
next i |
|||
return True |
|||
end function</syntaxhighlight> |
|||
{{out}} |
|||
<pre> |
|||
Igual que la entrada de FreeBASIC. |
|||
</pre> |
|||
==={{header|FreeBASIC}}=== |
|||
<syntaxhighlight lang="freebasic">Function isPrime(Byval ValorEval As Integer) As Boolean |
|||
If ValorEval <= 1 Then Return False |
|||
For i As Integer = 2 To Int(Sqr(ValorEval)) |
|||
If ValorEval Mod i = 0 Then Return False |
|||
Next i |
|||
Return True |
|||
End Function |
|||
Dim As Integer m = 0 |
|||
While True |
|||
m += 1 |
|||
If isPrime((2^m)-1) Then Print m |
|||
Wend |
|||
Sleep</syntaxhighlight> |
|||
{{out}} |
|||
<pre>2 |
|||
3 |
|||
5 |
|||
7 |
|||
13 |
|||
17 |
|||
19 |
|||
31</pre> |
|||
==={{header|Yabasic}}=== |
|||
{{trans|FreeBASIC}} |
|||
<syntaxhighlight lang="yabasic">m = 0 |
|||
while True |
|||
m = m + 1 |
|||
if isPrime((2^m)-1) = True then print m : fi |
|||
wend |
|||
end |
|||
sub isPrime(v) |
|||
if v < 2 then return False : fi |
|||
if mod(v, 2) = 0 then return v = 2 : fi |
|||
if mod(v, 3) = 0 then return v = 3 : fi |
|||
d = 5 |
|||
while d * d <= v |
|||
if mod(v, d) = 0 then return False else d = d + 2 : fi |
|||
end while |
|||
return True |
|||
end sub</syntaxhighlight> |
|||
{{out}} |
|||
<pre> |
|||
Igual que la entrada de FreeBASIC. |
|||
</pre> |
|||
=={{header|C}}== |
=={{header|C}}== |
||
< |
<syntaxhighlight lang="c">#include <stdbool.h> |
||
#include <stdint.h> |
#include <stdint.h> |
||
#include <stdio.h> |
#include <stdio.h> |
||
Line 120: | Line 455: | ||
return 0; |
return 0; |
||
}</ |
}</syntaxhighlight> |
||
{{out}} |
{{out}} |
||
<pre>2 ^ 2 - 1 |
<pre>2 ^ 2 - 1 |
||
Line 130: | Line 465: | ||
2 ^ 19 - 1 |
2 ^ 19 - 1 |
||
2 ^ 31 - 1</pre> |
2 ^ 31 - 1</pre> |
||
{{libheader|GMP}} |
|||
Alternatively, we can use GMP to find the first 23 Mersenne primes in about the same time as the corresponding Wren entry. |
|||
<syntaxhighlight lang="c">#include <stdio.h> |
|||
#include <stdbool.h> |
|||
#include <stdint.h> |
|||
#include <gmp.h> |
|||
#define MAX 23 |
|||
bool isPrime(uint64_t n) { |
|||
uint64_t test; |
|||
if (n < 2) return false; |
|||
if (n % 2 == 0) return n == 2; |
|||
if (n % 3 == 0) return n == 3; |
|||
test = 5; |
|||
while (test * test < n) { |
|||
if (n % test == 0) return false; |
|||
test += 2; |
|||
if (n % test == 0) return false; |
|||
test += 4; |
|||
} |
|||
return true; |
|||
} |
|||
int main() { |
|||
uint64_t p = 2; |
|||
int count = 0; |
|||
mpz_t m, one; |
|||
mpz_init(m); |
|||
mpz_init_set_ui(one, 1); |
|||
while (true) { |
|||
mpz_mul_2exp(m, one, p); |
|||
mpz_sub_ui(m, m, 1); |
|||
if (mpz_probab_prime_p(m, 15) > 0) { |
|||
printf("2 ^ %ld - 1\n", p); |
|||
if (++count == MAX) break; |
|||
} |
|||
while (true) { |
|||
p = (p > 2) ? p + 2 : 3; |
|||
if (isPrime(p)) break; |
|||
} |
|||
} |
|||
mpz_clear(m); |
|||
mpz_clear(one); |
|||
return 0; |
|||
} </syntaxhighlight> |
|||
{{out}} |
|||
<pre> |
|||
Same as Wren example. |
|||
</pre> |
|||
=={{header|C++}}== |
=={{header|C++}}== |
||
{{trans|C}} |
{{trans|C}} |
||
< |
<syntaxhighlight lang="cpp">#include <iostream> |
||
bool isPrime(uint64_t n) { |
bool isPrime(uint64_t n) { |
||
Line 162: | Line 550: | ||
return 0; |
return 0; |
||
}</ |
}</syntaxhighlight> |
||
{{out}} |
{{out}} |
||
<pre>2 ^ 2 - 1 |
<pre>2 ^ 2 - 1 |
||
Line 175: | Line 563: | ||
=={{header|C sharp|C#}}== |
=={{header|C sharp|C#}}== |
||
Needs a better primality checking algorithm to do really large prime numbers. |
Needs a better primality checking algorithm to do really large prime numbers. |
||
< |
<syntaxhighlight lang="csharp">using System; |
||
using System.Numerics; |
using System.Numerics; |
||
Line 233: | Line 621: | ||
} |
} |
||
} |
} |
||
}</ |
}</syntaxhighlight> |
||
{{out}} |
{{out}} |
||
<pre>2 ^ 2 - 1 |
<pre>2 ^ 2 - 1 |
||
Line 247: | Line 635: | ||
=={{header|D}}== |
=={{header|D}}== |
||
Simplest thing that could possibly work. Using better primality tests will allow for more results to be calculated in a reasonable amount of time. |
Simplest thing that could possibly work. Using better primality tests will allow for more results to be calculated in a reasonable amount of time. |
||
< |
<syntaxhighlight lang="d">import std.bigint; |
||
import std.stdio; |
import std.stdio; |
||
Line 275: | Line 663: | ||
base *= 2; |
base *= 2; |
||
} |
} |
||
}</ |
}</syntaxhighlight> |
||
{{out}} |
{{out}} |
||
<pre>2 ^ 2 - 1 |
<pre>2 ^ 2 - 1 |
||
Line 286: | Line 674: | ||
2 ^ 31 - 1</pre> |
2 ^ 31 - 1</pre> |
||
=={{header| |
=={{header|Delphi}}== |
||
{{works with|Delphi|6.0}} |
|||
{{libheader|SysUtils,StdCtrls}} |
|||
<syntaxhighlight lang="Delphi"> |
|||
function IsPrime(N: int64): boolean; |
|||
{Fast, optimised prime test} |
|||
var I,Stop: int64; |
|||
begin |
|||
if (N = 2) or (N=3) then Result:=true |
|||
else if (n <= 1) or ((n mod 2) = 0) or ((n mod 3) = 0) then Result:= false |
|||
else |
|||
begin |
|||
I:=5; |
|||
Stop:=Trunc(sqrt(N+0.0)); |
|||
Result:=False; |
|||
while I<=Stop do |
|||
begin |
|||
if ((N mod I) = 0) or ((N mod (I + 2)) = 0) then exit; |
|||
Inc(I,6); |
|||
end; |
|||
Result:=True; |
|||
end; |
|||
end; |
|||
procedure MersennePrimes(Memo: TMemo); |
|||
var N: integer; |
|||
var Mn: int64; |
|||
begin |
|||
Memo.Lines.Add('2^N-1 Prime'); |
|||
Memo.Lines.Add('--------------------' ); |
|||
for N:=1 to 32 do |
|||
begin |
|||
Mn:=(1 shl N)-1; |
|||
if IsPrime(Mn) then |
|||
Memo.Lines.Add(Format('2^%d-1 %14.0n',[N,Mn+0.0])); |
|||
end; |
|||
end; |
|||
</syntaxhighlight> |
|||
{{out}} |
|||
<pre> |
|||
2^N-1 Prime |
|||
-------------------- |
|||
2^2-1 3 |
|||
2^3-1 7 |
|||
2^5-1 31 |
|||
2^7-1 127 |
|||
2^13-1 8,191 |
|||
2^17-1 131,071 |
|||
2^19-1 524,287 |
|||
2^31-1 2,147,483,647 |
|||
</pre> |
|||
=={{header|EasyLang}}== |
|||
<syntaxhighlight> |
|||
fastfunc isprim num . |
|||
if num mod 2 = 0 |
|||
if num = 2 |
|||
return 1 |
|||
. |
|||
return 0 |
|||
. |
|||
i = 3 |
|||
while i <= sqrt num |
|||
if num mod i = 0 |
|||
return 0 |
|||
. |
|||
i += 2 |
|||
. |
|||
return 1 |
|||
. |
|||
base = 4 |
|||
for p = 2 to 52 |
|||
if isprim (base - 1) = 1 |
|||
print "2 ^ " & p & " - 1" |
|||
. |
|||
base *= 2 |
|||
. |
|||
</syntaxhighlight> |
|||
{{out}} |
|||
<pre> |
|||
2 ^ 2 - 1 |
|||
2 ^ 3 - 1 |
|||
2 ^ 5 - 1 |
|||
2 ^ 7 - 1 |
|||
2 ^ 13 - 1 |
|||
2 ^ 17 - 1 |
|||
2 ^ 19 - 1 |
|||
2 ^ 31 - 1 |
|||
</pre> |
|||
=={{header|F_Sharp|F#}}== |
|||
{{trans|C#}} |
{{trans|C#}} |
||
< |
<syntaxhighlight lang="fsharp">open System |
||
open System.Numerics |
open System.Numerics |
||
Line 341: | Line 826: | ||
loop 2 0 |> ignore |
loop 2 0 |> ignore |
||
0 // return an integer exit code</ |
0 // return an integer exit code</syntaxhighlight> |
||
{{out}} |
{{out}} |
||
<pre>2 ^ 2 - 1 |
<pre>2 ^ 2 - 1 |
||
Line 355: | Line 840: | ||
=={{header|Factor}}== |
=={{header|Factor}}== |
||
Factor comes with a Lucas-Lehmer primality test. |
Factor comes with a Lucas-Lehmer primality test. |
||
< |
<syntaxhighlight lang="factor">USING: formatting math.primes.lucas-lehmer math.ranges sequences ; |
||
: mersennes-upto ( n -- seq ) [1,b] [ lucas-lehmer ] filter ; |
: mersennes-upto ( n -- seq ) [1,b] [ lucas-lehmer ] filter ; |
||
3500 mersennes-upto [ "2 ^ %d - 1\n" printf ] each</ |
3500 mersennes-upto [ "2 ^ %d - 1\n" printf ] each</syntaxhighlight> |
||
{{out}} |
{{out}} |
||
<pre> |
<pre> |
||
Line 380: | Line 865: | ||
2 ^ 2281 - 1 |
2 ^ 2281 - 1 |
||
2 ^ 3217 - 1 |
2 ^ 3217 - 1 |
||
</pre> |
|||
=={{header|Fortran}}== |
|||
{{Trans|C}} |
|||
<syntaxhighlight lang="fortran"> |
|||
program mersenne |
|||
use iso_fortran_env, only: output_unit, INT64 |
|||
implicit none |
|||
integer, parameter :: l=INT64 |
|||
integer(kind=l) :: base |
|||
integer :: pow |
|||
base = 2 |
|||
do pow = 1, 32 |
|||
if (is_prime(base-1)) then |
|||
write(output_unit,'(A2,x,I0,x,A3)') "2^", pow, "- 1" |
|||
end if |
|||
base = base * 2 |
|||
end do |
|||
contains |
|||
pure function is_prime(n) |
|||
integer(kind=l), intent(in) :: n |
|||
logical :: is_prime |
|||
integer(kind=l) :: test |
|||
is_prime = .false. |
|||
if (n < 2) return |
|||
if (modulo(n, 2) == 0) then |
|||
is_prime = n==2 |
|||
return |
|||
end if |
|||
if (modulo(n, 3) == 0) then |
|||
is_prime = n==3 |
|||
return |
|||
end if |
|||
test = 5 |
|||
do |
|||
if (test**2 >= n) then |
|||
is_prime = .true. |
|||
return |
|||
end if |
|||
if (modulo(n, test) == 0) return |
|||
test = test + 2 |
|||
if (modulo(n, test) == 0) return |
|||
test = test + 4 |
|||
end do |
|||
end function is_prime |
|||
end program mersenne |
|||
</syntaxhighlight> |
|||
{{out}} |
|||
<pre> |
|||
2^ 2 - 1 |
|||
2^ 3 - 1 |
|||
2^ 5 - 1 |
|||
2^ 7 - 1 |
|||
2^ 13 - 1 |
|||
2^ 17 - 1 |
|||
2^ 19 - 1 |
|||
2^ 31 - 1 |
|||
</pre> |
|||
=={{header|Frink}}== |
|||
Frink has built-in routines for iterating through prime numbers. Frink's <CODE>isPrime[n]</CODE> function automatically detects numbers of the form 2<sup>n</sup>-1 and performs a more efficient Lucas-Lehmer primality test on the number. This works with arbitrarily large numbers. |
|||
Did you know that all Java-based JVM languages are many many orders of magnitude faster because Frink's developer contributed vastly faster BigInteger algorithms to Java? It took the Java developers 11 years to integrate them but they became part of 1.8 and later! Operations that used to take days now can be done in seconds thanks to Frink's contributions to Java. |
|||
<syntaxhighlight lang="frink">for n = primes[] |
|||
if isPrime[2^n - 1] |
|||
println["2^$n - 1"]</syntaxhighlight> |
|||
{{out}} |
|||
<pre> |
|||
2^2 - 1 |
|||
2^3 - 1 |
|||
2^5 - 1 |
|||
2^7 - 1 |
|||
2^13 - 1 |
|||
2^17 - 1 |
|||
2^19 - 1 |
|||
2^31 - 1 |
|||
2^61 - 1 |
|||
2^89 - 1 |
|||
2^107 - 1 |
|||
2^127 - 1 |
|||
2^521 - 1 |
|||
2^607 - 1 |
|||
2^1279 - 1 |
|||
2^2203 - 1 |
|||
2^2281 - 1 |
|||
2^3217 - 1 |
|||
2^4253 - 1 |
|||
2^4423 - 1 |
|||
2^9689 - 1 |
|||
2^9941 - 1 |
|||
2^11213 - 1 |
|||
2^19937 - 1 |
|||
2^21701 - 1 |
|||
... |
|||
</pre> |
</pre> |
||
Line 391: | Line 977: | ||
Note that the use of ProbablyPrime(0) requires Go 1.8 or later. When using the <code>math/big</code> package, passing a parameter of zero to this method forces it to apply only the Baillie-PSW test to check for primality. This is 100% accurate for numbers up to 2^64 and at the time of writing (June 2018) no known composite number above that bound passes the test. |
Note that the use of ProbablyPrime(0) requires Go 1.8 or later. When using the <code>math/big</code> package, passing a parameter of zero to this method forces it to apply only the Baillie-PSW test to check for primality. This is 100% accurate for numbers up to 2^64 and at the time of writing (June 2018) no known composite number above that bound passes the test. |
||
< |
<syntaxhighlight lang="go">package main |
||
import ( |
import ( |
||
Line 432: | Line 1,018: | ||
} |
} |
||
} |
} |
||
}</ |
}</syntaxhighlight> |
||
{{out|text=using the GMP package on a 3.4 GHz Xeon E3-1245:}} |
{{out|text=using the GMP package on a 3.4 GHz Xeon E3-1245:}} |
||
Line 462: | Line 1,048: | ||
This can be sped up quite a bit for modern multi-core CPUs by some simple changes to use goroutines. |
This can be sped up quite a bit for modern multi-core CPUs by some simple changes to use goroutines. |
||
< |
<syntaxhighlight lang="go">package main |
||
import ( |
import ( |
||
Line 527: | Line 1,113: | ||
} |
} |
||
} |
} |
||
}</ |
}</syntaxhighlight> |
||
<!-- The output and the previous should be done on the same hardware to keep the numbers comparible; if changing/updating one change the other too. --> |
<!-- The output and the previous should be done on the same hardware to keep the numbers comparible; if changing/updating one change the other too. --> |
||
{{out|text=using the GMP package on the same 3.4 GHz Xeon E3-1245 (4 core × 2 SMT threads) as above:}} |
{{out|text=using the GMP package on the same 3.4 GHz Xeon E3-1245 (4 core × 2 SMT threads) as above:}} |
||
Line 558: | Line 1,144: | ||
Using this approach, the Celeron machine (dual core) takes ~180 seconds to reach M<sub>9941</sub> and ~270 seconds to reach M<sub>11213</sub>. |
Using this approach, the Celeron machine (dual core) takes ~180 seconds to reach M<sub>9941</sub> and ~270 seconds to reach M<sub>11213</sub>. |
||
=={{header|Haskell}}== |
=={{header|Haskell}}== |
||
< |
<syntaxhighlight lang="haskell">import Data.Numbers.Primes (primes) |
||
import Data.Numbers.Primes (primes) |
|||
import Text.Printf (printf) |
import Text.Printf (printf) |
||
lucasLehmer :: Int -> Bool |
lucasLehmer :: Int -> Bool |
||
lucasLehmer p = |
lucasLehmer p = iterate f 4 !! p-2 == 0 |
||
where |
where |
||
f b = (b^2 - 2) `mod` m |
|||
m = 2^p - 1 |
m = 2^p - 1 |
||
main = mapM_ (printf "M %d\n") $ take 20 mersenne |
main = mapM_ (printf "M %d\n") $ take 20 mersenne |
||
where |
where |
||
mersenne = filter lucasLehmer primes</ |
mersenne = filter lucasLehmer primes</syntaxhighlight> |
||
{{out}} |
{{out}} |
||
<pre> |
<pre> |
||
Line 594: | Line 1,179: | ||
M 9689 |
M 9689 |
||
</pre> |
</pre> |
||
=={{header|J}}== |
|||
<syntaxhighlight lang="j"> I. 1 p: <: 2x ^ i. 1280 |
|||
2 3 5 7 13 17 19 31 61 89 107 127 521 607 1279</syntaxhighlight> |
|||
=={{header|Java}}== |
=={{header|Java}}== |
||
{{trans|Kotlin}} |
{{trans|Kotlin}} |
||
< |
<syntaxhighlight lang="java">import java.math.BigInteger; |
||
public class MersennePrimes { |
public class MersennePrimes { |
||
Line 634: | Line 1,223: | ||
} |
} |
||
} |
} |
||
}</ |
}</syntaxhighlight> |
||
{{out}} |
{{out}} |
||
<pre>2 ^ 2 - 1 |
<pre>2 ^ 2 - 1 |
||
Line 662: | Line 1,251: | ||
Julia module <code>Primes</code> uses Miller-Rabin primality test. |
Julia module <code>Primes</code> uses Miller-Rabin primality test. |
||
< |
<syntaxhighlight lang="julia">using Primes |
||
mersenne(n::Integer) = convert(typeof(n), 2) ^ n - one(n) |
mersenne(n::Integer) = convert(typeof(n), 2) ^ n - one(n) |
||
Line 676: | Line 1,265: | ||
end |
end |
||
main(big(20))</ |
main(big(20))</syntaxhighlight> |
||
{{out}} |
{{out}} |
||
Line 705: | Line 1,294: | ||
A 'certainty' parameter of 10 is enough to find the first 20 Mersenne primes but as even this takes about 90 seconds on my modest machine I've not bothered going beyond that. |
A 'certainty' parameter of 10 is enough to find the first 20 Mersenne primes but as even this takes about 90 seconds on my modest machine I've not bothered going beyond that. |
||
< |
<syntaxhighlight lang="scala">// version 1.2.10 |
||
import java.math.BigInteger |
import java.math.BigInteger |
||
Line 744: | Line 1,333: | ||
} |
} |
||
} |
} |
||
}</ |
}</syntaxhighlight> |
||
{{out}} |
{{out}} |
||
Line 772: | Line 1,361: | ||
=={{header|Lua}}== |
=={{header|Lua}}== |
||
This checks for primality using a trial division function. The limitation is 'until p == p + 1', meaning that the program will halt when Lua's number type (a 64-bit floating point value) no longer has enough precision to distiguish between one integer and the next. |
This checks for primality using a trial division function. The limitation is 'until p == p + 1', meaning that the program will halt when Lua's number type (a 64-bit floating point value) no longer has enough precision to distiguish between one integer and the next. |
||
< |
<syntaxhighlight lang="lua">-- Returns a boolean to show whether x is prime |
||
function isPrime (x) |
function isPrime (x) |
||
if x < 2 then return false end |
if x < 2 then return false end |
||
Line 791: | Line 1,380: | ||
print("2 ^ " .. i .. " - 1 = " .. p) |
print("2 ^ " .. i .. " - 1 = " .. p) |
||
end |
end |
||
until p == p + 1</ |
until p == p + 1</syntaxhighlight> |
||
{{out}} |
{{out}} |
||
<pre>2 ^ 2 - 1 = 3 |
<pre>2 ^ 2 - 1 = 3 |
||
Line 801: | Line 1,390: | ||
2 ^ 19 - 1 = 524287 |
2 ^ 19 - 1 = 524287 |
||
2 ^ 31 - 1 = 2147483647</pre> |
2 ^ 31 - 1 = 2147483647</pre> |
||
=={{header|Mathematica}}/{{header|Wolfram Language}}== |
|||
Mathematica has a built-in function to show the Mersenne primes: |
|||
<syntaxhighlight lang="mathematica">MersennePrimeExponent /@ Range[40]</syntaxhighlight> |
|||
{{out}} |
|||
<pre>{2, 3, 5, 7, 13, 17, 19, 31, 61, 89, 107, 127, 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}</pre> |
|||
Alternatively, we can calculate them: |
|||
<syntaxhighlight lang="mathematica">res = {}; |
|||
Do[ |
|||
If[PrimeQ[2^i - 1], |
|||
AppendTo[res, i]; |
|||
If[Length[res] == 20, |
|||
Break[] |
|||
] |
|||
] |
|||
, |
|||
{i, 1, \[Infinity]} |
|||
] |
|||
res</syntaxhighlight> |
|||
{{out}} |
|||
<pre>{2, 3, 5, 7, 13, 17, 19, 31, 61, 89, 107, 127, 521, 607, 1279, 2203, 2281, 3217, 4253, 4423}</pre> |
|||
=={{header|Nim}}== |
|||
===Using only standard library=== |
|||
If we want to use only the standard library, we are limited to 64 bits. So we used a simple primality test. |
|||
<syntaxhighlight lang="nim">func isOddPrime(n: uint64): bool = |
|||
if n == 1: return false |
|||
if n mod 3 == 0: return n == 3 |
|||
var d = 5u |
|||
while d * d <= n: |
|||
if n mod d == 0: return false |
|||
inc d, 2 |
|||
if n mod d == 0: return false |
|||
inc d, 4 |
|||
result = true |
|||
var p = 2u64 |
|||
for e in 1..63: |
|||
if isOddPrime(p - 1): |
|||
echo "2^", e, " - 1" |
|||
p *= 2</syntaxhighlight> |
|||
{{out}} |
|||
<pre>2^2 - 1 |
|||
2^3 - 1 |
|||
2^5 - 1 |
|||
2^7 - 1 |
|||
2^13 - 1 |
|||
2^17 - 1 |
|||
2^19 - 1 |
|||
2^31 - 1 |
|||
2^61 - 1</pre> |
|||
===Using big integers=== |
|||
{{libheader|bignum}} |
|||
The module <code>bignum</code> provides big integers and a probabilistic primality test. We searched the Mersenne numbers for exponents between 1 and 10_000. |
|||
<syntaxhighlight lang="nim">import bignum |
|||
var p = newInt(2) |
|||
for e in 1..10_000: |
|||
if probablyPrime(p - 1, 25) != 0: |
|||
echo "2^", e, " - 1" |
|||
p *= 2</syntaxhighlight> |
|||
{{out}} |
|||
<pre>2^2 - 1 |
|||
2^3 - 1 |
|||
2^5 - 1 |
|||
2^7 - 1 |
|||
2^13 - 1 |
|||
2^17 - 1 |
|||
2^19 - 1 |
|||
2^31 - 1 |
|||
2^61 - 1 |
|||
2^89 - 1 |
|||
2^107 - 1 |
|||
2^127 - 1 |
|||
2^521 - 1 |
|||
2^607 - 1 |
|||
2^1279 - 1 |
|||
2^2203 - 1 |
|||
2^2281 - 1 |
|||
2^3217 - 1 |
|||
2^4253 - 1 |
|||
2^4423 - 1 |
|||
2^9689 - 1 |
|||
2^9941 - 1</pre> |
|||
=={{header|PARI/GP}}== |
=={{header|PARI/GP}}== |
||
< |
<syntaxhighlight lang="parigp">LL(p)={ |
||
my(m=Mod(4,1<<p-1)); |
my(m=Mod(4,1<<p-1)); |
||
for(i=3,p,m=m^2-2); |
for(i=3,p,m=m^2-2); |
||
m==0 |
m==0 |
||
}; |
}; |
||
forprime(p=2,, if(LL(p), print("2^"p"-1")))</ |
forprime(p=2,, if(LL(p), print("2^"p"-1")))</syntaxhighlight> |
||
=={{header|Perl}}== |
=={{header|Perl}}== |
||
Line 816: | Line 1,493: | ||
{{libheader|ntheory}} |
{{libheader|ntheory}} |
||
< |
<syntaxhighlight lang="perl">use ntheory qw/forprimes is_mersenne_prime/; |
||
forprimes { is_mersenne_prime($_) && say } 1e9;</ |
forprimes { is_mersenne_prime($_) && say } 1e9;</syntaxhighlight> |
||
{{out}} |
{{out}} |
||
<pre> |
<pre> |
||
Line 833: | Line 1,510: | ||
=={{header|Phix}}== |
=={{header|Phix}}== |
||
{{libheader|mpfr}} |
{{libheader|Phix/mpfr}} |
||
<!--<syntaxhighlight lang="phix">(phixonline)--> |
|||
<lang Phix>include mpfr.e |
|||
<span style="color: #008080;">with</span> <span style="color: #008080;">javascript_semantics</span> |
|||
atom t0 = time() |
|||
<span style="color: #008080;">include</span> <span style="color: #004080;">mpfr</span><span style="color: #0000FF;">.</span><span style="color: #000000;">e</span> |
|||
mpz mp = mpz_init(), |
|||
<span style="color: #004080;">atom</span> <span style="color: #000000;">t0</span> <span style="color: #0000FF;">=</span> <span style="color: #7060A8;">time</span><span style="color: #0000FF;">()</span> |
|||
bp = mpz_init() |
|||
<span style="color: #004080;">mpz</span> <span style="color: #000000;">mp</span> <span style="color: #0000FF;">=</span> <span style="color: #7060A8;">mpz_init</span><span style="color: #0000FF;">(),</span> |
|||
randstate state = gmp_randinit_mt() |
|||
<span style="color: #000000;">bp</span> <span style="color: #0000FF;">=</span> <span style="color: #7060A8;">mpz_init</span><span style="color: #0000FF;">()</span> |
|||
integer p = 0, count = 0 |
|||
<span style="color: #004080;">integer</span> <span style="color: #000000;">p</span> <span style="color: #0000FF;">=</span> <span style="color: #000000;">0</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">count</span> <span style="color: #0000FF;">=</span> <span style="color: #000000;">0</span> |
|||
while true do |
|||
<span style="color: #008080;">constant</span> <span style="color: #000000;">lim</span> <span style="color: #0000FF;">=</span> <span style="color: #008080;">iff</span><span style="color: #0000FF;">(</span><span style="color: #7060A8;">platform</span><span style="color: #0000FF;">()=</span><span style="color: #004600;">JS</span><span style="color: #0000FF;">?</span><span style="color: #000000;">14</span><span style="color: #0000FF;">:</span><span style="color: #000000;">17</span><span style="color: #0000FF;">)</span> |
|||
mpz_ui_pow_ui(mp,2,p) |
|||
<span style="color: #008080;">while</span> <span style="color: #004600;">true</span> <span style="color: #008080;">do</span> |
|||
mpz_sub_ui(mp,mp,1) |
|||
<span style="color: #7060A8;">mpz_ui_pow_ui</span><span style="color: #0000FF;">(</span><span style="color: #000000;">mp</span><span style="color: #0000FF;">,</span><span style="color: #000000;">2</span><span style="color: #0000FF;">,</span><span style="color: #000000;">p</span><span style="color: #0000FF;">)</span> |
|||
if mpz_probable_prime_p(mp, state) then |
|||
<span style="color: #7060A8;">mpz_sub_ui</span><span style="color: #0000FF;">(</span><span style="color: #000000;">mp</span><span style="color: #0000FF;">,</span><span style="color: #000000;">mp</span><span style="color: #0000FF;">,</span><span style="color: #000000;">1</span><span style="color: #0000FF;">)</span> |
|||
string s = iff(time()-t0<1?"":", "&elapsed(time()-t0)) |
|||
<span style="color: #008080;">if</span> <span style="color: #7060A8;">mpz_prime</span><span style="color: #0000FF;">(</span><span style="color: #000000;">mp</span><span style="color: #0000FF;">)</span> <span style="color: #008080;">then</span> |
|||
printf(1, "2^%d-1%s\n",{p,s}) |
|||
<span style="color: #004080;">string</span> <span style="color: #000000;">s</span> <span style="color: #0000FF;">=</span> <span style="color: #008080;">iff</span><span style="color: #0000FF;">(</span><span style="color: #7060A8;">time</span><span style="color: #0000FF;">()-</span><span style="color: #000000;">t0</span><span style="color: #0000FF;"><</span><span style="color: #000000;">1</span><span style="color: #0000FF;">?</span><span style="color: #008000;">""</span><span style="color: #0000FF;">:</span><span style="color: #008000;">", "</span><span style="color: #0000FF;">&</span><span style="color: #7060A8;">elapsed</span><span style="color: #0000FF;">(</span><span style="color: #7060A8;">time</span><span style="color: #0000FF;">()-</span><span style="color: #000000;">t0</span><span style="color: #0000FF;">))</span> |
|||
count += 1 |
|||
<span style="color: #7060A8;">printf</span><span style="color: #0000FF;">(</span><span style="color: #000000;">1</span><span style="color: #0000FF;">,</span> <span style="color: #008000;">"2^%d-1%s\n"</span><span style="color: #0000FF;">,{</span><span style="color: #000000;">p</span><span style="color: #0000FF;">,</span><span style="color: #000000;">s</span><span style="color: #0000FF;">})</span> |
|||
if count>=17 then exit end if |
|||
<span style="color: #000000;">count</span> <span style="color: #0000FF;">+=</span> <span style="color: #000000;">1</span> |
|||
end if |
|||
<span style="color: #008080;">if</span> <span style="color: #000000;">count</span><span style="color: #0000FF;">>=</span><span style="color: #000000;">lim</span> <span style="color: #008080;">then</span> <span style="color: #008080;">exit</span> <span style="color: #008080;">end</span> <span style="color: #008080;">if</span> |
|||
while true do |
|||
<span style="color: #008080;">end</span> <span style="color: #008080;">if</span> |
|||
p = iff(p>2?p+2:3) |
|||
<span style="color: #008080;">while</span> <span style="color: #004600;">true</span> <span style="color: #008080;">do</span> |
|||
mpz_set_si(bp,p) |
|||
<span style="color: #000000;">p</span> <span style="color: #0000FF;">=</span> <span style="color: #008080;">iff</span><span style="color: #0000FF;">(</span><span style="color: #000000;">p</span><span style="color: #0000FF;">></span><span style="color: #000000;">2</span><span style="color: #0000FF;">?</span><span style="color: #000000;">p</span><span style="color: #0000FF;">+</span><span style="color: #000000;">2</span><span style="color: #0000FF;">:</span><span style="color: #000000;">3</span><span style="color: #0000FF;">)</span> |
|||
if mpz_probable_prime_p(bp, state) then exit end if |
|||
<span style="color: #7060A8;">mpz_set_si</span><span style="color: #0000FF;">(</span><span style="color: #000000;">bp</span><span style="color: #0000FF;">,</span><span style="color: #000000;">p</span><span style="color: #0000FF;">)</span> |
|||
end while |
|||
<span style="color: #008080;">if</span> <span style="color: #7060A8;">mpz_prime</span><span style="color: #0000FF;">(</span><span style="color: #000000;">bp</span><span style="color: #0000FF;">)</span> <span style="color: #008080;">then</span> <span style="color: #008080;">exit</span> <span style="color: #008080;">end</span> <span style="color: #008080;">if</span> |
|||
end while |
|||
<span style="color: #008080;">end</span> <span style="color: #008080;">while</span> |
|||
{mp,bp} = mpz_free({mp,bp}) |
|||
<span style="color: #008080;">end</span> <span style="color: #008080;">while</span> |
|||
state = gmp_randclear(state)</lang> |
|||
<span style="color: #0000FF;">{</span><span style="color: #000000;">mp</span><span style="color: #0000FF;">,</span><span style="color: #000000;">bp</span><span style="color: #0000FF;">}</span> <span style="color: #0000FF;">=</span> <span style="color: #7060A8;">mpz_free</span><span style="color: #0000FF;">({</span><span style="color: #000000;">mp</span><span style="color: #0000FF;">,</span><span style="color: #000000;">bp</span><span style="color: #0000FF;">})</span> |
|||
<!--</syntaxhighlight>--> |
|||
{{out}} |
{{out}} |
||
<pre> |
<pre> |
||
Line 880: | Line 1,559: | ||
=={{header|PHP}}== |
=={{header|PHP}}== |
||
< |
<syntaxhighlight lang="php"><?php |
||
function is_prime($n) { |
function is_prime($n) { |
||
Line 908: | Line 1,587: | ||
echo '2 ^ ', $i, ' - 1', PHP_EOL; |
echo '2 ^ ', $i, ' - 1', PHP_EOL; |
||
} |
} |
||
}</ |
}</syntaxhighlight> |
||
{{out}} |
{{out}} |
||
Line 920: | Line 1,599: | ||
2 ^ 31 - 1 |
2 ^ 31 - 1 |
||
2 ^ 61 - 1</pre> |
2 ^ 61 - 1</pre> |
||
=={{header|PicoLisp}}== |
|||
<syntaxhighlight lang="picolisp">(de **Mod (X Y N) |
|||
(let M 1 |
|||
(loop |
|||
(when (bit? 1 Y) |
|||
(setq M (% (* M X) N)) ) |
|||
(T (=0 (setq Y (>> 1 Y))) |
|||
M ) |
|||
(setq X (% (* X X) N)) ) ) ) |
|||
(de isprime (N) |
|||
(cache '(NIL) N |
|||
(if (== N 2) |
|||
T |
|||
(and |
|||
(> N 1) |
|||
(bit? 1 N) |
|||
(let (Q (dec N) N1 (dec N) K 0 X) |
|||
(until (bit? 1 Q) |
|||
(setq |
|||
Q (>> 1 Q) |
|||
K (inc K) ) ) |
|||
(catch 'composite |
|||
(do 16 |
|||
(loop |
|||
(setq X |
|||
(**Mod |
|||
(rand 2 (min (dec N) 1000000000000)) |
|||
Q |
|||
N ) ) |
|||
(T (or (=1 X) (= X N1))) |
|||
(T |
|||
(do K |
|||
(setq X (**Mod X 2 N)) |
|||
(when (=1 X) (throw 'composite)) |
|||
(T (= X N1) T) ) ) |
|||
(throw 'composite) ) ) |
|||
(throw 'composite T) ) ) ) ) ) ) |
|||
(for N 1000 |
|||
(and |
|||
(isprime (dec (** 2 N))) |
|||
(prinl "2 \^ " N " - 1") ) )</syntaxhighlight> |
|||
{{out}} |
|||
<pre> |
|||
2 ^ 2 - 1 |
|||
2 ^ 3 - 1 |
|||
2 ^ 5 - 1 |
|||
2 ^ 7 - 1 |
|||
2 ^ 13 - 1 |
|||
2 ^ 17 - 1 |
|||
2 ^ 19 - 1 |
|||
2 ^ 31 - 1 |
|||
2 ^ 61 - 1 |
|||
2 ^ 89 - 1 |
|||
2 ^ 107 - 1 |
|||
2 ^ 127 - 1 |
|||
2 ^ 521 - 1 |
|||
2 ^ 607 - 1 |
|||
</pre> |
|||
=={{header|Pike}}== |
=={{header|Pike}}== |
||
< |
<syntaxhighlight lang="pike">int power = 1; |
||
while(power++) { |
while(power++) { |
||
int candidate = 2->pow(power)-1; |
int candidate = 2->pow(power)-1; |
||
if( candidate->probably_prime_p() ) |
if( candidate->probably_prime_p() ) |
||
write("2 ^ %d - 1\n", power); |
write("2 ^ %d - 1\n", power); |
||
}</ |
}</syntaxhighlight> |
||
Output: |
Output: |
||
<pre> |
<pre> |
||
Line 943: | Line 1,681: | ||
=={{header|Prolog}}== |
=={{header|Prolog}}== |
||
Lucas-Lehmer test, works with SWI Prolog |
Lucas-Lehmer test, works with SWI Prolog |
||
< |
<syntaxhighlight lang="prolog"> |
||
lucas_lehmer_seq(M, L) :- |
lucas_lehmer_seq(M, L) :- |
||
lazy_list(ll_iter, 4-M, L). |
lazy_list(ll_iter, 4-M, L). |
||
Line 960: | Line 1,698: | ||
Skip is P - 3, drop(Skip, Residues, [R|_]), |
Skip is P - 3, drop(Skip, Residues, [R|_]), |
||
R =:= 0. |
R =:= 0. |
||
</syntaxhighlight> |
|||
</lang> |
|||
{{out}} |
{{out}} |
||
<pre> |
<pre> |
||
Line 970: | Line 1,708: | ||
=={{header|Python}}== |
=={{header|Python}}== |
||
{{trans|Java}} |
{{trans|Java}} |
||
< |
<syntaxhighlight lang="python">import random |
||
#Take from https://www.codeproject.com/Articles/691200/%2FArticles%2F691200%2FPrimality-test-algorithms-Prime-test-The-fastest-w |
#Take from https://www.codeproject.com/Articles/691200/%2FArticles%2F691200%2FPrimality-test-algorithms-Prime-test-The-fastest-w |
||
Line 1,053: | Line 1,791: | ||
if MillerRabinPrimalityTest(p): |
if MillerRabinPrimalityTest(p): |
||
break |
break |
||
print "done"</ |
print "done"</syntaxhighlight> |
||
{{out}} |
{{out}} |
||
<pre>2 ^ 2 - 1 |
<pre>2 ^ 2 - 1 |
||
Line 1,088: | Line 1,826: | ||
It doesn't specify to ''calculate'' them, only to ''list'' them; why throw away all of the computer '''millenia''' of processing power that the GIMPS has invested? |
It doesn't specify to ''calculate'' them, only to ''list'' them; why throw away all of the computer '''millenia''' of processing power that the GIMPS has invested? |
||
<lang |
<syntaxhighlight lang="raku" line>use HTTP::UserAgent; |
||
use Gumbo; |
use Gumbo; |
||
Line 1,098: | Line 1,836: | ||
for $table[1]».[*][0][*].comb(/'exp_lo='\d+/)».subst(/\D/, '',:g) |
for $table[1]».[*][0][*].comb(/'exp_lo='\d+/)».subst(/\D/, '',:g) |
||
.trans([<0123456789>.comb] => [<⁰¹²³⁴⁵⁶⁷⁸⁹>.comb]).words; |
.trans([<0123456789>.comb] => [<⁰¹²³⁴⁵⁶⁷⁸⁹>.comb]).words; |
||
</syntaxhighlight> |
|||
</lang> |
|||
{{out}} |
{{out}} |
||
<pre>All known Mersenne primes as of 2018-12-21 |
<pre>All known Mersenne primes as of 2018-12-21 |
||
Line 1,156: | Line 1,894: | ||
This REXX version (using a 32-bit Regina REXX interpreter) will find those Mersenne primes which are less than |
This REXX version (using a 32-bit Regina REXX interpreter) will find those Mersenne primes which are less than |
||
<br>8 million decimal digits (which would be '''M43'''). |
<br>8 million decimal digits (which would be '''M43'''). |
||
< |
<syntaxhighlight lang="rexx">/*REXX program uses exponent─and─mod operator to test possible Mersenne numbers. */ |
||
do j=1; /*process a range, or run out of time.*/ |
do j=1; /*process a range, or run out of time.*/ |
||
if \isPrime(j) then iterate /*if J isn't a prime, keep plugging.*/ |
if \isPrime(j) then iterate /*if J isn't a prime, keep plugging.*/ |
||
Line 1,196: | Line 1,934: | ||
end /*until*/ |
end /*until*/ |
||
if sq==1 then return q /*Not a prime? Return a factor.*/ |
if sq==1 then return q /*Not a prime? Return a factor.*/ |
||
end /*k*/</ |
end /*k*/</syntaxhighlight> |
||
<br><br> |
<br><br> |
||
=={{header|Ring}}== |
=={{header|Ring}}== |
||
< |
<syntaxhighlight lang="ring"> |
||
# Project : Mersenne primes |
# Project : Mersenne primes |
||
Line 1,218: | Line 1,956: | ||
next |
next |
||
return 1 |
return 1 |
||
</syntaxhighlight> |
|||
</lang> |
|||
Output: |
Output: |
||
<pre> |
<pre> |
||
Line 1,229: | Line 1,967: | ||
19 |
19 |
||
</pre> |
</pre> |
||
=={{header|RPL}}== |
|||
This version use binary integers |
|||
{{works with|Halcyon Calc|4.2.7}} |
|||
{| class="wikitable" |
|||
! RPL code |
|||
! Comment |
|||
|- |
|||
| |
|||
≪ / LAST ROT * - #0 == ≫ ''''BDIV?'''' STO |
|||
≪ |
|||
'''IF''' DUP #3 ≤ '''THEN''' #2 / B→R |
|||
'''ELSE''' |
|||
'''IF''' DUP #2 '''BDIV?''' OVER #3 '''BDIV?''' OR |
|||
'''THEN''' DROP 0 |
|||
'''ELSE''' |
|||
DUP B→R √ R→B → a maxd |
|||
≪ a #2 #5 |
|||
'''WHILE''' a OVER '''BDIV?''' NOT OVER maxd ≤ AND |
|||
'''REPEAT''' OVER + #6 ROT - SWAP '''END''' |
|||
≫ |
|||
SWAP DROP '''BDIV?''' NOT |
|||
'''END''' |
|||
'''END''' |
|||
≫ ''''BPRIM?'''' STO |
|||
≪ {} 1 32 '''FOR''' n |
|||
#1 1 n '''START''' SL '''NEXT''' #1 - |
|||
'''IF BPRIM? THEN''' + '''END''' |
|||
'''NEXT''' |
|||
≫ EVAL |
|||
| |
|||
''( #a #b -- boolean )'' |
|||
''( #a -- boolean )'' |
|||
return 1 if a is 2 or 3 |
|||
if 2 or 3 divides a |
|||
return 0 |
|||
else |
|||
store a and root(a) |
|||
initialize stack with a i d |
|||
while d does not divide a and d <= root(a) |
|||
i = 6 - i which modifies 2 into 4 and viceversa |
|||
convert stack status into result |
|||
Let's loop |
|||
Generate M(n) |
|||
test primality |
|||
|} |
|||
{{out}} |
|||
<pre> |
|||
{ 2 3 5 7 13 17 19 31 } |
|||
</pre> |
|||
=={{header|Ruby}}== |
|||
<syntaxhighlight lang="ruby">require 'openssl' |
|||
(0..).each{|n| puts "2**#{n} - 1" if OpenSSL::BN.new(2**n -1).prime? } |
|||
</syntaxhighlight> |
|||
{{out}} |
|||
<pre>2**2 - 1 |
|||
2**3 - 1 |
|||
2**5 - 1 |
|||
2**7 - 1 |
|||
2**13 - 1 |
|||
2**17 - 1 |
|||
2**19 - 1 |
|||
2**31 - 1 |
|||
2**61 - 1 |
|||
2**89 - 1 |
|||
2**107 - 1 |
|||
2**127 - 1 |
|||
2**521 - 1 |
|||
2**607 - 1 |
|||
2**1279 - 1 |
|||
2**2203 - 1 |
|||
2**2281 - 1 |
|||
2**3217 - 1 |
|||
2**4253 - 1 |
|||
^Ctest2.rb:7:in `prime?': Interrupt |
|||
</pre> |
|||
=={{header|Scala}}== |
=={{header|Scala}}== |
||
<syntaxhighlight lang="scala"> |
|||
<lang Scala> |
|||
object MersennePrimes extends App { |
object MersennePrimes extends App { |
||
Line 1,259: | Line 2,087: | ||
println("That's All Folks!") |
println("That's All Folks!") |
||
} |
} |
||
</syntaxhighlight> |
|||
</lang> |
|||
=={{header|Sidef}}== |
=={{header|Sidef}}== |
||
Uses the ''is_mersenne_prime()'' function from [https://metacpan.org/pod/Math::Prime::Util::GMP Math::Prime::Util::GMP]. |
Uses the ''is_mersenne_prime()'' function from [https://metacpan.org/pod/Math::Prime::Util::GMP Math::Prime::Util::GMP]. |
||
< |
<syntaxhighlight lang="ruby">for p in (^Inf -> lazy.grep { .is_mersenne_prime }) { |
||
say "2^#{p} - 1" |
say "2^#{p} - 1" |
||
}</ |
}</syntaxhighlight> |
||
{{out}} |
{{out}} |
||
<pre> |
<pre> |
||
Line 1,292: | Line 2,120: | ||
^C |
^C |
||
sidef mersenne.sf 12.47s user 0.02s system 99% cpu 12.495 total |
sidef mersenne.sf 12.47s user 0.02s system 99% cpu 12.495 total |
||
</pre> |
|||
=={{header|Transd}}== |
|||
<syntaxhighlight lang="Scheme">#lang transd |
|||
MainModule: { |
|||
_start: (λ locals: curexp 1 maxexp 1000 base BigLong(2) |
|||
(while (< curexp maxexp) |
|||
(+= curexp 1) |
|||
(if (not (is-prime BigLong(curexp))) continue) |
|||
(with cand (- (pow base curexp) 1) |
|||
(if (is-probable-prime cand 10) |
|||
(lout curexp " : " cand )) |
|||
) ) ) |
|||
}</syntaxhighlight> |
|||
{{out}} |
|||
<pre> |
|||
2 : 3 |
|||
3 : 7 |
|||
5 : 31 |
|||
7 : 127 |
|||
13 : 8191 |
|||
17 : 131071 |
|||
19 : 524287 |
|||
31 : 2147483647 |
|||
61 : 2305843009213693951 |
|||
89 : 618970019642690137449562111 |
|||
107 : 162259276829213363391578010288127 |
|||
127 : 170141183460469231731687303715884105727 |
|||
521 : 6864797660130609714981900799081393217269435300143305409394463459185543183397656052122559640661454554977296311391480858037121987999716643812574028291115057151 |
|||
607 : 531137992816767098689588206552468627329593117727031923199444138200403559860852242739162502265229285668889329486246501015346579337652707239409519978766587351943831270835393219031728127 |
|||
</pre> |
|||
=={{header|Visual Basic .NET}}== |
|||
{{trans|C#}} |
|||
<syntaxhighlight lang="vbnet">Imports System.Numerics |
|||
Module Module1 |
|||
Function Sqrt(x As BigInteger) As BigInteger |
|||
If x < 0 Then |
|||
Throw New ArgumentException("Negative argument.") |
|||
End If |
|||
If x < 2 Then |
|||
Return x |
|||
End If |
|||
Dim y = x / 2 |
|||
While y > (x / y) |
|||
y = ((x / y) + y) / 2 |
|||
End While |
|||
Return y |
|||
End Function |
|||
Function IsPrime(bi As BigInteger) As Boolean |
|||
If bi < 2 Then |
|||
Return False |
|||
End If |
|||
If bi Mod 2 = 0 Then |
|||
Return bi = 2 |
|||
End If |
|||
If bi Mod 3 = 0 Then |
|||
Return bi = 3 |
|||
End If |
|||
If bi Mod 5 = 0 Then |
|||
Return bi = 5 |
|||
End If |
|||
If bi Mod 7 = 0 Then |
|||
Return bi = 7 |
|||
End If |
|||
If bi Mod 11 = 0 Then |
|||
Return bi = 11 |
|||
End If |
|||
If bi Mod 13 = 0 Then |
|||
Return bi = 13 |
|||
End If |
|||
If bi Mod 17 = 0 Then |
|||
Return bi = 17 |
|||
End If |
|||
If bi Mod 19 = 0 Then |
|||
Return bi = 19 |
|||
End If |
|||
Dim limit = Sqrt(bi) |
|||
Dim test As BigInteger = 23 |
|||
While test < limit |
|||
If bi Mod test = 0 Then |
|||
Return False |
|||
End If |
|||
test += 2 |
|||
If bi Mod test = 0 Then |
|||
Return False |
|||
End If |
|||
test += 4 |
|||
End While |
|||
Return True |
|||
End Function |
|||
Sub Main() |
|||
Const MAX = 9 |
|||
Dim pow = 2 |
|||
Dim count = 0 |
|||
While True |
|||
If IsPrime(pow) Then |
|||
Dim p = BigInteger.Pow(2, pow) - 1 |
|||
If IsPrime(p) Then |
|||
Console.WriteLine("2 ^ {0} - 1", pow) |
|||
count += 1 |
|||
If count >= MAX Then |
|||
Exit While |
|||
End If |
|||
End If |
|||
End If |
|||
pow += 1 |
|||
End While |
|||
End Sub |
|||
End Module</syntaxhighlight> |
|||
{{out}} |
|||
<pre>2 ^ 2 - 1 |
|||
2 ^ 3 - 1 |
|||
2 ^ 5 - 1 |
|||
2 ^ 7 - 1 |
|||
2 ^ 13 - 1 |
|||
2 ^ 17 - 1 |
|||
2 ^ 19 - 1 |
|||
2 ^ 31 - 1 |
|||
2 ^ 61 - 1</pre> |
|||
=={{header|Wren}}== |
|||
===Wren-CLI=== |
|||
{{libheader|Wren-math}} |
|||
{{libheader|Wren-big}} |
|||
A bit slow so limited to first 14 Mersenne primes. |
|||
<syntaxhighlight lang="wren">import "./math" for Int |
|||
import "./big" for BigInt |
|||
var MAX = 14 |
|||
System.print("The first %(MAX) Mersenne primes are:") |
|||
var count = 0 |
|||
var p = 2 |
|||
while (true) { |
|||
var m = (BigInt.one << p) - 1 |
|||
if (m.isProbablePrime(10)) { |
|||
System.print("2 ^ %(p) - 1") |
|||
count = count + 1 |
|||
if (count == MAX) break |
|||
} |
|||
while (true) { |
|||
p = (p > 2) ? p + 2 : 3 |
|||
if (Int.isPrime(p)) break |
|||
} |
|||
}</syntaxhighlight> |
|||
{{out}} |
|||
<pre> |
|||
The first 14 Mersenne primes are: |
|||
2 ^ 2 - 1 |
|||
2 ^ 3 - 1 |
|||
2 ^ 5 - 1 |
|||
2 ^ 7 - 1 |
|||
2 ^ 13 - 1 |
|||
2 ^ 17 - 1 |
|||
2 ^ 19 - 1 |
|||
2 ^ 31 - 1 |
|||
2 ^ 61 - 1 |
|||
2 ^ 89 - 1 |
|||
2 ^ 107 - 1 |
|||
2 ^ 127 - 1 |
|||
2 ^ 521 - 1 |
|||
2 ^ 607 - 1 |
|||
</pre> |
|||
===Embedded (GMP)=== |
|||
{{libheader|Wren-gmp}} |
|||
This finds the first 23 Mersenne primes in about 172 seconds which is virtually the same as the Go non-concurrent version using their GMP plug-in when run on my machine. |
|||
<syntaxhighlight lang="wren">import "./math" for Int |
|||
import "./gmp" for Mpz |
|||
var MAX = 23 |
|||
System.print("The first %(MAX) Mersenne primes are:") |
|||
var count = 0 |
|||
var p = 2 |
|||
while (true) { |
|||
var m = Mpz.one.lsh(p).sub(1) |
|||
if (m.probPrime(15) > 0) { |
|||
System.print("2 ^ %(p) - 1") |
|||
count = count + 1 |
|||
if (count == MAX) break |
|||
} |
|||
while (true) { |
|||
p = (p > 2) ? p + 2 : 3 |
|||
if (Int.isPrime(p)) break |
|||
} |
|||
}</syntaxhighlight> |
|||
{{out}} |
|||
<pre> |
|||
The first 23 Mersenne primes are: |
|||
2 ^ 2 - 1 |
|||
2 ^ 3 - 1 |
|||
2 ^ 5 - 1 |
|||
2 ^ 7 - 1 |
|||
2 ^ 13 - 1 |
|||
2 ^ 17 - 1 |
|||
2 ^ 19 - 1 |
|||
2 ^ 31 - 1 |
|||
2 ^ 61 - 1 |
|||
2 ^ 89 - 1 |
|||
2 ^ 107 - 1 |
|||
2 ^ 127 - 1 |
|||
2 ^ 521 - 1 |
|||
2 ^ 607 - 1 |
|||
2 ^ 1279 - 1 |
|||
2 ^ 2203 - 1 |
|||
2 ^ 2281 - 1 |
|||
2 ^ 3217 - 1 |
|||
2 ^ 4253 - 1 |
|||
2 ^ 4423 - 1 |
|||
2 ^ 9689 - 1 |
|||
2 ^ 9941 - 1 |
|||
2 ^ 11213 - 1 |
|||
</pre> |
|||
=={{header|XPL0}}== |
|||
<syntaxhighlight lang="xpl0">func IsPrime(N); \Return 'true' if N is prime |
|||
int N, I; |
|||
[if N <= 2 then return N = 2; |
|||
if (N&1) = 0 then \even >2\ return false; |
|||
for I:= 3 to sqrt(N) do |
|||
[if rem(N/I) = 0 then return false; |
|||
I:= I+1; |
|||
]; |
|||
return true; |
|||
]; |
|||
int N; |
|||
[for N:= 1 to 31 do |
|||
if IsPrime(1<<N-1) then |
|||
[Text(0, "2^^"); IntOut(0, N); Text(0, " - 1"); |
|||
CrLf(0); |
|||
]; |
|||
]</syntaxhighlight> |
|||
{{out}} |
|||
<pre> |
|||
2^2 - 1 |
|||
2^3 - 1 |
|||
2^5 - 1 |
|||
2^7 - 1 |
|||
2^13 - 1 |
|||
2^17 - 1 |
|||
2^19 - 1 |
|||
2^31 - 1 |
|||
</pre> |
</pre> |
||
Line 1,297: | Line 2,381: | ||
{{libheader|GMP}} |
{{libheader|GMP}} |
||
Uses libGMP (GNU MP Bignum Library) and its Miller-Rabin probabilistic primality testing algorithm. |
Uses libGMP (GNU MP Bignum Library) and its Miller-Rabin probabilistic primality testing algorithm. |
||
< |
<syntaxhighlight lang="zkl">var [const] BN=Import.lib("zklBigNum"); // libGMP |
||
fcn mprimes{ |
fcn mprimes{ |
||
n,m := BN(2),0; |
n,m := BN(2),0; |
||
Line 1,305: | Line 2,389: | ||
} |
} |
||
}() |
}() |
||
// gets rather slow after M(4423)</ |
// gets rather slow after M(4423)</syntaxhighlight> |
||
{{out}} |
{{out}} |
||
<pre style="height:40ex"> |
<pre style="height:40ex"> |
Latest revision as of 21:48, 17 March 2024
- Task
Create code that will list (preferably calculate) all of the Mersenne primes until some limitation is reached.
The number of known Mersenne primes is 51 (as of June, 2020), and the largest known Mersenne prime contains contains 24,862,048 decimal digits.
- Also see
- the Wikipedia entry: Mersenne prime.
- the MathWorld entry; Mersenne prime.
- For a list of all the know Mersenne primes: [https://primes.utm.edu/mersenne/index.html#known list of Mersenne
- For a general website about primes: prime pages.
- the OEIS entry: Mersenne primes.
- the OEIS entry: A000043 Mersenne exponents: primes p such that 2^p - 1 is prime. Then 2^p - 1 is called a Mersenne prime.
11l
F is_prime(BigInt bi)
I bi < 2 {R 0B}
I bi % 2 == 0 {R bi == 2}
I bi % 3 == 0 {R bi == 3}
V test = BigInt(5)
L test * test < bi
I bi % test == 0
R 0B
test += 2
I bi % test == 0
R 0B
test += 4
R 1B
V base = BigInt(2)
L(p) 1..31
I is_prime(base - 1)
print(‘2 ^ ’p‘ - 1’)
base *= 2
- Output:
2 ^ 2 - 1 2 ^ 3 - 1 2 ^ 5 - 1 2 ^ 7 - 1 2 ^ 13 - 1 2 ^ 17 - 1 2 ^ 19 - 1 2 ^ 31 - 1
ALGOL 60
begin
integer procedure mersenne(n);
value n; integer n;
begin
integer i, m;
m := 1;
for i := 1 step 1 until n do
m := m * 2;
mersenne := m - 1;
end;
boolean procedure isprime(n);
value n; integer n;
begin
if n < 2 then
isprime := false
else if entier(n / 2) * 2 = n then
isprime := (n = 2)
else
begin
comment - check odd divisors up to sqrt(n);
integer i, limit;
boolean divisible;
i := 3;
limit := entier(sqrt(n));
divisible := false;
for i := i while i <= limit and not divisible do
begin
if entier(n / i) * i = n then
divisible := true;
i := i + 2
end;
isprime := if divisible then false else true;
end;
end;
comment - main code begins here;
integer i, m;
outstring(1,"Searching to M(31) for Mersenne primes\n");
for i := 1 step 1 until 31 do
begin
m := mersenne(i);
if isprime(m) then
begin
outstring(1,"M(");
outinteger(1,i);
outstring(1,") : ");
outinteger(1,m);
outstring(1,"\n");
end;
end;
end
- Output:
Searching to M(31) for Mersenne primes M( 2 ) : 3 M( 3 ) : 7 M( 5 ) : 31 M( 7 ) : 127 M( 13 ) : 8191 M( 17 ) : 131071 M( 19 ) : 524287 M( 31 ) : 2147483647
ALGOL 68
BEGIN # find some Mersenne Primrs - primes of the form 2^n - 1, n must be prime #
# This assumes LONG INT is at least 64 bits (as in e.g. Algol 68G) #
# we handle 2 as a special case and then the odd numbers starting at 3 #
PR read "primes.incl.a68" PR # include prime utilities #
# 2^0 - 1 = 0 and 2^1 - 1 = 1, neither of which is prime #
# start from 2^2 #
INT n := 2;
LONG INT p2 := 4;
IF is probably prime( p2 - 1 ) THEN print( ( " ", whole( n, 0 ) ) ) FI;
# 2^3, 2^5, etc. #
n +:= 1;
p2 *:= 2;
WHILE n < 61 DO
IF is probably prime( p2 - 1 ) THEN print( ( " ", whole( n, 0 ) ) ) FI;
n +:= 2;
p2 *:= 4
OD;
IF is probably prime( p2 - 1 ) THEN print( ( " ", whole( n, 0 ) ) ) FI
END
- Output:
2 3 5 7 13 17 19 31 61
ALGOL W
begin % find some Mersenne Primrs - primes of the form 2^n - 1, n must be prime %
% integers are 32 bit in Algol W, so there won't be many numbers to check %
% we handle 2 as a special case and then the odd numbers starting at 3 %
% simplified primality by trial division test - no need to check for even %
% numbers as 2^n - 1 is always odd for n >= 2 %
logical procedure oddOnlyPrimalityTest ( integer value n ) ;
begin
logical isPrime;
isPrime := true;
for i := 3 step 2 until entier( sqrt( n ) ) do begin;
isPrime := n rem i not = 0;
if not isPrime then goto endPrimalityTest
end for_i;
endPrimalityTest: isPrime
end oddOnlyPrimalityTest ;
integer p2, n;
n := 1;
p2 := 2;
while
begin
if n < 3 then begin
n := n + 1;
p2 := p2 * 2
end
else begin
n := n + 2;
p2 := p2 * 4
end if_n_le_3__ ;;
if oddOnlyPrimalityTest( p2 - 1 ) then writeon( i_w := 1, s_w := 0, " ", n );
n < 29
end
do begin end ;
% MAXINTEGER is 2**31 - 1 %
if oddOnlyPrimalityTest( MAXINTEGER ) then writeon( i_w := 1, s_w := 0, " ", 31 );
end.
- Output:
2 3 5 7 13 17 19 31
AppleScript
on isPrime(integ)
set isComposite to ""
if (integ / 2) = (integ / 2 div 1) then
log integ & " is composite because 2 is a factor" as string --buttons {"OK", "Cancel"} default button 1 cancel button 2
else
set x to 2
set sqrtOfInteg to integ ^ 0.5
repeat until x = integ ^ 0.5 + 1 as integer
if (integ / x) = integ / x div 1 then
log integ & " is composite because " & x & " & " & (integ / x div 1) & " are factors" as string --buttons {"OK", "Cancel"} default button 1 cancel button 2
set isComposite to 1
set x to x + 1
else
set x to x + 1
end if
end repeat
log integ & " is prime" as string --buttons {"OK", "Cancel"} default button 1 cancel button 2
if isComposite = 1 then
log integ & "is composite"
else
display dialog integ
end if
end if
end isPrime
set x to 2
repeat
isPrime(((2 ^ x) - 1) div 1)
set x to x + 1
end repeat
Or more efficiently:
on isPrime(n)
if (n < 4) then return (n > 1)
if ((n mod 2 is 0) or (n mod 3 is 0)) then return false
set {i, max} to {5, (n ^ 0.5) div 1}
repeat until (i > max)
if ((n mod i is 0) or (n mod (i + 2) is 0)) then return false
set i to i + 6
end repeat
return true
end isPrime
on join(lst, delim)
set astid to AppleScript's text item delimiters
set AppleScript's text item delimiters to delim
set txt to lst as text
set AppleScript's text item delimiters to astid
return txt
end join
on mersennePrimes()
set output to {"Mersenne primes within AppleScript's number precision:"}
-- Special-case 2 ^ 2 - 1.
set end of output to "2 ^ 2 - 1 = 3"
set p to 1 -- Otherwise test odd-numbered powers of 2.
try -- Survive the "numeric operation too large" error when it occurs.
repeat
set p to p + 2
if ((isPrime(p)) and (isPrime(2 ^ p - 1))) then ¬
set end of output to "2 ^ " & p & " - 1 = " & (2 ^ p div 1 - 1)
end repeat
end try
return join(output, linefeed)
end mersennePrimes
mersennePrimes()
- Output:
"Mersenne primes within AppleScript's number precision:
2 ^ 2 - 1 = 3
2 ^ 3 - 1 = 7
2 ^ 5 - 1 = 31
2 ^ 7 - 1 = 127
2 ^ 13 - 1 = 8191
2 ^ 17 - 1 = 131071
2 ^ 19 - 1 = 524287
2 ^ 31 - 1 = 2.147483647E+9"
Arturo
mersenne?: function [n][
prime? dec 2 ^ n
]
1..31 | select => mersenne?
| loop 'x -> print ["M (" x ") = 2 ^" x " - 1 =" (2^x)-1 "is a prime"]
- Output:
M ( 2 ) = 2 ^ 2 - 1 = 3 is a prime M ( 3 ) = 2 ^ 3 - 1 = 7 is a prime M ( 5 ) = 2 ^ 5 - 1 = 31 is a prime M ( 7 ) = 2 ^ 7 - 1 = 127 is a prime M ( 13 ) = 2 ^ 13 - 1 = 8191 is a prime M ( 17 ) = 2 ^ 17 - 1 = 131071 is a prime M ( 19 ) = 2 ^ 19 - 1 = 524287 is a prime M ( 31 ) = 2 ^ 31 - 1 = 2147483647 is a prime
AWK
# syntax: GAWK --bignum -f MERSENNE_PRIMES.AWK
BEGIN {
base = 2
for (i=1; i<62; i++) {
if (is_prime(base-1)) {
printf("2 ^ %d - 1\n",i)
}
base *= 2
}
exit(0)
}
function is_prime(n, d) {
d = 5
if (n < 2) { return(0) }
if (n % 2 == 0) { return(n == 2) }
if (n % 3 == 0) { return(n == 3) }
while (d*d <= n) {
if (n % d == 0) { return(0) }
d += 2
if (n % d == 0) { return(0) }
d += 4
}
return(1)
}
- Output:
2 ^ 2 - 1 2 ^ 3 - 1 2 ^ 5 - 1 2 ^ 7 - 1 2 ^ 13 - 1 2 ^ 17 - 1 2 ^ 19 - 1 2 ^ 31 - 1 2 ^ 61 - 1
BASIC
BASIC256
m = 0
while True
m += 1
if isPrime((2^m)-1) = True then print m
end while
end
function isPrime(v)
if v <= 1 then return False
for i = 2 To int(sqr(v))
if v % i = 0 then return False
next i
return True
end function
- Output:
Igual que la entrada de FreeBASIC.
FreeBASIC
Function isPrime(Byval ValorEval As Integer) As Boolean
If ValorEval <= 1 Then Return False
For i As Integer = 2 To Int(Sqr(ValorEval))
If ValorEval Mod i = 0 Then Return False
Next i
Return True
End Function
Dim As Integer m = 0
While True
m += 1
If isPrime((2^m)-1) Then Print m
Wend
Sleep
- Output:
2 3 5 7 13 17 19 31
Yabasic
m = 0
while True
m = m + 1
if isPrime((2^m)-1) = True then print m : fi
wend
end
sub isPrime(v)
if v < 2 then return False : fi
if mod(v, 2) = 0 then return v = 2 : fi
if mod(v, 3) = 0 then return v = 3 : fi
d = 5
while d * d <= v
if mod(v, d) = 0 then return False else d = d + 2 : fi
end while
return True
end sub
- Output:
Igual que la entrada de FreeBASIC.
C
#include <stdbool.h>
#include <stdint.h>
#include <stdio.h>
bool isPrime(uint64_t n) {
uint64_t test;
if (n < 2) return false;
if (n % 2 == 0) return n == 2;
if (n % 3 == 0) return n == 3;
test = 5;
while (test * test < n) {
if (n % test == 0) return false;
test += 2;
if (n % test == 0) return false;
test += 4;
}
return true;
}
int main() {
uint64_t base = 2;
int pow;
for (pow = 1; pow <= 32; pow++) {
if (isPrime(base - 1)) {
printf("2 ^ %d - 1\n", pow);
}
base *= 2;
}
return 0;
}
- Output:
2 ^ 2 - 1 2 ^ 3 - 1 2 ^ 5 - 1 2 ^ 7 - 1 2 ^ 13 - 1 2 ^ 17 - 1 2 ^ 19 - 1 2 ^ 31 - 1
Alternatively, we can use GMP to find the first 23 Mersenne primes in about the same time as the corresponding Wren entry.
#include <stdio.h>
#include <stdbool.h>
#include <stdint.h>
#include <gmp.h>
#define MAX 23
bool isPrime(uint64_t n) {
uint64_t test;
if (n < 2) return false;
if (n % 2 == 0) return n == 2;
if (n % 3 == 0) return n == 3;
test = 5;
while (test * test < n) {
if (n % test == 0) return false;
test += 2;
if (n % test == 0) return false;
test += 4;
}
return true;
}
int main() {
uint64_t p = 2;
int count = 0;
mpz_t m, one;
mpz_init(m);
mpz_init_set_ui(one, 1);
while (true) {
mpz_mul_2exp(m, one, p);
mpz_sub_ui(m, m, 1);
if (mpz_probab_prime_p(m, 15) > 0) {
printf("2 ^ %ld - 1\n", p);
if (++count == MAX) break;
}
while (true) {
p = (p > 2) ? p + 2 : 3;
if (isPrime(p)) break;
}
}
mpz_clear(m);
mpz_clear(one);
return 0;
}
- Output:
Same as Wren example.
C++
#include <iostream>
bool isPrime(uint64_t n) {
if (n < 2) return false;
if (n % 2 == 0) return n == 2;
if (n % 3 == 0) return n == 3;
uint64_t test = 5;
while (test * test < n) {
if (n % test == 0) return false;
test += 2;
if (n % test == 0) return false;
test += 4;
}
return true;
}
int main() {
uint64_t base = 2;
for (int pow = 1; pow <= 32; pow++) {
if (isPrime(base - 1)) {
std::cout << "2 ^ " << pow << " - 1\n";
}
base *= 2;
}
return 0;
}
- Output:
2 ^ 2 - 1 2 ^ 3 - 1 2 ^ 5 - 1 2 ^ 7 - 1 2 ^ 13 - 1 2 ^ 17 - 1 2 ^ 19 - 1 2 ^ 31 - 1
C#
Needs a better primality checking algorithm to do really large prime numbers.
using System;
using System.Numerics;
namespace MersennePrimes {
class Program {
static BigInteger Sqrt(BigInteger x) {
if (x < 0) throw new ArgumentException("Negative argument.");
if (x < 2) return x;
BigInteger y = x / 2;
while (y > x / y) {
y = ((x / y) + y) / 2;
}
return y;
}
static bool IsPrime(BigInteger bi) {
if (bi < 2) return false;
if (bi % 2 == 0) return bi == 2;
if (bi % 3 == 0) return bi == 3;
if (bi % 5 == 0) return bi == 5;
if (bi % 7 == 0) return bi == 7;
if (bi % 11 == 0) return bi == 11;
if (bi % 13 == 0) return bi == 13;
if (bi % 17 == 0) return bi == 17;
if (bi % 19 == 0) return bi == 19;
BigInteger limit = Sqrt(bi);
BigInteger test = 23;
while (test < limit) {
if (bi % test == 0) return false;
test += 2;
if (bi % test == 0) return false;
test += 4;
}
return true;
}
static void Main(string[] args) {
const int MAX = 9;
int pow = 2;
int count = 0;
while (true) {
if (IsPrime(pow)) {
BigInteger p = BigInteger.Pow(2, pow) - 1;
if (IsPrime(p)) {
Console.WriteLine("2 ^ {0} - 1", pow);
if (++count >= MAX) {
break;
}
}
}
pow++;
}
}
}
}
- Output:
2 ^ 2 - 1 2 ^ 3 - 1 2 ^ 5 - 1 2 ^ 7 - 1 2 ^ 13 - 1 2 ^ 17 - 1 2 ^ 19 - 1 2 ^ 31 - 1 2 ^ 61 - 1
D
Simplest thing that could possibly work. Using better primality tests will allow for more results to be calculated in a reasonable amount of time.
import std.bigint;
import std.stdio;
bool isPrime(BigInt bi) {
if (bi < 2) return false;
if (bi % 2 == 0) return bi == 2;
if (bi % 3 == 0) return bi == 3;
auto test = BigInt(5);
while (test * test < bi) {
if (bi % test == 0) return false;
test += 2;
if (bi % test == 0) return false;
test += 4;
}
return true;
}
void main() {
auto base = BigInt(2);
for (int pow=1; pow<32; pow++) {
if (isPrime(base-1)) {
writeln("2 ^ ", pow, " - 1");
}
base *= 2;
}
}
- Output:
2 ^ 2 - 1 2 ^ 3 - 1 2 ^ 5 - 1 2 ^ 7 - 1 2 ^ 13 - 1 2 ^ 17 - 1 2 ^ 19 - 1 2 ^ 31 - 1
Delphi
function IsPrime(N: int64): boolean;
{Fast, optimised prime test}
var I,Stop: int64;
begin
if (N = 2) or (N=3) then Result:=true
else if (n <= 1) or ((n mod 2) = 0) or ((n mod 3) = 0) then Result:= false
else
begin
I:=5;
Stop:=Trunc(sqrt(N+0.0));
Result:=False;
while I<=Stop do
begin
if ((N mod I) = 0) or ((N mod (I + 2)) = 0) then exit;
Inc(I,6);
end;
Result:=True;
end;
end;
procedure MersennePrimes(Memo: TMemo);
var N: integer;
var Mn: int64;
begin
Memo.Lines.Add('2^N-1 Prime');
Memo.Lines.Add('--------------------' );
for N:=1 to 32 do
begin
Mn:=(1 shl N)-1;
if IsPrime(Mn) then
Memo.Lines.Add(Format('2^%d-1 %14.0n',[N,Mn+0.0]));
end;
end;
- Output:
2^N-1 Prime -------------------- 2^2-1 3 2^3-1 7 2^5-1 31 2^7-1 127 2^13-1 8,191 2^17-1 131,071 2^19-1 524,287 2^31-1 2,147,483,647
EasyLang
fastfunc isprim num .
if num mod 2 = 0
if num = 2
return 1
.
return 0
.
i = 3
while i <= sqrt num
if num mod i = 0
return 0
.
i += 2
.
return 1
.
base = 4
for p = 2 to 52
if isprim (base - 1) = 1
print "2 ^ " & p & " - 1"
.
base *= 2
.
- Output:
2 ^ 2 - 1 2 ^ 3 - 1 2 ^ 5 - 1 2 ^ 7 - 1 2 ^ 13 - 1 2 ^ 17 - 1 2 ^ 19 - 1 2 ^ 31 - 1
F#
open System
open System.Numerics
let Sqrt (n:BigInteger) =
if n < (BigInteger 0) then raise (ArgumentException "Negative argument.")
if n < (BigInteger 2) then n
else
let rec H v r s =
if v < s then
r
else
H (v - s) (r + (BigInteger 1)) (s + (BigInteger 2))
H n (BigInteger 0) (BigInteger 1)
let IsPrime (n:BigInteger) =
if n < (BigInteger 2) then false
elif n % (BigInteger 2) = (BigInteger 0) then n = (BigInteger 2)
elif n % (BigInteger 3) = (BigInteger 0) then n = (BigInteger 3)
elif n % (BigInteger 5) = (BigInteger 0) then n = (BigInteger 5)
elif n % (BigInteger 7) = (BigInteger 0) then n = (BigInteger 7)
elif n % (BigInteger 11) = (BigInteger 0) then n = (BigInteger 11)
elif n % (BigInteger 13) = (BigInteger 0) then n = (BigInteger 13)
elif n % (BigInteger 17) = (BigInteger 0) then n = (BigInteger 17)
elif n % (BigInteger 19) = (BigInteger 0) then n = (BigInteger 19)
else
let limit = (Sqrt n)
let rec H t =
if t <= limit then
if n % t = (BigInteger 0) then false
else
let t2 = t + (BigInteger 2)
if n % t2 = (BigInteger 0) then false
else H (t2 + (BigInteger 4))
else
true
H (BigInteger 23)
[<EntryPoint>]
let main _ =
let MAX = BigInteger 9
let rec loop (pow:int) (count:int) =
if IsPrime (BigInteger pow) then
let p = BigInteger.Pow((BigInteger 2), pow) - (BigInteger 1)
if IsPrime p then
printfn "2 ^ %A - 1" pow
if (BigInteger (count + 1)) >= MAX then count
else loop (pow + 1) (count + 1)
else loop (pow + 1) count
else loop (pow + 1) count
loop 2 0 |> ignore
0 // return an integer exit code
- Output:
2 ^ 2 - 1 2 ^ 3 - 1 2 ^ 5 - 1 2 ^ 7 - 1 2 ^ 13 - 1 2 ^ 17 - 1 2 ^ 19 - 1 2 ^ 31 - 1 2 ^ 61 - 1
Factor
Factor comes with a Lucas-Lehmer primality test.
USING: formatting math.primes.lucas-lehmer math.ranges sequences ;
: mersennes-upto ( n -- seq ) [1,b] [ lucas-lehmer ] filter ;
3500 mersennes-upto [ "2 ^ %d - 1\n" printf ] each
- Output:
2 ^ 2 - 1 2 ^ 3 - 1 2 ^ 5 - 1 2 ^ 7 - 1 2 ^ 13 - 1 2 ^ 17 - 1 2 ^ 19 - 1 2 ^ 31 - 1 2 ^ 61 - 1 2 ^ 89 - 1 2 ^ 107 - 1 2 ^ 127 - 1 2 ^ 521 - 1 2 ^ 607 - 1 2 ^ 1279 - 1 2 ^ 2203 - 1 2 ^ 2281 - 1 2 ^ 3217 - 1
Fortran
program mersenne
use iso_fortran_env, only: output_unit, INT64
implicit none
integer, parameter :: l=INT64
integer(kind=l) :: base
integer :: pow
base = 2
do pow = 1, 32
if (is_prime(base-1)) then
write(output_unit,'(A2,x,I0,x,A3)') "2^", pow, "- 1"
end if
base = base * 2
end do
contains
pure function is_prime(n)
integer(kind=l), intent(in) :: n
logical :: is_prime
integer(kind=l) :: test
is_prime = .false.
if (n < 2) return
if (modulo(n, 2) == 0) then
is_prime = n==2
return
end if
if (modulo(n, 3) == 0) then
is_prime = n==3
return
end if
test = 5
do
if (test**2 >= n) then
is_prime = .true.
return
end if
if (modulo(n, test) == 0) return
test = test + 2
if (modulo(n, test) == 0) return
test = test + 4
end do
end function is_prime
end program mersenne
- Output:
2^ 2 - 1 2^ 3 - 1 2^ 5 - 1 2^ 7 - 1 2^ 13 - 1 2^ 17 - 1 2^ 19 - 1 2^ 31 - 1
Frink
Frink has built-in routines for iterating through prime numbers. Frink's isPrime[n]
function automatically detects numbers of the form 2n-1 and performs a more efficient Lucas-Lehmer primality test on the number. This works with arbitrarily large numbers.
Did you know that all Java-based JVM languages are many many orders of magnitude faster because Frink's developer contributed vastly faster BigInteger algorithms to Java? It took the Java developers 11 years to integrate them but they became part of 1.8 and later! Operations that used to take days now can be done in seconds thanks to Frink's contributions to Java.
for n = primes[]
if isPrime[2^n - 1]
println["2^$n - 1"]
- Output:
2^2 - 1 2^3 - 1 2^5 - 1 2^7 - 1 2^13 - 1 2^17 - 1 2^19 - 1 2^31 - 1 2^61 - 1 2^89 - 1 2^107 - 1 2^127 - 1 2^521 - 1 2^607 - 1 2^1279 - 1 2^2203 - 1 2^2281 - 1 2^3217 - 1 2^4253 - 1 2^4423 - 1 2^9689 - 1 2^9941 - 1 2^11213 - 1 2^19937 - 1 2^21701 - 1 ...
Go
The github.com/ncw/gmp
package is a drop-in replacement for Go's math/big
package.
It's a CGo wrapper around the C GMP library and under these circumstances is two to four times as fast as the native Go package.
Editing just the import line you can use whichever is more convenient for you
(CGo has drawbacks, including limited portability).
Normally build tags would be used to control this instead of editing imports in the source, but this keeps the example simpler.
Note that the use of ProbablyPrime(0) requires Go 1.8 or later. When using the math/big
package, passing a parameter of zero to this method forces it to apply only the Baillie-PSW test to check for primality. This is 100% accurate for numbers up to 2^64 and at the time of writing (June 2018) no known composite number above that bound passes the test.
package main
import (
"fmt"
"time"
// Use one or the other of these:
"math/big"
//big "github.com/ncw/gmp"
)
func main() {
start := time.Now()
one := big.NewInt(1)
mp := big.NewInt(0)
bp := big.NewInt(0)
const max = 22
for count, p := 0, uint(2); count < max; {
mp.Lsh(one, p)
mp.Sub(mp, one)
if mp.ProbablyPrime(0) {
elapsed := time.Since(start).Seconds()
if elapsed >= 0.01 {
fmt.Printf("2 ^ %-4d - 1 took %6.2f secs\n", p, elapsed)
} else {
fmt.Printf("2 ^ %-4d - 1\n", p)
}
count++
}
for {
if p > 2 {
p += 2
} else {
p = 3
}
bp.SetUint64(uint64(p))
if bp.ProbablyPrime(0) {
break
}
}
}
}
- Output using the GMP package on a 3.4 GHz Xeon E3-1245:
2 ^ 2 - 1 2 ^ 3 - 1 2 ^ 5 - 1 2 ^ 7 - 1 2 ^ 13 - 1 2 ^ 17 - 1 2 ^ 19 - 1 2 ^ 31 - 1 2 ^ 61 - 1 2 ^ 89 - 1 2 ^ 107 - 1 2 ^ 127 - 1 2 ^ 521 - 1 2 ^ 607 - 1 2 ^ 1279 - 1 took 0.05 secs 2 ^ 2203 - 1 took 0.38 secs 2 ^ 2281 - 1 took 0.44 secs 2 ^ 3217 - 1 took 1.53 secs 2 ^ 4253 - 1 took 4.39 secs 2 ^ 4423 - 1 took 5.02 secs 2 ^ 9689 - 1 took 73.78 secs 2 ^ 9941 - 1 took 81.24 secs
(A previous run on more modest hardware - Celeron N3050 @ 1.60GHz × 2 - was ~365 seconds for M9941.)
This can be sped up quite a bit for modern multi-core CPUs by some simple changes to use goroutines.
package main
import (
"fmt"
"runtime"
"time"
// Use one or the other of these:
"math/big"
//big "github.com/ncw/gmp"
)
func main() {
start := time.Now()
nworkers := runtime.GOMAXPROCS(0)
fmt.Println("Using", nworkers, "workers.")
workC := make(chan uint, 1)
resultC := make(chan uint, nworkers)
// Generate possible Mersenne exponents and send them to workC.
go func() {
workC <- 2
bp := big.NewInt(0)
for p := uint(3); ; p += 2 {
// Possible exponents must be prime.
bp.SetUint64(uint64(p))
if bp.ProbablyPrime(0) {
workC <- p
}
}
}()
// Start up worker go routines, each takes
// possible Mersenne exponents from workC as `p`
// and if 2^p-1 is prime sends `p` to resultC.
one := big.NewInt(1)
for i := 0; i < nworkers; i++ {
go func() {
mp := big.NewInt(0)
for p := range workC {
mp.Lsh(one, p)
mp.Sub(mp, one)
if mp.ProbablyPrime(0) {
resultC <- p
}
}
}()
}
// Receive some maximum number of Mersenne prime exponents
// from resultC and show the Mersenne primes.
const max = 24
for count := 0; count < max; count++ {
// Note: these could come back out of order, although usually
// only the first few. If that is an issue, correcting it is
// left as an excercise to the reader :).
p := <-resultC
elapsed := time.Since(start).Seconds()
if elapsed >= 0.01 {
fmt.Printf("2 ^ %-5d - 1 took %6.2f secs\n", p, elapsed)
} else {
fmt.Printf("2 ^ %-5d - 1\n", p)
}
}
}
- Output using the GMP package on the same 3.4 GHz Xeon E3-1245 (4 core × 2 SMT threads) as above:
Using 8 workers. 2 ^ 2 - 1 2 ^ 5 - 1 2 ^ 3 - 1 2 ^ 7 - 1 2 ^ 13 - 1 2 ^ 19 - 1 2 ^ 61 - 1 2 ^ 31 - 1 2 ^ 107 - 1 2 ^ 17 - 1 2 ^ 127 - 1 2 ^ 89 - 1 2 ^ 521 - 1 2 ^ 607 - 1 2 ^ 1279 - 1 took 0.01 secs 2 ^ 2203 - 1 took 0.09 secs 2 ^ 2281 - 1 took 0.12 secs 2 ^ 3217 - 1 took 0.36 secs 2 ^ 4253 - 1 took 0.94 secs 2 ^ 4423 - 1 took 1.06 secs 2 ^ 9689 - 1 took 16.28 secs 2 ^ 9941 - 1 took 18.02 secs 2 ^ 11213 - 1 took 26.76 secs 2 ^ 19937 - 1 took 194.16 secs
Using this approach, the Celeron machine (dual core) takes ~180 seconds to reach M9941 and ~270 seconds to reach M11213.
Haskell
import Data.Numbers.Primes (primes)
import Text.Printf (printf)
lucasLehmer :: Int -> Bool
lucasLehmer p = iterate f 4 !! p-2 == 0
where
f b = (b^2 - 2) `mod` m
m = 2^p - 1
main = mapM_ (printf "M %d\n") $ take 20 mersenne
where
mersenne = filter lucasLehmer primes
- Output:
M 3 M 5 M 7 M 13 M 17 M 19 M 31 M 61 M 89 M 107 M 127 M 521 M 607 M 1279 M 2203 M 2281 M 3217 M 4253 M 4423 M 9689
J
I. 1 p: <: 2x ^ i. 1280
2 3 5 7 13 17 19 31 61 89 107 127 521 607 1279
Java
import java.math.BigInteger;
public class MersennePrimes {
private static final int MAX = 20;
private static final BigInteger ONE = BigInteger.ONE;
private static final BigInteger TWO = BigInteger.valueOf(2);
private static boolean isPrime(int n) {
if (n < 2) return false;
if (n % 2 == 0) return n == 2;
if (n % 3 == 0) return n == 3;
int d = 5;
while (d * d <= n) {
if (n % d == 0) return false;
d += 2;
if (n % d == 0) return false;
d += 4;
}
return true;
}
public static void main(String[] args) {
int count = 0;
int p = 2;
while (true) {
BigInteger m = TWO.shiftLeft(p - 1).subtract(ONE);
if (m.isProbablePrime(10)) {
System.out.printf("2 ^ %d - 1\n", p);
if (++count == MAX) break;
}
// obtain next prime, p
do {
p = (p > 2) ? p + 2 : 3;
} while (!isPrime(p));
}
}
}
- Output:
2 ^ 2 - 1 2 ^ 3 - 1 2 ^ 5 - 1 2 ^ 7 - 1 2 ^ 13 - 1 2 ^ 17 - 1 2 ^ 19 - 1 2 ^ 31 - 1 2 ^ 61 - 1 2 ^ 89 - 1 2 ^ 107 - 1 2 ^ 127 - 1 2 ^ 521 - 1 2 ^ 607 - 1 2 ^ 1279 - 1 2 ^ 2203 - 1 2 ^ 2281 - 1 2 ^ 3217 - 1 2 ^ 4253 - 1 2 ^ 4423 - 1
Julia
Julia module Primes
uses Miller-Rabin primality test.
using Primes
mersenne(n::Integer) = convert(typeof(n), 2) ^ n - one(n)
function main(nmax::Integer)
n = ith = zero(nmax)
while ith ≤ nmax
if isprime(mersenne(n))
println("M$n")
ith += 1
end
n += 1
end
end
main(big(20))
- Output:
M2 M3 M5 M7 M13 M17 M19 M31 M61 M89 M107 M127 M521 M607 M1279 M2203 M2281 M3217 M4253 M4423 M9689
Kotlin
This task is similar to the Lucas-Lehmer test task except that you can use whatever method you like to test the primality of the Mersenne numbers. Here, I've chosen to use the JDK's BigInteger.isProbablePrime(certainty) method. The exact algorithm is implementation dependent --- GNU classpath uses only Miller-Rabin, while Oracle JDK uses Miller-Rabin and sometimes adds a Lucas test (this is not the Lucas-Lehmer test).
A 'certainty' parameter of 10 is enough to find the first 20 Mersenne primes but as even this takes about 90 seconds on my modest machine I've not bothered going beyond that.
// version 1.2.10
import java.math.BigInteger
const val MAX = 20
val bigOne = BigInteger.ONE
val bigTwo = 2.toBigInteger()
/* for checking 'small' primes */
fun isPrime(n: Int): Boolean {
if (n < 2) return false
if (n % 2 == 0) return n == 2
if (n % 3 == 0) return n == 3
var d : Int = 5
while (d * d <= n) {
if (n % d == 0) return false
d += 2
if (n % d == 0) return false
d += 4
}
return true
}
fun main(args: Array<String>) {
var count = 0
var p = 2
while (true) {
val m = (bigTwo shl (p - 1)) - bigOne
if (m.isProbablePrime(10)) {
println("2 ^ $p - 1")
if (++count == MAX) break
}
// obtain next prime, p
while(true) {
p = if (p > 2) p + 2 else 3
if (isPrime(p)) break
}
}
}
- Output:
2 ^ 2 - 1 2 ^ 3 - 1 2 ^ 5 - 1 2 ^ 7 - 1 2 ^ 13 - 1 2 ^ 17 - 1 2 ^ 19 - 1 2 ^ 31 - 1 2 ^ 61 - 1 2 ^ 89 - 1 2 ^ 107 - 1 2 ^ 127 - 1 2 ^ 521 - 1 2 ^ 607 - 1 2 ^ 1279 - 1 2 ^ 2203 - 1 2 ^ 2281 - 1 2 ^ 3217 - 1 2 ^ 4253 - 1 2 ^ 4423 - 1
Lua
This checks for primality using a trial division function. The limitation is 'until p == p + 1', meaning that the program will halt when Lua's number type (a 64-bit floating point value) no longer has enough precision to distiguish between one integer and the next.
-- Returns a boolean to show whether x is prime
function isPrime (x)
if x < 2 then return false end
if x < 4 then return true end
if x % 2 == 0 then return false end
for d = 3, math.sqrt(x), 2 do
if x % d == 0 then return false end
end
return true
end
-- Main procedure
local i, p = 0
repeat
i = i + 1
p = 2 ^ i - 1
if isPrime(p) then
print("2 ^ " .. i .. " - 1 = " .. p)
end
until p == p + 1
- Output:
2 ^ 2 - 1 = 3 2 ^ 3 - 1 = 7 2 ^ 5 - 1 = 31 2 ^ 7 - 1 = 127 2 ^ 13 - 1 = 8191 2 ^ 17 - 1 = 131071 2 ^ 19 - 1 = 524287 2 ^ 31 - 1 = 2147483647
Mathematica/Wolfram Language
Mathematica has a built-in function to show the Mersenne primes:
MersennePrimeExponent /@ Range[40]
- Output:
{2, 3, 5, 7, 13, 17, 19, 31, 61, 89, 107, 127, 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}
Alternatively, we can calculate them:
res = {};
Do[
If[PrimeQ[2^i - 1],
AppendTo[res, i];
If[Length[res] == 20,
Break[]
]
]
,
{i, 1, \[Infinity]}
]
res
- Output:
{2, 3, 5, 7, 13, 17, 19, 31, 61, 89, 107, 127, 521, 607, 1279, 2203, 2281, 3217, 4253, 4423}
Nim
Using only standard library
If we want to use only the standard library, we are limited to 64 bits. So we used a simple primality test.
func isOddPrime(n: uint64): bool =
if n == 1: return false
if n mod 3 == 0: return n == 3
var d = 5u
while d * d <= n:
if n mod d == 0: return false
inc d, 2
if n mod d == 0: return false
inc d, 4
result = true
var p = 2u64
for e in 1..63:
if isOddPrime(p - 1):
echo "2^", e, " - 1"
p *= 2
- Output:
2^2 - 1 2^3 - 1 2^5 - 1 2^7 - 1 2^13 - 1 2^17 - 1 2^19 - 1 2^31 - 1 2^61 - 1
Using big integers
The module bignum
provides big integers and a probabilistic primality test. We searched the Mersenne numbers for exponents between 1 and 10_000.
import bignum
var p = newInt(2)
for e in 1..10_000:
if probablyPrime(p - 1, 25) != 0:
echo "2^", e, " - 1"
p *= 2
- Output:
2^2 - 1 2^3 - 1 2^5 - 1 2^7 - 1 2^13 - 1 2^17 - 1 2^19 - 1 2^31 - 1 2^61 - 1 2^89 - 1 2^107 - 1 2^127 - 1 2^521 - 1 2^607 - 1 2^1279 - 1 2^2203 - 1 2^2281 - 1 2^3217 - 1 2^4253 - 1 2^4423 - 1 2^9689 - 1 2^9941 - 1
PARI/GP
LL(p)={
my(m=Mod(4,1<<p-1));
for(i=3,p,m=m^2-2);
m==0
};
forprime(p=2,, if(LL(p), print("2^"p"-1")))
Perl
Since GIMPS went to the trouble of dedicating thousands of CPU years to finding Mersenne primes, we should be kind enough to use the results. The ntheory module front end does this, so the results up to 43 million is extremely fast (4 seconds), and we can reduce this another 10x by only checking primes. After the GIMPS double-checked mark, a Lucas-Lehmer test is done using code similar to Rosetta Code Lucas-Lehmer in C+GMP.
If this is too contrived, we can use Math::Prime::Util::GMP::is_mersenne_prime
instead, which will run the Lucas-Lehmer test on each input. The first 23 Mersenne primes are found in under 15 seconds.
use ntheory qw/forprimes is_mersenne_prime/;
forprimes { is_mersenne_prime($_) && say } 1e9;
- Output:
2 3 5 7 13 17 19 31 61 ...
Phix
with javascript_semantics include mpfr.e atom t0 = time() mpz mp = mpz_init(), bp = mpz_init() integer p = 0, count = 0 constant lim = iff(platform()=JS?14:17) while true do mpz_ui_pow_ui(mp,2,p) mpz_sub_ui(mp,mp,1) if mpz_prime(mp) then string s = iff(time()-t0<1?"":", "&elapsed(time()-t0)) printf(1, "2^%d-1%s\n",{p,s}) count += 1 if count>=lim then exit end if end if while true do p = iff(p>2?p+2:3) mpz_set_si(bp,p) if mpz_prime(bp) then exit end if end while end while {mp,bp} = mpz_free({mp,bp})
- Output:
2^3-1 2^5-1 2^7-1 2^13-1 2^17-1 2^19-1 2^31-1 2^61-1 2^89-1 2^107-1 2^127-1 2^521-1 2^607-1 2^1279-1 2^2203-1, 2.5s 2^2281-1, 2.9s 2^3217-1, 9.5s
PHP
<?php
function is_prime($n) {
if ($n <= 3) {
return $n > 1;
} elseif (($n % 2 == 0) or ($n % 3 == 0)) {
return false;
}
$i = 5;
while ($i * $i <= $n) {
if ($n % $i == 0) {
return false;
}
$i += 2;
if ($n % $i == 0) {
return false;
}
$i += 4;
}
return true;
}
for ($i = 0 ; $i <= 63 ; $i++) {
$pow = pow(2, $i) - 1;
$mersenne = is_prime($pow);
if ($mersenne) {
echo '2 ^ ', $i, ' - 1', PHP_EOL;
}
}
- Output:
2 ^ 2 - 1 2 ^ 3 - 1 2 ^ 5 - 1 2 ^ 7 - 1 2 ^ 13 - 1 2 ^ 17 - 1 2 ^ 19 - 1 2 ^ 31 - 1 2 ^ 61 - 1
PicoLisp
(de **Mod (X Y N)
(let M 1
(loop
(when (bit? 1 Y)
(setq M (% (* M X) N)) )
(T (=0 (setq Y (>> 1 Y)))
M )
(setq X (% (* X X) N)) ) ) )
(de isprime (N)
(cache '(NIL) N
(if (== N 2)
T
(and
(> N 1)
(bit? 1 N)
(let (Q (dec N) N1 (dec N) K 0 X)
(until (bit? 1 Q)
(setq
Q (>> 1 Q)
K (inc K) ) )
(catch 'composite
(do 16
(loop
(setq X
(**Mod
(rand 2 (min (dec N) 1000000000000))
Q
N ) )
(T (or (=1 X) (= X N1)))
(T
(do K
(setq X (**Mod X 2 N))
(when (=1 X) (throw 'composite))
(T (= X N1) T) ) )
(throw 'composite) ) )
(throw 'composite T) ) ) ) ) ) )
(for N 1000
(and
(isprime (dec (** 2 N)))
(prinl "2 \^ " N " - 1") ) )
- Output:
2 ^ 2 - 1 2 ^ 3 - 1 2 ^ 5 - 1 2 ^ 7 - 1 2 ^ 13 - 1 2 ^ 17 - 1 2 ^ 19 - 1 2 ^ 31 - 1 2 ^ 61 - 1 2 ^ 89 - 1 2 ^ 107 - 1 2 ^ 127 - 1 2 ^ 521 - 1 2 ^ 607 - 1
Pike
int power = 1;
while(power++) {
int candidate = 2->pow(power)-1;
if( candidate->probably_prime_p() )
write("2 ^ %d - 1\n", power);
}
Output:
2 ^ 2 - 1 2 ^ 3 - 1 2 ^ 5 - 1 2 ^ 7 - 1 2 ^ 13 - 1 2 ^ 17 - 1 2 ^ 19 - 1 2 ^ 31 - 1 2 ^ 61 - 1
Prolog
Lucas-Lehmer test, works with SWI Prolog
lucas_lehmer_seq(M, L) :-
lazy_list(ll_iter, 4-M, L).
ll_iter(S-M, T-M, T) :-
T is ((S*S) - 2) mod M.
drop(N, Lz1, Lz2) :-
append(Pfx, Lz2, Lz1), length(Pfx, N), !.
mersenne_prime(2).
mersenne_prime(P) :-
P > 2,
M is (1 << P) - 1,
lucas_lehmer_seq(M, Residues),
Skip is P - 3, drop(Skip, Residues, [R|_]),
R =:= 0.
- Output:
?- findall(X, (between(1, 1000, X), mersenne_prime(X)), L), write(L). [2,3,5,7,13,17,19,31,61,89,107,127,521,607] L = [2, 3, 5, 7, 13, 17, 19, 31, 61|...].
Python
import random
#Take from https://www.codeproject.com/Articles/691200/%2FArticles%2F691200%2FPrimality-test-algorithms-Prime-test-The-fastest-w
def MillerRabinPrimalityTest(number):
'''
because the algorithm input is ODD number than if we get
even and it is the number 2 we return TRUE ( spcial case )
if we get the number 1 we return false and any other even
number we will return false.
'''
if number == 2:
return True
elif number == 1 or number % 2 == 0:
return False
''' first we want to express n as : 2^s * r ( were r is odd ) '''
''' the odd part of the number '''
oddPartOfNumber = number - 1
''' The number of time that the number is divided by two '''
timesTwoDividNumber = 0
''' while r is even divid by 2 to find the odd part '''
while oddPartOfNumber % 2 == 0:
oddPartOfNumber = oddPartOfNumber / 2
timesTwoDividNumber = timesTwoDividNumber + 1
'''
since there are number that are cases of "strong liar" we
need to check more then one number
'''
for time in range(3):
''' choose "Good" random number '''
while True:
''' Draw a RANDOM number in range of number ( Z_number ) '''
randomNumber = random.randint(2, number)-1
if randomNumber != 0 and randomNumber != 1:
break
''' randomNumberWithPower = randomNumber^oddPartOfNumber mod number '''
randomNumberWithPower = pow(randomNumber, oddPartOfNumber, number)
''' if random number is not 1 and not -1 ( in mod n ) '''
if (randomNumberWithPower != 1) and (randomNumberWithPower != number - 1):
# number of iteration
iterationNumber = 1
''' while we can squre the number and the squered number is not -1 mod number'''
while (iterationNumber <= timesTwoDividNumber - 1) and (randomNumberWithPower != number - 1):
''' squre the number '''
randomNumberWithPower = pow(randomNumberWithPower, 2, number)
# inc the number of iteration
iterationNumber = iterationNumber + 1
'''
if x != -1 mod number then it because we did not found strong witnesses
hence 1 have more then two roots in mod n ==>
n is composite ==> return false for primality
'''
if (randomNumberWithPower != (number - 1)):
return False
''' well the number pass the tests ==> it is probably prime ==> return true for primality '''
return True
# Main
MAX = 20
p = 2
count = 0
while True:
m = (2 << (p - 1)) - 1
if MillerRabinPrimalityTest(m):
print "2 ^ {} - 1".format(p)
count = count + 1
if count == MAX:
break
# obtain next prime, p
while True:
p = p + 2 if (p > 2) else 3
if MillerRabinPrimalityTest(p):
break
print "done"
- Output:
2 ^ 2 - 1 2 ^ 3 - 1 2 ^ 5 - 1 2 ^ 7 - 1 2 ^ 13 - 1 2 ^ 17 - 1 2 ^ 19 - 1 2 ^ 31 - 1 2 ^ 61 - 1 2 ^ 89 - 1 2 ^ 107 - 1 2 ^ 127 - 1 2 ^ 521 - 1 2 ^ 607 - 1 2 ^ 1279 - 1 2 ^ 2203 - 1 2 ^ 2281 - 1 2 ^ 3217 - 1 2 ^ 4253 - 1 2 ^ 4423 - 1 done
Raku
(formerly Perl 6)
We already have a multitude of tasks that demonstrate how to find Mersenne primes; Prime decomposition, Primality by trial division, Trial factoring of a Mersenne number, Lucas-Lehmer test, Miller–Rabin primality_test, etc. that all have Raku entries. I'm not sure what I could add here that would be useful.
Hmmm.
Create code that will list all of the Mersenne primes until some limitation is reached.
It doesn't specify to calculate them, only to list them; why throw away all of the computer millenia of processing power that the GIMPS has invested?
use HTTP::UserAgent;
use Gumbo;
my $table = parse-html(HTTP::UserAgent.new.get('https://www.mersenne.org/primes/').content, :TAG<table>);
say 'All known Mersenne primes as of ', Date(now);
say 'M', ++$, ": 2$_ - 1"
for $table[1]».[*][0][*].comb(/'exp_lo='\d+/)».subst(/\D/, '',:g)
.trans([<0123456789>.comb] => [<⁰¹²³⁴⁵⁶⁷⁸⁹>.comb]).words;
- Output:
All known Mersenne primes as of 2018-12-21 M1: 2² - 1 M2: 2³ - 1 M3: 2⁵ - 1 M4: 2⁷ - 1 M5: 2¹³ - 1 M6: 2¹⁷ - 1 M7: 2¹⁹ - 1 M8: 2³¹ - 1 M9: 2⁶¹ - 1 M10: 2⁸⁹ - 1 M11: 2¹⁰⁷ - 1 M12: 2¹²⁷ - 1 M13: 2⁵²¹ - 1 M14: 2⁶⁰⁷ - 1 M15: 2¹²⁷⁹ - 1 M16: 2²²⁰³ - 1 M17: 2²²⁸¹ - 1 M18: 2³²¹⁷ - 1 M19: 2⁴²⁵³ - 1 M20: 2⁴⁴²³ - 1 M21: 2⁹⁶⁸⁹ - 1 M22: 2⁹⁹⁴¹ - 1 M23: 2¹¹²¹³ - 1 M24: 2¹⁹⁹³⁷ - 1 M25: 2²¹⁷⁰¹ - 1 M26: 2²³²⁰⁹ - 1 M27: 2⁴⁴⁴⁹⁷ - 1 M28: 2⁸⁶²⁴³ - 1 M29: 2¹¹⁰⁵⁰³ - 1 M30: 2¹³²⁰⁴⁹ - 1 M31: 2²¹⁶⁰⁹¹ - 1 M32: 2⁷⁵⁶⁸³⁹ - 1 M33: 2⁸⁵⁹⁴³³ - 1 M34: 2¹²⁵⁷⁷⁸⁷ - 1 M35: 2¹³⁹⁸²⁶⁹ - 1 M36: 2²⁹⁷⁶²²¹ - 1 M37: 2³⁰²¹³⁷⁷ - 1 M38: 2⁶⁹⁷²⁵⁹³ - 1 M39: 2¹³⁴⁶⁶⁹¹⁷ - 1 M40: 2²⁰⁹⁹⁶⁰¹¹ - 1 M41: 2²⁴⁰³⁶⁵⁸³ - 1 M42: 2²⁵⁹⁶⁴⁹⁵¹ - 1 M43: 2³⁰⁴⁰²⁴⁵⁷ - 1 M44: 2³²⁵⁸²⁶⁵⁷ - 1 M45: 2³⁷¹⁵⁶⁶⁶⁷ - 1 M46: 2⁴²⁶⁴³⁸⁰¹ - 1 M47: 2⁴³¹¹²⁶⁰⁹ - 1 M48: 2⁵⁷⁸⁸⁵¹⁶¹ - 1 M49: 2⁷⁴²⁰⁷²⁸¹ - 1 M50: 2⁷⁷²³²⁹¹⁷ - 1 M51: 2⁸²⁵⁸⁹⁹³³ - 1
REXX
This REXX version (using a 32-bit Regina REXX interpreter) will find those Mersenne primes which are less than
8 million decimal digits (which would be M43).
/*REXX program uses exponent─and─mod operator to test possible Mersenne numbers. */
do j=1; /*process a range, or run out of time.*/
if \isPrime(j) then iterate /*if J isn't a prime, keep plugging.*/
r= testMer(j) /*If J is prime, give J the 3rd degree.*/
if r==0 then say right('M'j, 10) "──────── is a Mersenne prime."
else say right('M'j, 50) "is composite, a factor:" r
end /*j*/
exit /*stick a fork in it, we're all done. */
/*──────────────────────────────────────────────────────────────────────────────────────*/
isPrime: procedure; parse arg x; if wordpos(x, '2 3 5 7') \== 0 then return 1
if x<11 then return 0; if x//2 == 0 | x//3 == 0 then return 0
do j=5 by 6; if x//j == 0 | x//(j+2) == 0 then return 0
if j*j>x then return 1 /*◄─┐ ___ */
end /*j*/ /* └─◄ Is j>√ x ? Then return 1*/
/*──────────────────────────────────────────────────────────────────────────────────────*/
iSqrt: procedure; parse arg x; #= 1; r= 0; do while #<=x; #=#*4; end
do while #>1; #=#%4; _= x-r-#; r= r%2; if _>=0 then do; x=_; r=r+#; end
end /*while*/ /*iSqrt ≡ integer square root.*/
return r /*───── ─ ── ─ ─ */
/*──────────────────────────────────────────────────────────────────────────────────────*/
testMer: procedure; parse arg x; p =2**x /* [↓] do we have enough digits?*/
$$=x2b( d2x(x) ) + 0
if pos('E',p)\==0 then do; parse var p "E" _; numeric digits _+2; p=2**x; end
!.=1; !.1=0; !.7=0 /*array used for a quicker test. */
R=iSqrt(p) /*obtain integer square root of P*/
do k=2 by 2; q=k*x + 1 /*(shortcut) compute value of Q. */
m=q // 8 /*obtain the remainder when ÷ 8. */
if !.m then iterate /*M must be either one or seven.*/
parse var q '' -1 _; if _==5 then iterate /*last digit a five?*/
if q// 3==0 then iterate /* ÷ by three? */
if q// 7==0 then iterate /* " " seven? */
if q//11==0 then iterate /* " " eleven?*/
/* ____ */
if q>R then return 0 /*Is q>√2**x ? A Mersenne prime*/
sq=1; $=$$ /*obtain binary version from $. */
do until $==''; sq=sq*sq
parse var $ _ 2 $ /*obtain 1st digit and the rest. */
if _ then sq=(sq+sq) // q
end /*until*/
if sq==1 then return q /*Not a prime? Return a factor.*/
end /*k*/
Ring
# Project : Mersenne primes
n = 0
while true
n = n +1
if isprime(pow(2,n)-1) = 1
see n + nl
ok
end
func isprime num
if (num <= 1) return 0 ok
if (num % 2 = 0) and num != 2 return 0 ok
for i = 3 to floor(num / 2) -1 step 2
if (num % i = 0) return 0 ok
next
return 1
Output:
2 3 5 7 13 17 19
RPL
This version use binary integers
RPL code | Comment |
---|---|
≪ / LAST ROT * - #0 == ≫ 'BDIV?' STO ≪ IF DUP #3 ≤ THEN #2 / B→R ELSE IF DUP #2 BDIV? OVER #3 BDIV? OR THEN DROP 0 ELSE DUP B→R √ R→B → a maxd ≪ a #2 #5 WHILE a OVER BDIV? NOT OVER maxd ≤ AND REPEAT OVER + #6 ROT - SWAP END ≫ SWAP DROP BDIV? NOT END END ≫ 'BPRIM?' STO ≪ {} 1 32 FOR n #1 1 n START SL NEXT #1 - IF BPRIM? THEN + END NEXT ≫ EVAL |
( #a #b -- boolean ) ( #a -- boolean ) return 1 if a is 2 or 3 if 2 or 3 divides a return 0 else store a and root(a) initialize stack with a i d while d does not divide a and d <= root(a) i = 6 - i which modifies 2 into 4 and viceversa convert stack status into result Let's loop Generate M(n) test primality |
- Output:
{ 2 3 5 7 13 17 19 31 }
Ruby
require 'openssl'
(0..).each{|n| puts "2**#{n} - 1" if OpenSSL::BN.new(2**n -1).prime? }
- Output:
2**2 - 1 2**3 - 1 2**5 - 1 2**7 - 1 2**13 - 1 2**17 - 1 2**19 - 1 2**31 - 1 2**61 - 1 2**89 - 1 2**107 - 1 2**127 - 1 2**521 - 1 2**607 - 1 2**1279 - 1 2**2203 - 1 2**2281 - 1 2**3217 - 1 2**4253 - 1 ^Ctest2.rb:7:in `prime?': Interrupt
Scala
object MersennePrimes extends App {
val primes = primeSieve(LazyList.from(2))
val upbPrime = 9941
def primeSieve(s: LazyList[Int]): LazyList[Int] =
s.head #:: primeSieve(s.tail filter {
_ % s.head != 0
})
def mersenne(p: Int): BigInt = (BigInt(2) pow p) - 1
def s(mp: BigInt, p: Int): BigInt = {
if (p == 1) 4 else ((s(mp, p - 1) pow 2) - 2) % mp
}
println(s"Finding Mersenne primes in M[2..$upbPrime]")
((primes takeWhile (_ <= upbPrime)).map { p => (p, mersenne(p)) }
map { p => if (p._1 == 2) (p, 0) else (p, s(p._2, p._1 - 1)) } filter {
_._2 == 0
})
.foreach { p =>
println(s"prime M${(p._1)._1}: " + {
if ((p._1)._1 < 200) (p._1)._2 else s"(${(p._1)._2.toString.length} digits)"
})
}
println("That's All Folks!")
}
Sidef
Uses the is_mersenne_prime() function from Math::Prime::Util::GMP.
for p in (^Inf -> lazy.grep { .is_mersenne_prime }) {
say "2^#{p} - 1"
}
- Output:
2^2 - 1 2^3 - 1 2^5 - 1 2^7 - 1 2^13 - 1 2^17 - 1 2^19 - 1 2^31 - 1 2^61 - 1 2^89 - 1 2^107 - 1 2^127 - 1 2^521 - 1 2^607 - 1 2^1279 - 1 2^2203 - 1 2^2281 - 1 2^3217 - 1 2^4253 - 1 2^4423 - 1 2^9689 - 1 2^9941 - 1 ^C sidef mersenne.sf 12.47s user 0.02s system 99% cpu 12.495 total
Transd
#lang transd
MainModule: {
_start: (λ locals: curexp 1 maxexp 1000 base BigLong(2)
(while (< curexp maxexp)
(+= curexp 1)
(if (not (is-prime BigLong(curexp))) continue)
(with cand (- (pow base curexp) 1)
(if (is-probable-prime cand 10)
(lout curexp " : " cand ))
) ) )
}
- Output:
2 : 3 3 : 7 5 : 31 7 : 127 13 : 8191 17 : 131071 19 : 524287 31 : 2147483647 61 : 2305843009213693951 89 : 618970019642690137449562111 107 : 162259276829213363391578010288127 127 : 170141183460469231731687303715884105727 521 : 6864797660130609714981900799081393217269435300143305409394463459185543183397656052122559640661454554977296311391480858037121987999716643812574028291115057151 607 : 531137992816767098689588206552468627329593117727031923199444138200403559860852242739162502265229285668889329486246501015346579337652707239409519978766587351943831270835393219031728127
Visual Basic .NET
Imports System.Numerics
Module Module1
Function Sqrt(x As BigInteger) As BigInteger
If x < 0 Then
Throw New ArgumentException("Negative argument.")
End If
If x < 2 Then
Return x
End If
Dim y = x / 2
While y > (x / y)
y = ((x / y) + y) / 2
End While
Return y
End Function
Function IsPrime(bi As BigInteger) As Boolean
If bi < 2 Then
Return False
End If
If bi Mod 2 = 0 Then
Return bi = 2
End If
If bi Mod 3 = 0 Then
Return bi = 3
End If
If bi Mod 5 = 0 Then
Return bi = 5
End If
If bi Mod 7 = 0 Then
Return bi = 7
End If
If bi Mod 11 = 0 Then
Return bi = 11
End If
If bi Mod 13 = 0 Then
Return bi = 13
End If
If bi Mod 17 = 0 Then
Return bi = 17
End If
If bi Mod 19 = 0 Then
Return bi = 19
End If
Dim limit = Sqrt(bi)
Dim test As BigInteger = 23
While test < limit
If bi Mod test = 0 Then
Return False
End If
test += 2
If bi Mod test = 0 Then
Return False
End If
test += 4
End While
Return True
End Function
Sub Main()
Const MAX = 9
Dim pow = 2
Dim count = 0
While True
If IsPrime(pow) Then
Dim p = BigInteger.Pow(2, pow) - 1
If IsPrime(p) Then
Console.WriteLine("2 ^ {0} - 1", pow)
count += 1
If count >= MAX Then
Exit While
End If
End If
End If
pow += 1
End While
End Sub
End Module
- Output:
2 ^ 2 - 1 2 ^ 3 - 1 2 ^ 5 - 1 2 ^ 7 - 1 2 ^ 13 - 1 2 ^ 17 - 1 2 ^ 19 - 1 2 ^ 31 - 1 2 ^ 61 - 1
Wren
Wren-CLI
A bit slow so limited to first 14 Mersenne primes.
import "./math" for Int
import "./big" for BigInt
var MAX = 14
System.print("The first %(MAX) Mersenne primes are:")
var count = 0
var p = 2
while (true) {
var m = (BigInt.one << p) - 1
if (m.isProbablePrime(10)) {
System.print("2 ^ %(p) - 1")
count = count + 1
if (count == MAX) break
}
while (true) {
p = (p > 2) ? p + 2 : 3
if (Int.isPrime(p)) break
}
}
- Output:
The first 14 Mersenne primes are: 2 ^ 2 - 1 2 ^ 3 - 1 2 ^ 5 - 1 2 ^ 7 - 1 2 ^ 13 - 1 2 ^ 17 - 1 2 ^ 19 - 1 2 ^ 31 - 1 2 ^ 61 - 1 2 ^ 89 - 1 2 ^ 107 - 1 2 ^ 127 - 1 2 ^ 521 - 1 2 ^ 607 - 1
Embedded (GMP)
This finds the first 23 Mersenne primes in about 172 seconds which is virtually the same as the Go non-concurrent version using their GMP plug-in when run on my machine.
import "./math" for Int
import "./gmp" for Mpz
var MAX = 23
System.print("The first %(MAX) Mersenne primes are:")
var count = 0
var p = 2
while (true) {
var m = Mpz.one.lsh(p).sub(1)
if (m.probPrime(15) > 0) {
System.print("2 ^ %(p) - 1")
count = count + 1
if (count == MAX) break
}
while (true) {
p = (p > 2) ? p + 2 : 3
if (Int.isPrime(p)) break
}
}
- Output:
The first 23 Mersenne primes are: 2 ^ 2 - 1 2 ^ 3 - 1 2 ^ 5 - 1 2 ^ 7 - 1 2 ^ 13 - 1 2 ^ 17 - 1 2 ^ 19 - 1 2 ^ 31 - 1 2 ^ 61 - 1 2 ^ 89 - 1 2 ^ 107 - 1 2 ^ 127 - 1 2 ^ 521 - 1 2 ^ 607 - 1 2 ^ 1279 - 1 2 ^ 2203 - 1 2 ^ 2281 - 1 2 ^ 3217 - 1 2 ^ 4253 - 1 2 ^ 4423 - 1 2 ^ 9689 - 1 2 ^ 9941 - 1 2 ^ 11213 - 1
XPL0
func IsPrime(N); \Return 'true' if N is prime
int N, I;
[if N <= 2 then return N = 2;
if (N&1) = 0 then \even >2\ return false;
for I:= 3 to sqrt(N) do
[if rem(N/I) = 0 then return false;
I:= I+1;
];
return true;
];
int N;
[for N:= 1 to 31 do
if IsPrime(1<<N-1) then
[Text(0, "2^^"); IntOut(0, N); Text(0, " - 1");
CrLf(0);
];
]
- Output:
2^2 - 1 2^3 - 1 2^5 - 1 2^7 - 1 2^13 - 1 2^17 - 1 2^19 - 1 2^31 - 1
zkl
Uses libGMP (GNU MP Bignum Library) and its Miller-Rabin probabilistic primality testing algorithm.
var [const] BN=Import.lib("zklBigNum"); // libGMP
fcn mprimes{
n,m := BN(2),0;
foreach e in ([2..]){
n,m = n.shiftLeft(1), n-1;
if(m.probablyPrime()) println("2^%d - 1".fmt(e));
}
}()
// gets rather slow after M(4423)
- Output:
2^2 - 1 2^3 - 1 2^5 - 1 2^7 - 1 2^13 - 1 2^17 - 1 2^19 - 1 2^31 - 1 2^61 - 1 2^89 - 1 2^107 - 1 2^127 - 1 2^521 - 1 2^607 - 1 2^1279 - 1 2^2203 - 1 2^2281 - 1 2^3217 - 1 2^4253 - 1 2^4423 - 1 2^9689 - 1 2^9941 - 1 2^11213 - 1 ...
- Draft Programming Tasks
- Prime Numbers
- 11l
- ALGOL 60
- ALGOL 68
- ALGOL W
- AppleScript
- Arturo
- AWK
- BASIC
- BASIC256
- FreeBASIC
- Yabasic
- C
- GMP
- C++
- C sharp
- D
- Delphi
- SysUtils,StdCtrls
- EasyLang
- F Sharp
- Factor
- Fortran
- Frink
- Go
- Haskell
- J
- Java
- Julia
- Kotlin
- Lua
- Mathematica
- Wolfram Language
- Nim
- Bignum
- PARI/GP
- Perl
- Ntheory
- Phix
- Phix/mpfr
- PHP
- PicoLisp
- Pike
- Prolog
- Python
- Raku
- REXX
- Ring
- RPL
- Ruby
- Scala
- Sidef
- Transd
- Visual Basic .NET
- Wren
- Wren-math
- Wren-big
- Wren-gmp
- XPL0
- Zkl