Lucas-Lehmer test: Difference between revisions

{{trans|D}}

 

I p < 2 | p % 2 == 0

R p == 2

L(p) 2..2299

I isMersennePrime(p)

 

{{out}}

=={{header|360 Assembly}}==

For maximum compatibility, this program uses only the basic instruction set.

LUCASLEH CSECT

USING LUCASLEH,R12

LTORG

YREGS

{{out}}

<pre>M002

 

=={{header|Ada}}==

with Ada.Integer_Text_Io; use Ada.Integer_Text_Io;

 

end if;

end loop;

{{Out}}

Mersenne primes:

Because of the very large numbers computed,

the mapm binding is used to calculate with arbitrary precision.

 

mapm.xdigits(100);

write('M' & i & ' ')

fi

produces:

 

=={{header|ALGOL 68}}==

{{works with|ALGOL 68G|Any - tested with release [http://sourceforge.net/projects/algol68/files/algol68g/algol68g-1.18.0/algol68g-1.18.0-9h.tiny.el5.centos.fc11.i386.rpm/download 1.18.0-9h.tiny]}}

{{wont work with|ELLA ALGOL 68|Any (with appropriate job cards) - tested with release [http://sourceforge.net/projects/algol68/files/algol68toc/algol68toc-1.8.8d/algol68toc-1.8-8d.fc9.i386.rpm/download 1.8-8d] - due to extensive use of '''format'''[ted] ''transput''}}

 

PROC is prime = ( INT p )BOOL:

count <= upb count

DO SKIP OD

{{Out}}

<pre>

=={{header|ARM Assembly}}==

{{works with|as|Raspberry Pi}}

/* ARM assembly Raspberry PI */

/* program lucaslehmer.s */

iMagicNumber: .int 0xCCCCCCCD

 

{{Out}}

<pre>

=={{header|Arturo}}==

 

if p=2 -> return true

mp: dec shl 1 p

print "Mersenne primes:"

mersennes: select 2..32 'x -> and? prime? x mersenne? x

 

{{out}}

 

=={{header|AWK}}==

# syntax: GAWK -f LUCAS-LEHMER_TEST.AWK

# converted from Pascal

exit(0)

}

{{Out}}

<pre>

</pre>

 

{{works with|BBC BASIC for Windows}}

Using its native arithmetic BBC BASIC can only test up to M23.

PRINT "Mersenne Primes:"

FOR p% = 2 TO 23

sn -= (mp * INT(sn / mp))

NEXT

{{Out}}

<pre>

M19

</pre>

 

=={{header|BCPL}}==

{{trans|Zig}}

Uses the 64 bit version

GET "libhdr"

 

RESULTIS 0

}

{{Out}}

<pre>

If a number cannot be factorized, (either because it is prime or because it is to great to fit in a computer word) the root expression doesn't change much.

For example, the expression <code>13^(13^-1)</code> becomes <code>13^1/13</code>, and this matches the pattern <code>13^%</code>.

& ( time

= ( print

)

& done

{{Out}} (after 4.5 hours):

<pre>0,00 0,00 s: M3 is PRIME!

1101,12 10491,08 s: M9941 is PRIME!

5618,45 16109,53 s: M11213 is PRIME!</pre>

 

=={{header|C}}==

 

{{libheader|GMP}}

#include <stdlib.h>

#include <limits.h>

printf("\n");

return 0;

 

{{out}}

{{works with|C99}}

Compiler options: gcc -std=c99 -lm Lucas-Lehmer_test.c -o Lucas-Lehmer_test<br>

#include <stdio.h>

#include <limits.h>

}

printf("\n");

 

{{Out}}

{{works with|Visual Studio|2010}}

{{works with|.NET Framework|4.0}}

using System.Collections.Generic;

using System.Numerics;

}

}

 

{{out}} (Run only to 11213)

The mod function, (<code>%</code>) has a computation cost equivalent to the divide operation. In this case, a combination of ands, shifts and adds can replace the mod function. Another change is creating the list of candidate Mersenne numbers in descending order, the point being to start the more time consuming calculations first. This avoids a long calculation occurring by itself at the end of the <code>Parallel.For</code> queue. Also added trial division step, translated from the '''Rust''' and '''C''' versions.

 

using System.Collections.Generic;

using System.Numerics;

if (System.Diagnostics.Debugger.IsAttached) Console.ReadLine();

}

{{out}}

<pre>M2 M3 M5 M7 M13 M17 M19 M31 M61 M89 M107 M127 M521 M607 M1279 M2203 M2281 M3217 M4253 M4423 M9689 M9941 M11213

{{libheader|GMP}}

 

#include <gmpxx.h>

 

static bool is_mersenne_prime(mpz_class p)

{

return true;

{

 

}

int main()

{

}

}

 

{{out}} (Incomplete; It takes a long time.)

 

=={{header|Clojure}}==

(cond (< i 4) (>= i 2)

(zero? (rem i 2)) false

(recur (inc n) (rem (- (* s s) 2) mp)))))))

 

