Rosetta Code

Lucas-Lehmer test: Difference between revisions

Page Discussion

Lucas-Lehmer test (view source)

Revision as of 15:20, 5 May 2024

7,656 bytes added ,  14 days ago
→‎{{header|langur}}
m (→‎{{header|langur}}: change where() to filter())
 
(24 intermediate revisions by 13 users not shown)
Line 530:
</pre>
 
=={{header|BBC BASIC}}==
==={{header|BASIC256}}===
BASIC256 has no large integer support. Calculations are limited to the range of a integer type.
<syntaxhighlight lang="basic256">print "Mersenne Primes :"
for p = 2 to 18
if lucasLehmer(p) then print "M"; p
next p
end
 
function lucasLehmer (p)
mp = (2 ^ p) - 1
sn = 4
for i = 2 to p-1
sn = (sn ^ 2) - 2
sn = sn - (mp * floor(sn / mp))
next
return sn = 0
end function</syntaxhighlight>
 
==={{header|BBC BASIC}}===
{{works with|BBC BASIC for Windows}}
Using its native arithmetic BBC BASIC can only test up to M23.
Line 562 ⟶ 581:
M19
</pre>
 
==={{header|Craft Basic}}===
<syntaxhighlight lang="basic">let m = 7
let n = 1
 
for e = 2 to m
 
if e = 2 then
 
let s = 0
 
else
 
let s = 4
 
endif
 
let n = (n + 1) * 2 - 1
 
for i = 1 to e - 2
 
let s = (s * s - 2) % n
 
next i
 
if s = 0 then
print e, " ",
 
endif
 
next e</syntaxhighlight>
{{Out}}
<pre>2 3 5 7 </pre>
 
=={{header|BCPL}}==
Line 960 ⟶ 1,013:
static bool is_mersenne_prime(mpz_class p)
{
if( 2 == p ) {
return true;
else}
 
mpz_class s(4);
mpz_class div( (mpz_class(1) << p.get_ui()) - 1 );
for( mpz_class i(3); i <= p; ++i )
{
mpz_classs = (s * s - mpz_class(42)) % div ;
}
mpz_class div( (mpz_class(1) << p.get_ui()) - 1 );
 
for( mpz_class i(3); i <= p; ++i )
return ( s == mpz_class(0) {);
s = (s * s - mpz_class(2)) % div ;
}
 
return ( s == mpz_class(0) );
}
}
 
int main()
{
Line 1,152 ⟶ 1,204:
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|Delphi}}==
{{works with|Delphi|6.0}}
{{libheader|SysUtils,StdCtrls}}
 
 
<syntaxhighlight lang="Delphi">
 
function IsMersennePrime(P: int64): boolean;
{Test for Mersenne Prime - P cannot be bigger than 63}
{Because (1 shl 64) would be bigger than in64}
var S,MP: int64;
var I: integer;
begin
if P= 2 then Result:=true
else
begin
MP:=(1 shl P) - 1;
S:=4;
for I:=3 to P do
begin
S:=(S * S - 2) mod MP;
end;
Result:=S = 0;
end;
end;
 
 
procedure ShowMersennePrime(Memo: TMemo);
var Sieve: TPrimeSieve;
var I: integer;
begin
{Create Sieve}
Sieve:=TPrimeSieve.Create;
{Test cannot handle values bigger than 64}
Sieve.Intialize(64);
for I:=0 to Sieve.PrimeCount-1 do
if IsMersennePrime(Sieve.Primes[I]) then
begin
Memo.Lines.Add(IntToStr(Sieve.Primes[I]));
end;
Sieve.Free;
end;
 
 
 
</syntaxhighlight>
{{out}}
<pre>
2
3
5
7
13
17
19
31
 
Elapsed Time: 10.167 ms.
</pre>
 
 
=={{header|DWScript}}==
Line 1,183 ⟶ 1,296:
{{Out}}
<pre> M2 M3 M5 M7 M13 M17 M19 M31</pre>
 
=={{header|EasyLang}}==
{{trans|BASIC256}}
<syntaxhighlight>
write "Mersenne Primes: "
func lulehm p .
mp = bitshift 1 p - 1
sn = 4
for i = 2 to p - 1
sn = sn * sn - 2
sn = sn - (mp * (sn div mp))
.
return if sn = 0
.
for p = 2 to 23
if lulehm p = 1
write "M" & p & " "
.
.
</syntaxhighlight>
 
{{out}}
<pre>
Mersenne Primes: M3 M5 M7 M13 M17 M19
</pre>
 
