88 lines
2.3 KiB
Forth
88 lines
2.3 KiB
Forth
\ 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 ;
|
|
|
|
\ 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
|
|
π-128 bounds do
|
|
i c@ dup * over > if leave then
|
|
dup i c@ mod 0= if 2drop false unloop exit then
|
|
loop drop ;
|
|
|
|
\ actual Miller-Rabin test
|
|
\
|
|
: factor-2s ( n -- s d )
|
|
0 swap
|
|
begin dup 1 and 0= while
|
|
swap 1+ swap 2/
|
|
repeat ;
|
|
|
|
: fermat-square-test ( n m s -- ? ) \ perform n = n^2 (mod m), s-1 times
|
|
1- 0 ?do
|
|
2dup - -1 =
|
|
if leave
|
|
then >r dup r@ mod* r>
|
|
loop
|
|
- -1 = ;
|
|
|
|
: strong-fermat-pseudoprime? ( n a -- ? )
|
|
over >r \ keep the modulus on the return stack
|
|
>r 1- factor-2s r> \ -- s d a
|
|
swap r@ mod^ \ s d a -- s, a^d (mod n)
|
|
dup 1 = \ a^d == 1 (mod n) => Fermat pseudoprime
|
|
if 2drop rdrop true
|
|
else r> rot fermat-square-test
|
|
then ;
|
|
|
|
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 cells bounds do
|
|
dup i @ strong-fermat-pseudoprime? invert
|
|
if drop false unloop exit then
|
|
cell +loop drop 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 ;
|