{{Out}}

<pre>

=={{header|CoffeeScript}}==

Coffee Script is really hampered by its lack of full syntactic support for JavaScript generators. The loop to collect Mersenne numbers must be done in imperative style, rather than a more functional style, when using the infinite lazy prime generator.

sorenson = require('sieve').primes # Sorenson's extensible sieve from task: Extensible Prime Number Generator

 

# Test if 2^n-1 is a Mersenne prime.

#

isMersennePrime = (p) ->

 

console.log "Some Mersenne primes: #{"M" + String p for p in mersennes}"

{{Out}}

<pre>

=={{header|Common Lisp}}==

{{trans|Clojure}}

(defun or-f (&optional a b) (or a b));necessary for reduce, as 'or' is implemented as a macro

 

 

(princ (remove-if-not #'mersenne-p (remove-if-not #'prime-p (loop for i to 5000 collect i))))

{{out}}

<pre>

=={{header|Crystal}}==

{{trans|Ruby}}

 

def is_prime(n) # P3 Prime Generator primality test

break if count >= upb_count

end

 

{{out}}

=={{header|D}}==

{{trans|Python}}

 

bool isPrime(in uint p) pure nothrow @safe @nogc {

stdout.flush;

}

{{out}}

<pre>M2 M3 M5 M7 M13 M17 M19 M31 M61 M89 M107 M127 M521 M607 M1279 M2203 M2281 </pre>

With p up to 10_000 it prints:

<pre>M2 M3 M5 M7 M13 M17 M19 M31 M61 M89 M107 M127 M521 M607 M1279 M2203 M2281 M3217 M4253 M4423 M9941 </pre>

 

=={{header|DWScript}}==

Using Integer type, which is 64bit, limits the search to M31.

var

i, s, m_p : Integer;

Print(' M'+IntToStr(p));

end;

{{Out}}

<pre> M2 M3 M5 M7 M13 M17 M19 M31</pre>

 

=={{header|EchoLisp}}==

(require 'bigint)

(define (mersenne-prime? odd-prime: p)

 

→ M3 M5 M7 M13 M17 M19 M31 M61 M89 M107 M127 M521 M607 M1279 M2203 M2281

 

=={{header|Elixir}}==

{{trans|Erlang}}

use Bitwise

def test do

end

 

 

{{out}}

 

=={{header|Erlang}}==

-export([main/0]).

 

 

s(MP,1) -> 4 rem MP;

 

In 3 seconds will print

=={{header|ERRE}}==

With native arithmetic up to 23: for bigger numbers you must use MULPREC program.

 

!$DOUBLE

END FOR

END PROGRAM

{{out}}

<pre>Mersenne Primes:

=={{header|F Sharp|F#}}==

Simple arbitrary-precision version:

if n = 1 then 4I % mp else ((s mp (n - 1)) ** 2 - 2I) % mp

 

[ for p in 2..47 do

if p = 2 || s ((1I <<< p) - 1I) (p - 1) = 0I then

 

Tail-recursive version:

if exponent <= 1 then failwith "Exponent must be >= 2"

let prime = 2I ** exponent - 1I;

LucasLehmer 0 4I = 0I

 

Version using library folding function (way shorter and faster than the above):

if exponent <= 1 then failwith "Exponent must be >= 2"

let prime = 2I ** exponent - 1I;

LucasLehmer = 0I

 

=={{header|Factor}}==

sequences ;

 

47 [1,b] [ lucas-lehmer ] filter

"Mersenne primes:" print

{{out}}

<pre>

 

=={{header|Forth}}==

1+ 2 do

4 i 2 <> * abs swap 1+ dup + 1- swap

;

 

==Alternate version to handle 64 and 128 bit integers.==

18 constant π-64 \ count of primes < 64

31 constant π-128 \ count of primes < 128

if 'M emit i c@ . then

loop ;

{{Out}}

<pre>

{{works with|Fortran|90 and later}}

Only Mersenne number with prime exponent can be themselves prime but for the small numbers used in this example it was not worth the effort to include this check. As the size of the exponent increases this becomes more important.

IMPLICIT NONE

 

END DO

 

=={{header|FreeBASIC}}==

===Native types for Mersenne primes <= M63===

' compile with: fbc -s console

 

Print : Print "hit any key to end program"

Sleep

{{out}}

<pre> M3 M5 M7 M13 M17 M19 M31 M61</pre>

==={{libheader|GMP}}===

Uses the trick from the '''C''' entry to avoid the slow Mod

' compile with: fbc -s console

 

Print : Print "hit any key to end program"

Sleep

{{out}}

<pre>M3 M5 M7 M13 M17

Frink's <CODE>isPrime</CODE> function automatically detects numbers of the form 2ⁿ-1 and performs a Lucas-Lehmer test on them, including testing if n is prime, which is sufficient to prove primality for this form.

 

for n = primes[]

if isPrime[2^n-1]

println[n]

 

=={{header|FunL}}==

if p == 2 then return true

 

for p <- primes().filter( mersenne ).take( 20 )

 

{{out}}

 

=={{header|GAP}}==

local i, m, s;

if n = 2 then

 

Filtered([1 .. 2000], LucasLehmer);

 

=={{header|Go}}==

Processing the first list indicates that the test works. Processing the second shows it working on some larger numbers.

 

import (

}

}

{{out}}

<pre>

{{works with|GHC|6.8.2}}

 

where

 

lucasLehmer p = s (2^p-1) (p-1) == 0

 

It is pointed out on the [[Sieve of Eratosthenes]] page that the following "sieve" is inefficient. Nonetheless it takes very little time compared to the Lucas-Lehmer test itself.

It takes about 30 minutes to get up to:

<pre>

 

=={{header|HicEst}}==

DO exponent = 2, 31

IF(exponent > 2) s = 4

ENDDO

 

 

=={{header|J}}==

We use arbitrary-precision integers in order to be able to test any arbitrary prime.

 

public class Mersenne

{

System.out.println();

}

{{Out}} (after about eight hours):

<pre>

In JavaScript we using BigInt ( numbers with 'n' suffix ) - so we can use really big numbers

 

////////// In JavaScript we don't have sqrt for BigInt - so here is implementation

function newtonIteration(n, x0) {

}

console.log(`... Took: ${Date.now()-tm} ms`);

{{Out}}

<pre>

Output includes the length of the decimal representation of the Mersenne prime.

 

def power($b): . as $in | reduce range(0;$b) as $i (1; . * $in);

 

| "M\($i):\($mp|tostring|length)" );

 

{{out}}

<pre>

</pre>

=={{header|Julia}}==

using Primes

 

 

getmersenneprimes(50)

{{output}}<pre>

M2, cumulative time elapsed: 0.019999980926513672 seconds

=={{header|Kotlin}}==

In view of the Java result, I've set the program to stop at M4423 so it will run in a reasonable time (about 85 seconds) on a typical laptop:

 

import java.math.BigInteger

}

}

 

{{out}}

=={{header|langur}}==

{{trans|D}}

 

 

if .p == 2: return true

if not .isPrime(.p): return false

 

val .mp = 2 ^ .p - 1

.s = (.s ^ 2 - 2) rem .mp

} == 0

}

 

 

=={{header|Mathematica}}/{{header|Wolfram Language}}==

This version is very speedy and is bounded.

For[i = 1; s = 4, i <= p - 2, i++, s = Mod[s^2 - 2, M]];

If[s == 0, "M" <> ToString@p, p], {p,

Prime /@ Range[300]}], StringQ]

 

 