=={{header|EchoLisp}}==
Line 1,338 ⟶ 1,476:
1 15 lucas-lehmer</syntaxhighlight>
==Alternate version to handle 64 and 128 bit integers.==
Forth supports modular multiplication without overflow by having mixed mode operations that work on 128 bit intermediate results. It's also possible to do do the Lucas-Lehmer test using double-precision (128 bit) integers, though support for that is more limited in the Forth language, so it requires writing more support code. This version contains the code for both 64 bit (mixed mode) and 128 bit (double precision mode)
<syntaxhighlight lang="forth">
18 constant π-64 \ count of primes < 64
Line 1,442 ⟶ 1,580:
END PROGRAM LUCAS_LEHMER</syntaxhighlight>
===128 Bit Version===
This version can find all Mersenne Primes up to M127. Its based on the Zig code but written in Fortran 77 style (fixed format, unstructured loops.) Works with GNU Fortran which has 128 bit integer support.
 
<syntaxhighlight lang="fortran">
PROGRAM Mersenne Primes
IMPLICIT INTEGER (a-z)
LOGICAL is mersenne prime
 
PARAMETER (sz primes = 31)
INTEGER*1 primes(sz primes)
 
DATA primes
& /2, 3, 5, 7, 11, 13, 17, 19, 23, 29,
& 31, 37, 41, 43, 47, 53, 59, 61, 67, 71,
& 73, 79, 83, 89, 97, 101, 103, 107, 109, 113,
& 127/
 
PRINT *, 'These Mersenne numbers are prime:'
 
DO 10 i = 1, sz primes
p = primes(i)
10 IF (is mersenne prime(p))
& WRITE (*, '(I5)', ADVANCE = 'NO'), p
 
PRINT *
END
 
 
FUNCTION is mersenne prime(p)
IMPLICIT NONE
LOGICAL is mersenne prime
INTEGER*4 p, i
INTEGER*16 n, s, modmul
 
IF (p .LT. 3) THEN
is mersenne prime = p .EQ. 2
ELSE
n = 2_16**p - 1
s = 4
DO 10 i = 1, p - 2
s = modmul(s, s, n) - 2
10 IF (s .LT. 0) s = s + n
is mersenne prime = s .EQ. 0
END IF
END
 
 
FUNCTION modmul(a0, b0, m)
IMPLICIT INTEGER*16 (a-z)
 
modmul = 0
a = MODULO(a0, m)
b = MODULO(b0, m)
 
10 IF (b .EQ. 0) RETURN
IF (MOD(b, 2) .EQ. 1) THEN
IF (a .LT. m - modmul) THEN
modmul = modmul + a
ELSE
modmul = modmul - m + a
END IF
END IF
b = b / 2
IF (a .LT. m - a) THEN
a = a * 2
ELSE
a = a - m + a
END IF
GO TO 10
END
</syntaxhighlight>
{{Out}}
<pre>
These Mersenne numbers are prime:
2 3 5 7 13 17 19 31 61 89 107 127
</pre>
 
=={{header|FreeBASIC}}==
Line 2,067 ⟶ 2,281:
=={{header|langur}}==
{{trans|D}}
It is theoretically possible to test to the 47th Mersenne prime, as stated in the task description, but it could take a while. As for the limit, it would be extremely high.
{{works with|langur|0.11}}
With slight syntactic differences, this would also work with earlier versions of langur if you limit it to smaller numbers. 0.8 uses arbitrary precision.
 
<syntaxhighlight lang="langur">val .isPrime = f fn(.i) == 2 or{
.i >== 2 and not any f(.x)or .i div .x, pseries> 2 to .i ^/ 2and
not any(fn .x: .i div .x, pseries 2 .. .i ^/ 2)
}
 
val .isMersennePrime = ffn(.p) {
if .p == 2: return true
if not .isPrime(.p): return false
 
val .mp = 2 ^ .p - 1
for[.s=4] of 3 to.. .p {
.s = (.s ^ 2 - 2) rem .mp
} == 0
}
 
writeln join " ", map ffn .x: $"M\{{.x;}}", filter .isMersennePrime, series 2300</syntaxhighlight>
</syntaxhighlight>
 
{{out}}
<pre>M2 M3 M5 M7 M13 M17 M19 M31 M61 M89 M107 M127 M521 M607 M1279 M2203 M2281
</pre>
 
=={{header|Mathematica}}/{{header|Wolfram Language}}==
Line 2,313 ⟶ 2,533:
 
