Miller–Rabin primality test: Difference between revisions
Forth version of Miller Rabin test
(Added Wren) |
(Forth version of Miller Rabin test) |
||
Line 1,508:
power(B, E - 1, B * Acc).</lang>
=={{header|Forth}}==
Forth only supports native ints (e.g. 64 bits on most modern machines), so this version uses a set of known deterministic bases that is known to be correct for 64 bit integers (and possibly greater). Prior to the Miller Rabin check, the "prime?" word checks for divisibility by some small primes.
<lang Forth>
\ small divisor check: true => possibly prime; false => definitely not prime.
\
31 constant π-128
create maybe-prime?
2 c, 3 c, 5 c, 7 c, 11 c, 13 c, 17 c, 19 c, 23 c, 29 c,
31 c, 37 c, 41 c, 43 c, 47 c, 53 c, 59 c, 61 c, 67 c, 71 c,
73 c, 79 c, 83 c, 89 c, 97 c, 101 c, 103 c, 107 c, 109 c, 113 c,
127 c,
does>
true -rot \ default: number not prime (not proven composite)
π-128 bounds do
i c@ dup * over > if leave then
dup i c@ mod 0= if 2drop false unloop exit then
loop drop ;
\ modular multiplication and exponentiation
\
: 3rd s" 2 pick" evaluate ; immediate
: mod* ( a b m -- a*b {mod m} )
>r um* r> ud/mod 2drop ;
: mod^ ( x n m -- x^n {mod m} )
>r 1 swap
begin ?dup while
dup 1 and 1 =
if
swap 3rd r@ mod* swap 1-
then dup 0>
if
rot dup r@ mod* -rot 2/
then
repeat nip rdrop ;
\ actual Miller-Rabin test
\
: factor-2s ( n -- s d )
0 swap
begin dup 1 and 0= while
swap 1+ swap 2/
repeat ;
4.759.123.141 drop constant mr-det-3 \ Deterministic threshold; 3 bases
create small-prime-bases 2 , 7 , 61 , \ deterministic up to mr-det-3
create large-prime-bases 2 , 325 , 9375 , 28178 , 450775 , 9780504 , 1795265022 , \ known to be deterministic for 64 bit integers.
: miler-rabin-bases ( n -- addr n )
mr-det-3 <
if small-prime-bases 3
else large-prime-bases 7
then ;
: miller-rabin-primality-test ( n -- f )
dup miler-rabin-bases
rot dup 1- factor-2s 0 locals| x d s n |
cells bounds DO
i @ d n mod^ dup to x 1 <>
if -1
begin 1+ dup s < x n 1- <> and while
x dup n mod* to x
repeat drop
x n 1- <>
if false unloop exit
then
then
cell +loop
true ;
: prime? ( n -- f )
dup 2 <
if drop false
else
dup maybe-prime?
if dup [ 127 dup * 1+ ] literal <
if drop true
else miller-rabin-primality-test
then
else drop false
then
then ;
</lang>
{{Out}}
Test on some Fermat numbers and some Mersenne numbers
<pre>
16 2^ 1+ dup . prime? . 65537 -1 ok
32 2^ 1+ dup . prime? . 4294967297 0 ok
53 2^ 1- dup . prime? . 9007199254740991 0 ok
61 2^ 1- dup . prime? . 2305843009213693951 -1 ok
</pre>
=={{header|Fortran}}==
===Direct translation===
| |||