This version is unbounded (and timed):

For[p = 2, True, p = NextPrime[p], M = 2^p - 1;

For[i = 1; s = 4, i <= p - 2, i++, s = Mod[s^2 - 2, M]];

If[s == 0, Print["M" <> ToString@p]]]

 

I'll see what this gets.

MATLAB suffers from a lack of an arbitrary precision math (bignums) library.

It also doesn't have great support for 64-bit integer arithmetic...or at least MATLAB 2007 doesn't. So, the best precision we have is doubles; therefore, this script can only find up to M19 and no greater.

 

function isPrime = lucasLehmerTest(thePrime)

mersennesPrime = (2.^mNumber) - 1;

 

{{Out}}

 

mNumber =

mersennesPrime =

 

 

=={{header|Maxima}}==

if not primep(p) then false elseif p = 2 then true else

(s: 4,

 

sublist(makelist(i, i, 1, 200), lucas_lehmer);

 

=={{header|Modula-3}}==

Modula-3 uses L as the literal for <tt>LONGINT</tt>.

 

IMPORT IO, Fmt, Long;

END;

IO.Put("\n");

{{Out}}

<pre>M2 M3 M5 M7 M13 M17 M19 M31 </pre>

 

=={{header|Nim}}==

 

proc isPrime(a: int): bool =

if isPrime(p) and isMersennePrime(p):

stdout.write " M",p

{{Out}}

<pre> Mersenne primes:

=={{header|Oz}}==

Oz's multiple precision number system use GMP core.

functor

import

{Application.exit 0}

 

=={{header|PARI/GP}}==

my(m=Mod(4,1<<p-1));

for(i=3,p,m=m^2-2);

if(LL(p), print("2^"p"-1"))

)

 

