Miller–Rabin primality test: Difference between revisions
Refactored the code in a way that is more idiomatic to Forth.
(Refactored the code in a way that is more idiomatic to Forth.) |
|||
Line 1,510:
=={{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▼
π-128 bounds do▼
i c@ dup * over > if leave then▼
dup i c@ mod 0= if 2drop false unloop exit then▼
loop drop ;▼
▲<lang Forth>
\ modular multiplication and exponentiation
\
Line 1,544 ⟶ 1,530:
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
Line 1,554 ⟶ 1,555:
4.759.123.141 drop constant mr-det-3 \ Deterministic threshold; 3 bases
: fermat-square-test ( n m s -- ? ) \ perform n = n^2 (mod m), s-1 times
1- 0 ?do
: strong-fermat-pseudoprime? ( n a -- ? )
over >r \ keep the modulus on the return stack
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 ;
create small-prime-bases 2 , 7 , 61 , \ deterministic up to mr-det-3
Line 1,565 ⟶ 1,583:
: miller-rabin-primality-test ( n -- f )
dup miler-rabin-bases cells bounds do
dup i @ strong-fermat-pseudoprime? invert
▲ rot dup 1- factor-2s 0 locals| x d s n |
cell +loop drop true ;
▲ if -1
▲ x dup n mod* to x
▲ then
▲ cell +loop
▲ true ;
: prime? ( n -- f )
| |||