Rosetta Code

Lucas-Lehmer test: Difference between revisions

Page Discussion

Lucas-Lehmer test (view source)

Revision as of 11:36, 21 February 2023

142 bytes added ,  1 year ago
PARI/GP →‎Version with trial division and fast modular reduction: improved trial division
m (Fix minor grammar errors.)
Line 2,337:
 
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) = { ▼
/* do some 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 s>> p;▼
if(s >= m, s -= n, s -= 2)▼
);▼
!s▼
}; /* end ll */</syntaxhighlight>▼
 
▲ /* 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=1, p=2, thr=default(nbthreads));▼
▲ for if(ia <= 3, p2,
printf("#%d\tM%d\t%3dh, %2dmin, %2d,%03d ms\n", c, p, t\3600000, t\60000%60, t\1000%60, t%1000);▼
c++ a = 3;▼
);▼
▲ s = sqrgettime(s);
parforprime(p=3 a, 132049b, 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 lang="parigp">
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,394 ⟶ 2,395:
#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}}==
Nebulus
40

edits

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