Initial data commit

Task/Miller-Rabin-primality-test/Forth/miller-rabin-primality-test.fth

@@ -0,0 +1,88 @@
+\ 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 ;