Miller–Rabin primality test: Difference between revisions
m
→{{header|REXX}}
VincentARM (talk | contribs) (add task to aarch64 assembly raspberry pi 3) |
|||
| (9 intermediate revisions by 3 users not shown) | |||
Line 5,680:
False
>>> </pre>
=={{header|Quackery}}==
Translated from the current (22 Sept 2023) pseudocode from [[wp:Miller-Rabin primality test#Miller–Rabin_test|Wikipedia]], which is:
'''Input #1''': ''n'' > 2, an odd integer to be tested for primality
'''Input #2''': ''k'', the number of rounds of testing to perform
'''Output''': “''composite''” if ''n'' is found to be composite, “''probably prime''” otherwise
'''let''' ''s'' > 0 and ''d'' odd > 0 such that ''n'' − 1 = 2<sup>''s''</sup>''d'' <u># by factoring out powers of 2 from ''n'' − 1</u>
'''repeat''' ''k'' '''times''':
''a'' ← random(2, ''n'' − 2) # ''n'' is always a probable prime to base 1 and ''n'' − 1
''x'' ← ''a''<sup>''d''</sup> mod ''n''
'''repeat''' ''s'' '''times''':
''y'' ← ''x''<sup>2</sup> mod ''n''
'''if''' ''y'' = 1 and ''x'' ≠ 1 and ''x'' ≠ ''n'' − 1 '''then''' <u># nontrivial square root of 1 modulo ''n''</u>
'''return''' “''composite''”
''x'' ← ''y''
'''if''' ''y'' ≠ 1 '''then'''
'''return''' “''composite''”
'''return''' “''probably prime''”
As Quackery does not have variables, I have included a "builder" (i.e. a compiler extension) to provide basic global variables, in order to facilitate comparison to the pseudocode.
<code>var myvar</code> will create a variable called <code>myvar</code> initialised with the value <code>0</code>, and additionally the words <code>->myvar</code> and <code>myvar-></code>, which assign a value to <code>myvar</code> and return the value of <code>myvar</code> respectively.
The variable <code>c</code> is set to <code>true</code> if the number being tested is composite. This is included as Quackery does not have a <code>break</code> that can exit nested iterative loops.
<code>eratosthenes</code> and <code>isprime</code> are defined at [[Sieve of Eratosthenes#Quackery]].
<code>**mod</code> is defined at [[Modular exponentiation#Quackery]].
<syntaxhighlight lang="Quackery"> [ nextword
dup $ "" = if
[ $ "var needs a name"
message put bail ]
$ " [ stack 0 ] is "
over join
$ " [ "
join over join
$ " replace ] is ->"
join over join
$ " [ "
join over join
$ " share ] is "
join swap join
$ "-> " join
swap join ] builds var ( [ $ --> [ $ )
10000 eratosthenes
[ 1 & not ] is even ( n --> b )
var n var d var s var a
var t var x var y var c
[ dup 10000 < iff isprime done
dup even iff
[ drop false ] done
dup ->n
[ 1 - 0 swap
[ dup even while
1 >> dip 1+ again ] ]
->d ->s
false ->c
20 times
[ n-> 2 - random 2 + ->a
s-> ->t
a-> d-> n-> **mod ->x
s-> times
[ x-> 2 n-> **mod ->y
y-> 1 =
x-> 1 !=
x-> n-> 1 - !=
and and iff
[ true ->c conclude ]
done
y-> ->x ]
c-> iff conclude done
y-> 1 != if
[ true ->c conclude ] ]
c-> not ] is prime ( n --> b )
say "Primes < 100:" cr
100 times [ i^ prime if [ i^ echo sp ] ]
cr cr
[] 1000 times
[ i^ 9223372036854774808 + dup
prime iff join else drop ]
say "There are "
dup size echo
say " primes between 9223372036854774808 and 9223372036854775808:"
cr cr
witheach [ echo i^ 1+ 4 mod iff sp else cr ]
cr cr
4547337172376300111955330758342147474062293202868155909489 dup echo
prime iff
[ say " is prime." ]
else
[ say " is not prime." ]
cr cr
4547337172376300111955330758342147474062293202868155909393 dup echo
prime iff
[ say "is prime." ]
else
[ say " is not prime." ]</syntaxhighlight>
{{out}}
<pre>Primes < 100:
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
There are 23 primes between 9223372036854774808 and 9223372036854775808:
9223372036854774893 9223372036854774917 9223372036854774937 9223372036854774959
9223372036854775057 9223372036854775073 9223372036854775097 9223372036854775139
9223372036854775159 9223372036854775181 9223372036854775259 9223372036854775279
9223372036854775291 9223372036854775337 9223372036854775351 9223372036854775399
9223372036854775417 9223372036854775421 9223372036854775433 9223372036854775507
9223372036854775549 9223372036854775643 9223372036854775783
4547337172376300111955330758342147474062293202868155909489 is prime.
4547337172376300111955330758342147474062293202868155909393 is not prime.</pre>
=={{header|Racket}}==
Line 5,845 ⟶ 5,969:
for 1──►10000, K=10, Miller─Rabin primality test found 1229 primes {out of 1229}.
</pre>
Above program has 2 issues:
1) The REXX BIF '''random()''' is limited to the range 1 to 99999. Thus the statement '''?= random(2, nL)''' runs into trouble very soon.
2) In the statement '''x= ?**d // n''' (meaning ?^d mod n) the part '''?**d''' overflows already on modest values of ? and d. REXX does not have 'modular exponentation'.
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}}==
Line 6,710 ⟶ 7,052:
I've therefore used this method to check the same numbers as in my Kotlin entry.
<syntaxhighlight lang="
var iters = 10
| |||