=={{header|Pascal}}==

int64 is good enough up to M31:

var

s, n: int64;

writeln('M', exponent, ' is PRIME!');

end;

{{Out}}

<pre>:> ./LucasLehmer

=={{header|Perl}}==

Using [https://metacpan.org/pod/Math::GMP Math::GMP]:

 

sub is_prime { Math::GMP->new(shift)->probab_prime(12); }

foreach my $p (2 .. 43_112_609) {

print "M$p\n" if is_prime($p) && is_mersenne_prime($p);

 

The ntheory module offers a couple options. This is direct:

{{libheader|ntheory}}

$|=1; # flush output on every print

my $n = 0;

print "M$n ";

}

However it uses knowledge from the thousands of CPU years spent by GIMPS to accelerate results for known values, so doesn't actually run the L-L test until after the 44th value, although code is included for C, Perl, and C+GMP. If we substitute <tt>Math::Prime::Util::GMP::is_mersenne_prime</tt> we can force the test to run.

 

A less opaque method uses the modular Lucas sequence, though it has no pretesting other than primality and calculates both <math>U_k</math> and <math>V_k</math> so won't be as fast:

use bigint try=>"GMP,Pari";

forprimes {

my $mp1 = 2**$p;

print "M$p\n" if $p == 2 || 0 == (lucas_sequence($mp1-1, 4, 1, $mp1))[0];

 

We can also use the core module <code>Math::BigInt</code>:

{{trans|Python}}

my $p = shift;

if ($p == 2) {

}

last if $count >= $upb_count;

 

=={{header|Phix}}==

{{libheader|Phix/mpfr}}

Native types work up to M31, after which inaccuracies mean that we need to wheel out gmp. Uses the mod replacement trick from C/FreeBASIC(gmp)

<span style="color: #008080;">with</span> <span style="color: #008080;">javascript_semantics</span>

<span style="color: #004080;">bool</span> <span style="color: #000000;">full</span> <span style="color: #0000FF;">=</span> <span style="color: #004600;">true</span> <span style="color: #000080;font-style:italic;">-- (see extended output below)</span>

<span style="color: #008080;">end</span> <span style="color: #008080;">while</span>

<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;">"completed in %s\n"</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>

{{out}}

<pre>

 

=={{header|PicoLisp}}==

(or

(= N 2)

(do (- P 2)

(setq S (% (- (* S S) 2) MP)) )

{{Out}}

<pre>: (for N 10000

be smaller than 1000.

 

lvars mp = 2**p - 1, sn = 4, i;

for i from 2 to p - 1 do

endif;

p + 2 -> p;

 

{{Out}} (obtained in few seconds)

M3

M5

M127

M521

 

=={{header|PowerShell}}==

This is just a translation of VBScript using [bigint], it could be optimized.

Flirt with the girl in the cubicle next door while it runs:

function Get-MersennePrime ([bigint]$Maximum = 4800)

{

}

}

Get-MersennePrime | Format-Wide {"{0,4}" -f $_} -Column 4 -Force

{{Out}}

<pre>

 

=={{header|Prolog}}==

show(Count) :-

findall(N, limit(Count, (between(2, infinite, N), mersenne_prime(N))), S),

N mod D =\= 0,

D2 is D + A, prime(N, D2, As).

{{Out}}

<pre>

=={{header|PureBasic}}==

PureBasic has no large integer support. Calculations are limited to the range of a signed quad integer type.

Protected mp.q = (1 << p) - 1, sn.q = 4, i

For i = 3 To p

Print(#CRLF$ + #CRLF$ + "Press ENTER to exit"): Input()

CloseConsole()

{{Out}}

<pre>M2

 

=={{header|Python}}==

from sys import stdout

from math import sqrt, log

if count >= upb_count: break

print

 

{{Out}}

 

=== Faster loop without division ===

def isqrt(n):

if n < 0:

if count >= upb_count: break

print

 

The main loop may be run much faster using [https://pypi.python.org/pypi/gmpy2 gmpy2] :

 

 

def lucas_lehmer(n):

s -= m

s -= two

 

With this, one can test all primes below 10^5 in around 24 hours on a Core i5 processor, with only one running thread.

 

Of course, they agree with [[oeis:A000043|OEIS A000043]].

 

=={{header|R}}==

# vectorized approach based on scalar code from primeSieve and mersenne in CRAN package `numbers`

require(gmp)

}

}

 

=={{header|Racket}}==

#lang racket

(require math)

 

(loop 3)

 

=={{header|Raku}}==

(formerly Perl 6)

multi is_mersenne_prime(Int $p) {

my $m_p = 2 ** $p - 1;

}

 

{{out|On my system}}

Letting it run for about a minute...

REXX won't have a problem with the large number of digits involved, but since it's an interpreted language,

<br>such massive number crunching isn't conducive in searching for large primes.

@.=0; @.2=1; @.3=1; @.5=1; @.7=1; @.11=1; @.13=1 /*a partial list of some low primes. */

!.=@.; !.0=1; !.2=1; !.4=1; !.5=1; !.6=1; !.8=1 /*#'s with these last digs aren't prime*/

return 'M'? /*return "modified" # (Mersenne index).*/

/*──────────────────────────────────────────────────────────────────────────────────────*/

{{out|output|text=&nbsp; when the following is used for input: &nbsp; &nbsp; <tt> 10000 </tt>}}

<pre>

 

=={{header|Ring}}==

see "Mersenne Primes :" + nl

for p = 2 to 18

next

return (sn=0)

 

=={{header|RPL}}==

%%HP: T(3)A(R)F(.); ; ASCII transfer header

\<< DUP LN DUP \pi * 4 SWAP / 1 + UNROT / * IP 2 { 2 } ROT 2 SWAP ; input n; n := Int(n/ln(n)*(1 + 4/(pi*ln(n)))), p:=2; (n ~ number of primes less then n, pi used here only as a convenience), 2 is assumed to be the 1st elemente in the list

NEXT NIP ; next i;

\>>

{{out}}

<pre>Outputs for arguments 130, 607 and 2281, respectively:

These take respectively 1m 22s on the real HP 50g, 4m 29s and 10h 29m 23s on the emulator (Debug4 running on PC under WinXP, Intel(R) Core(TM) Duo CPU T2350 @ 1.86GHz).

</pre>

 

=={{header|Ruby}}==

return true if p == 2

return false if p <= 1 || p.even?

break if count >= upb_count

end

 

{{out}}

 

=={{header|Rust}}==

 

extern crate rug;

}

}

with Intel(R) Core(TM) i7-5500U CPU @ 2.40GHz :

Less than 8,6 seconds to get the Mersenne primes up to 11213

{{libheader|Scala}}

In accordance with definition of Mersenne primes it could only be a Mersenne number with prime exponent.

import Stream._

 

}

