Miller–Rabin primality test: Difference between revisions
Content added Content deleted
(Added Wren) |
(Forth version of Miller Rabin test) |
||
Line 1,508: | Line 1,508: | ||
power(B, E - 1, B * Acc).</lang> |
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}}== |
=={{header|Fortran}}== |
||
===Direct translation=== |
===Direct translation=== |