+- HP Forums (https://www.hpmuseum.org/forum)
+-- Forum: HP Software Libraries (/forum-10.html)
+--- Forum: General Software Library (/forum-13.html)
+--- Thread: (48G) Miller-Rabin Primality Test (/thread-11026.html)
(48G) Miller-Rabin Primality Test - Thomas Klemm - 07-09-2018 03:12 PM
Code:
@PRIME? (n -- 1 if n is prime else 0 )
« DUP 1 -
@TIMES ( a b n -- a * b % n )
« 1000000 → a b n k
« a k MOD a OVER -
b k MOD b OVER -
→ y x v u
« x u * n NEG MOD
x v * n MOD
+ n NEG MOD
y u * n MOD
+ n NEG MOD
y v *
+ n MOD
»
»
»
→ n m TIMES
« 0 m
DO
SWAP 1 +
SWAP 2 /
UNTIL
DUP 2 MOD
END → s d
« { 2 13 23 1662803 }
IF DUP n POS @ one of the primes in the list
THEN DROP 1
ELSE
IFERR
WHILE DUP SIZE
REPEAT
@POWER ( a -- a ^ d (mod n) )
1 OVER HEAD d
WHILE DUP
REPEAT → r a b
« IF b 2 MOD
THEN a r n TIMES EVAL @ r → a * r (mod n)
ELSE r
END
a DUP n TIMES EVAL @ a → a ^ 2 (mod n)
b 2 / IP
»
END
DROP2
@ x = a ^ d (mod n)
IF DUP 1 ==
THEN DROP 0 @ x == 1 → prime
ELSE
IFERR s
WHILE DUP
REPEAT
SWAP
IF DUP m ==
THEN 0 DOERR
END
DUP n TIMES EVAL @ x → x ^ 2 (mod n)
SWAP 1 -
END
THEN DROP2 0 @ x == n - 1 → prime
ELSE DROP2 1 @ else → composite
END
END
IF @ composite?
THEN 0 DOERR
END
TAIL
END
THEN DROP 0
ELSE DROP 1
END
END
»
»
»Benchmark using HP-48GX
999,999,999,999 → 0 in ~ 8 seconds.
999,999,999,989 → 1 in ~ 31 seconds.
177,777,773 → 1 in ~ 21 seconds.