println("That's All Folks!")

{{out}} After approx 20 minutes (2.10 GHz dual core)

<pre style="height: 30ex; overflow: scroll">Finding Mersenne primes in M[2..9999]

 

=={{header|Scheme}}==

(define (S n)

(let ((m (- (expt 2 n) 1)))

(display "M") (display p) (display " ")

(loop (- i 1) (next-prime p)))

 

M2 M3 M5 M7 M13...

 

=={{header|Scilab}}==

n=1

for iexp=2:iexpmax

end

if s==0 then printf("M%d ",iexp); end

{{out}}

<pre>M2 M3 M5 M7 M13 M17 M19</pre>

To get maximum speed the program should be [http://seed7.sourceforge.net/scrshots/comp.htm compiled] with -O2.

 

include "bigint.s7i";

end while;

writeln;

 

Original source: [http://seed7.sourceforge.net/algorith/math.htm#lucasLehmerTest lucasLehmerTest],

=={{header|Sidef}}==

{{trans|Raku}}

return true if (p == 2)

var s = 4

say "M#{n}"

}

{{out}}

<pre>

Uses a sieve of Eratosthenes.

 

 

}

 

 

 

{{out}}

2^19 - 1 = 524287 is prime

2^31 - 1 = 2147483647 is prime

 

=={{header|Tcl}}==

{{trans|Pop11}}

set n 0

set t [clock seconds]

}

 