ll.gp
<syntaxhighlight lang="parigpc">/* ll(p): input odd prime 'p'. */
/* ll(p): input odd prime 'p'. */
/* returns '1' if 2^p-1 is a Mersenne prime. */
ll(p) = {
/* trial division up to a reasonable depth (time ratio tdiv/llt approx. 0.2) */
{
my(l=log(p), ld=log(l));
ll(p) =
forprimestep(q = 1, sqr(ld)^(l/log(2))\4, p+p,
my(s = 4, m = 2^p-1, n = m+2, d = 2*p, q = 1, r = exponent(p)+1);
if(Mod(2,q)^p == 1, return)
);
/* do some trial division up to a reasonable depth. */
/* Lucas-Lehmer test with fast modular reduction. */
for(i = 1, (r/4.)^r\p,
qmy(s=4, +m=2^p-1, dn=m+2);
iffor(Mod(2,q)^pi == 13, return)p,
s = sqr(s);
s = bitand(s,m)+ s>>p;
if(s >= m, s -= n, s -= 2)
);
!s
}; /* end ll */
 
/* Lucas-Lehmer test with fast modular reduction. */
for(i = 3, p,
s = sqr(s);
s = bitand(s,m)+s >> p;
if(s >= m, s -= n, s -= 2)
);
!s
}; /* end ll */</syntaxhighlight>
 
/* get Mersenne primes in range [a,b] */
llrun(a, b) = {
my(t=0, c=0, p=2, thr=default(nbthreads));
if(a <= 2,
c++;
printf("#%d\tM%d\t%3dh, %2dmin, %2d,%03d ms\n", c, p, t\3600000, t\60000%60, t\1000%60, t%1000);
a = 3;
);
gettime();
parforprime(p= a, b, ll(p), d, /* ll(p) -> d copy from parallel world into real world. */
if(d,
t += gettime()\thr;
c++;
printf("#%d\tM%d\t%3dh, %2dmin, %2d,%03d ms\n", c, p, t\3600000, t\60000%60, t\1000%60, t%1000)
)
)
}; /* end llrun */
 
<syntaxhighlight lang="parigp">
 
\\ export(ll); /* export 'll' to the parallel world/* if parallelrunning processingll isas enabledscript */
</syntaxhighlight>
Compiled with gp2c option: gp2c-run -g ll.gp.
 
<syntaxhighlight lang="c">llrun(2, 132049)</syntaxhighlight>
 
/* get the first 30 Mersenne primes */
{
t=0;c=1;p=2;gettime();
thr=default(nbthreads);
printf("#%d\tM%d\t%3dh, %2dmin, %2d,%03d ms\n", c, p, t\3600000, t\60000%60, t\1000%60, t%1000);
parforprime(p=3, 132049, ll(p), d, /* ll(p) -> d copy from parallel world into real world. */
if(d,
c++;
t += gettime()\thr;
printf("#%d\tM%d\t%3dh, %2dmin, %2d,%03d ms\n", c, p, t\3600000, t\60000%60, t\1000%60, t%1000)
)
)
};</syntaxhighlight>
{{out}}
Compiled with gp2c option: gp2c-run -g ll.gp.
 
Done on Intel(R) Core(TM) i5-8250U CPU @ 1.60GHz, 4 hyperthreaded cores.
Line 2,370 ⟶ 2,592:
#11 M107 0h, 0min, 0,000 ms
#12 M127 0h, 0min, 0,000 ms
#13 M607M521 0h, 0min, 0,002001 ms
#14 M521M607 0h, 0min, 0,002001 ms
#15 M1279 0h, 0min, 0,011007 ms
#16 M2203 0h, 0min, 0,040030 ms
#17 M2281 0h, 0min, 0,042033 ms
#18 M3217 0h, 0min, 0,096079 ms
#19 M4253 0h, 0min, 0,191163 ms
#20 M4423 0h, 0min, 0,213186 ms
#21 M9689 0h, 0min, 1,951789 ms
#22 M9941 0h, 0min, 2,187022 ms
#23 M11213 0h, 0min, 32,365835 ms
#24 M19937 0h, 0min, 2923,758858 ms
#25 M21701 0h, 0min, 4335,233268 ms
#26 M23209 0h, 0min, 5845,444233 ms
#27 M44497 0h, 9min6min, 4153,164051 ms
#28 M86243 1h, 17min 3min, 5841,437811 ms
#29 M110503 3h2h, 1min29min, 014,935055 ms
#30 M132049 5h4h, 38min42min, 4527,848694 ms
? ##
*** last result: cpu time 44h37h, 51min31min, 2841,870619 ms, real time 5h4h, 38min42min, 246,249515 ms.</pre>
?</pre>
 
=={{header|Pascal}}==
Line 2,979 ⟶ 3,200:
 
Of course, they agree with [[oeis:A000043|OEIS A000043]].
 
=={{header|Quackery}}==
 
