Miller–Rabin primality test: Difference between revisions
m
→{{header|REXX}}
| (5 intermediate revisions by the same user not shown) | |||
Line 5,977:
These two issues prevent the testing of bigger primes.
Below version implements the full Miller-Rabin primality test included witnesses for deterministic test for x < 3317044064679887385961981. Improved functions '''Rand''' and '''Powermod''' (see elsewhere on RosettaCode) are added. Also some small functions '''Floor''' and '''Whole''' are included.
<syntaxhighlight lang="rexx">
parse version version
say version
glob. = ''
numeric digits 1000
say '25 numbers of the form 2^n-1, mostly Mersenne primes'
say 'Up to about 25 digits deterministic, above probabilistic'
say
mm = '2 3 5 7 11 13 17 19 23 31 61 89 97 107 113 127 131 521 607 1279 2203 2281 2293 3217 3221'
do nn = 1 to 25
a = Word(mm,nn); b = 2**a-1; l = Length(b)
call time('r'); p = Prime(b); e = time('e')
if l > 25 then
b = Left(b,12)'...'Right(b,13)
select
when p = 0 then
p = 'not'
when l < 26 then
p = 'for sure'
otherwise
p = 'probable'
end
say '2^('a'-1)' '=' b '('l' digits) is' p 'prime' '('format(e,,3) 'seconds)'
end
return
Prime:
/* Is a number prime? */
procedure expose glob.
arg x
/* Validity */
if \ Whole(x) then
return 'X'
if x < 2 then
return 'X'
/* Low primes also used as deterministic witnesses */
w1 = ' 2 3 5 7 11 13 17 19 23 29 31 37 41 '
/* Fast values */
w = x
if Pos(' 'w' ',w1) > 0 then
return 1
if x//2 = 0 then
return 0
if x//3 = 0 then
return 0
if Right(x,1) = 5 then
return 0
/* Miller-Rabin primality test */
numeric digits 2*Length(x)
d = x-1; e = d
/* Reduce n-1 by factors of 2 */
do s = -1 while d//2 = 0
d = d%2
end
/* Thresholds deterministic witnesses */
w2 = ' 2047 1373653 25326001 3215031751 2152302898747 3474749660383 341550071728321 0 ',
|| ' 3825123056546413051 0 0 318665857834031151167461 3317044064679887385961981 '
w = Words(w2)
/* Up to 13 deterministic trials */
if x < Word(w2,w) then do
do k = 1 to w
if x < Word(w2,k) then
leave
end
end
/* or 3 probabilistic trials */
else do
w1 = ' '
do k = 1 to 3
r = Rand(2,e)/1; w1 = w1||r||' '
end
k = k-1
end
/* Algorithm using generated witnesses */
do k = 1 to k
a = Word(w1,k); y = powermod(a,d,x)
if y = 1 then
iterate
if y = e then
iterate
do s
y = (y*y)//x
if y = 1 then
return 0
if y = e then
leave
end
if y <> e then
return 0
end
return 1
Floor:
/* Floor */
procedure expose glob.
arg x
/* Formula */
if Whole(x) then
return x
else
return Trunc(x)-(x<0)
Powermod:
/* Power modulus function x^y mod z */
procedure expose glob.
arg x,y,z
/* Validity */
if \ Whole(x) then
return 'X'
if \ Whole(x) then
return 'X'
if \ Whole(x) then
return 'X'
if x < 0 then
return 'X'
if y < 0 then
return 'X'
if z < 1 then
return 'X'
/* Special values */
if z = 1 then
return 0
/* Binary method */
numeric digits Max(Length(Format(x,,,0)),Length(Format(y,,,0)),2*Length(Format(z,,,0)))
b = x//z; r = 1
do while y > 0
if y//2 then
r = r*x//z
x = x*x//z; y = y%2
end
return r
Rand:
/* Random number 12 digits precision */
procedure expose glob.
arg x,y
/* Validity */
if x <> '' then do
if \ Whole(x) then
return 'X'
if \ Whole(x) then
return 'X'
if x >= y then
return 'X'
end
/* Fixed precision */
p = Digits(); numeric digits 30
/* Get and save first seed from system Date and Time */
if glob.rand = '' then do
a = Date('s'); b = Time('l')
glob.rand = Substr(b,10,3)||Substr(b,7,2)||Substr(b,4,2)||Substr(b,1,2)||Substr(a,7,2)||Substr(a,5,2)||Substr(a,3,2)
end
/* Calculate next random number method HP calculators */
glob.rand = Right((glob.rand*2851130928467),15)
/* Uniform deviate [0,1) */
z = '0.'Left(glob.rand,Length(glob.rand)-3)
numeric digits 12
if x = '' then
return z/1
else
return Floor(z*(y+1-x)+x)
Whole:
/* Is a number integer? */
procedure expose glob.
arg x
/* Formula */
return Datatype(x,'w')
</syntaxhighlight>
As in the previous version, k = 3 was hard coded and seems good enough for this test. A higher k gives more confidence on the probable primes, but also a longer running time.
This version will produce output (Windows 11, Intel i7, 4.5Ghz, 16G):
<pre>
REXX-ooRexx_5.0.0(MT)_64-bit 6.05 23 Dec 2022
25 numbers of the form 2^n-1, mostly Mersenne primes
Up to about 25 digits deterministic, above probabilistic
2^2-1 = 3 (1 digits) is for sure prime (0.000 seconds)
2^3-1 = 7 (1 digits) is for sure prime (0.000 seconds)
2^5-1 = 31 (2 digits) is for sure prime (0.000 seconds)
2^7-1 = 127 (3 digits) is for sure prime (0.000 seconds)
2^11-1 = 2047 (4 digits) is not prime (0.000 seconds)
2^13-1 = 8191 (4 digits) is for sure prime (0.000 seconds)
2^17-1 = 131071 (6 digits) is for sure prime (0.000 seconds)
2^19-1 = 524287 (6 digits) is for sure prime (0.000 seconds)
2^23-1 = 8388607 (7 digits) is not prime (0.000 seconds)
2^31-1 = 2147483647 (10 digits) is for sure prime (0.000 seconds)
2^61-1 = 2305843009213693951 (19 digits) is for sure prime (0.000 seconds)
2^89-1 = 618970019642...0137449562111 (27 digits) is probable prime (0.016 seconds)
2^97-1 = 158456325028...5187087900671 (30 digits) is not prime (0.000 seconds)
2^107-1 = 162259276829...1578010288127 (33 digits) is probable prime (0.000 seconds)
2^113-1 = 103845937170...0992658440191 (35 digits) is not prime (0.000 seconds)
2^127-1 = 170141183460...3715884105727 (39 digits) is probable prime (0.016 seconds)
2^131-1 = 272225893536...9454145691647 (40 digits) is not prime (0.000 seconds)
2^521-1 = 686479766013...8291115057151 (157 digits) is probable prime (0.453 seconds)
2^607-1 = 531137992816...3219031728127 (183 digits) is probable prime (0.765 seconds)
2^1279-1 = 104079321946...5703168729087 (386 digits) is probable prime (9.407 seconds)
2^2203-1 = 147597991521...7686697771007 (664 digits) is probable prime (36.527 seconds)
2^2281-1 = 446087557183...2418132836351 (687 digits) is probable prime (37.575 seconds)
2^2293-1 = 182717463422...4672097697791 (691 digits) is not prime (16.324 seconds)
2^3217-1 = 259117086013...7362909315071 (969 digits) is probable prime (111.373 seconds)
2^3221-1 = 414587337621...7806549041151 (970 digits) is not prime (36.390 seconds)
</pre>
Above 1000 digits it becomes very slow.
=={{header|Ring}}==
| |||