{{Out}}

The program was still running, but as the next Mersenne prime is 19937

 

=={{header|TI-83 BASIC}}==

1→N

For(E,2,M)

Then:Disp E

End

{{out}}

<pre>2

=={{header|uBasic/4tH}}==

{{Trans|VBScript}}

n = 1

For j = 2 To m

Next i

If s = 0 Then Print "M";j

 

=={{header|VBScript}}==

n=1

out=""

End If

Next

{{Out}}

<pre>M2 M3 M5 M7 M13 </pre>

=={{header|Visual Basic .NET}}==

{{works with|Visual Basic .NET|2011}}

Private Sub btnGo_Click(sender As Object, e As EventArgs) Handles btnGo.Click

Const iexpmax = 31

Next iexp

End Sub

{{Out}}

<pre>

M2 M3 M5 M7 M13 M17 M19 M31 </pre>

 

=={{header|Wren}}==

{{libheader|Wren-big}}

{{libheader|wren-math}}

This follows the lines of my Kotlin entry but uses a table to quicken up the checking of the larger numbers. Despite this, it still takes just over 3 minutes to reach M4423. Surprisingly, using ''modPow'' rather than the simple ''%'' operator is even slower.

import "io" for Stdout

 

}

}

 

{{out}}

Took 181.271083 seconds.

</pre>

 

=={{header|Zig}}==

Zig supports 128 bit integer types natively, which means it's possible to find all Mersenne primes up to M127. (requires writing a modmul() routine to compute (a * b) % m for 128 bit integers without overflow.)

const std = @import("std");

const stdout = std.io.getStdOut().outStream();

return r;

}

{{Out}}

<pre>

=={{header|zkl}}==

Using [[Extensible prime generator#zkl]] and the GMP library.

primes:=Utils.Generator(Import("sieve").postponed_sieve);

fcn isMersennePrime(p){

s:=BN(4); do(p-2){ s.mul(s).sub(2).mod(mp) } // the % REALLY cuts down on mem usage

return(s==0);

Calculating S(n) is done in place (overwriting the value in the BigInt with the result); this really cuts down on memory usage.

println("Mersenne primes:");

This will "continuously" spew out Mersenne Primes.

 

 

Using five threads:

fcn(ps){ // a thread to generate primes, sleeps most of the time

Utils.Generator(Import("sieve").postponed_sieve).pump(ps)

}

println("Mersenne primes:");