Mersenne primes

From Rosetta Code
Mersenne primes is a draft programming task. It is not yet considered ready to be promoted as a complete task, for reasons that should be found in its talk page.
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



11l

Translation of: D
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

Works with: A60
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

Works with: ALGOL 68G version Any - tested with release 2.8.3.win32
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

Translation of: FreeBASIC
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

Translation of: FreeBASIC
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
Library: GMP

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++

Translation of: 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

Works with: Delphi version 6.0


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#

Translation of: C#
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

Translation of: C
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

Translation of: Kotlin
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

Works with: Julia version 0.6

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

Library: bignum

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.

Library: ntheory
use ntheory qw/forprimes is_mersenne_prime/;
forprimes { is_mersenne_prime($_) && say } 1e9;
Output:
2
3
5
7
13
17
19
31
61
...

Phix

Library: Phix/mpfr
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

Translation of: Java
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)

Works with: Rakudo version 2018.01

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

Works with: Halcyon Calc version 4.2.7
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

Translation of: C#
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

Library: Wren-math
Library: Wren-big

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)

Library: 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.

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

Library: GMP

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