<code>eratosthenes</code> and <code>isprime</code> are defined at [[Sieve of Eratosthenes#Quackery]].
 
<syntaxhighlight lang="Quackery"> [ dup temp put
dup bit 1 -
4
rot 2 - times
[ dup *
dup temp share >>
dip [ over & ] +
2dup > not if
[ over - ]
2 - ]
0 =
nip temp release ] is l-l ( n --> b )
 
25000 eratosthenes
[] 25000 times [ i^ isprime if [ i^ join ] ]
1 split
witheach
[ dup l-l iff join else drop ]
echo</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 ]</pre>
 
=={{header|R}}==
Line 3,186 ⟶ 3,435:
 
=={{header|RPL}}==
===RPL HP-50 series===
<syntaxhighlight lang="rpl">
%%HP: T(3)A(R)F(.); ; ASCII transfer header
Line 3,204 ⟶ 3,454:
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>
===RPL HP-28 series===
Unlike RPL implemented on HP-50 series, the first version of the language does not feature big integers, modular arithmetic operators, prime number test functions, nor even modulo operator for unsigned integers.
Let's build them all...
{{works with|Halcyon Calc|4.2.7}}
{| class="wikitable"
! RPL code
! Comment
|-
|
≪
'''IF''' { #1 #2 #3 #5 } OVER POS '''THEN''' #1 ≠
'''ELSE IF''' # 1d DUP2 AND ≠ OVER 3 DUP2 / * == OR
'''THEN''' DROP 0
'''ELSE''' DUP B→R √ → divm
≪ 1 SF 4
5 divm '''FOR''' n
'''IF''' OVER n DUP2 / * ==
'''THEN''' 1 CF divm 'n' STO '''END'''
6 SWAP - DUP '''STEP'''
DROP2 1 FS?
≫ '''END END'''
≫ '<span style="color:blue">bPRIM?</span>' STO
≪ → m
≪ #1
'''WHILE''' OVER #0 > '''REPEAT'''
IF OVER #1 AND #1 == '''THEN'''
3 PICK * m / LAST ROT * - '''END'''
SWAP SR SWAP
ROT DUP * m / LAST ROT * - ROT ROT
'''END'''
ROT ROT DROP2
≫ ≫ '<span style="color:blue">MODXP</span>' STO
≪ 2 OVER ^ R→B 1 - → mp
≪ #4
3 ROT '''FOR''' n
#2 mp <span style="color:blue">MODXP</span>
'''IF''' DUP #2 < '''THEN''' mp + '''END''' #2 -
'''NEXT'''
#0 ==
≫ ≫ '<span style="color:blue">MSNP?</span>' STO
≪ { 2 } 3 32 '''FOR''' j
'''IF''' j R→B <span style="color:blue">bPRIM?</span>
'''THEN IF''' j <span style="color:blue">MNSP?</span> '''THEN''' j + '''END END'''
'''NEXT'''
≫ '<span style="color:blue">TASK</span>' STO
|
<span style="color:blue">bPRIM?</span> ''( #a → boolean )''
return 1 if a is 2, 3 or 5 and 0 if a is 1
if 2 or 3 divides a
return 0
else store sqrt(a)
d = 4 ; flag 1 set while presumed prime
for n=5 to sqrt(a)
if d divides a
prepare loop exit
d = 6-d ; n += d
clean stack, return result
<span style="color:blue">MODXP</span> ''( #base #exp #m → #mod(base^exp,m) )''
result = 1;
while (exp > 0) {
if ((exp & 1) > 0)
result = (result * base) % m;
exp >>= 1;
base = (base * base) % m;
}
clean stack, return result
<span style="color:blue">MNSP?</span> ''( p → boolean )''
s0 = 4
loop p-2 times
r = mod(s(n-1)^2 mod Mp
s(n) = (r - 2) mod Mp
return 1 if s(p-2)=0 mod Mp
|}
{{out}}
<pre>
1: { 2 3 5 7 13 17 19 31 }
</pre>
Runs in 48 seconds on a standard HP-28S.
 
=={{header|Ruby}}==
Line 3,801 ⟶ 4,145:
{{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.
<syntaxhighlight lang="ecmascriptwren">import "./big" for BigInt
import "./math" for Int
import "io" for Stdout
 
Line 3,848 ⟶ 4,192:
{{libheader|Wren-gmp}}
Same approach as the CLI version but now uses GMP. Far quicker, of course, as we can now check up to M110503 in less time than before.
<syntaxhighlight lang="ecmascriptwren">import "./gmp" for Mpz
import "./math" for Int
 
Line 3,892 ⟶ 4,236:
Took 127.317323 seconds.
</pre>
 
=={{header|Yabasic}}==
<syntaxhighlight lang="yabasic">print "Mersenne Primes :"
for p = 2 to 20
if lucasLehmer(p) print "M",p
next p
end
 
sub lucasLehmer (p)
mp = (2 ^ p) - 1
sn = 4
for i = 2 to p-1
sn = (sn ^ 2) - 2
sn = sn - (mp * floor(sn / mp))
next
return sn = 0
end sub</syntaxhighlight>
 
=={{header|Zig}}==
Langurmonkey
885

edits

Retrieved from "https://rosettacode.org/wiki/Special:MobileDiff/332634...364251"