48 lines
1.4 KiB
Plaintext
48 lines
1.4 KiB
Plaintext
/* Miller-Rabin algorithm is builtin, see function primep. Here is another implementation */
|
|
|
|
|
|
/* find highest power of p, p^s, that divide n, and return s and n / p^s */
|
|
|
|
facpow(n, p) := block(
|
|
[s: 0],
|
|
while mod(n, p) = 0 do (s: s + 1, n: quotient(n, p)),
|
|
[s, n]
|
|
)$
|
|
|
|
/* check whether n is a strong pseudoprime to base a; s and d are given by facpow(n - 1, 2) */
|
|
|
|
sppp(n, a, s, d) := block(
|
|
[x: power_mod(a, d, n), q: false],
|
|
if x = 1 or x = n - 1 then true else (
|
|
from 2 thru s do (
|
|
x: mod(x * x, n),
|
|
if x = 1 then return(q: false) elseif x = n - 1 then return(q: true)
|
|
),
|
|
q
|
|
)
|
|
)$
|
|
|
|
/* Miller-Rabin primality test. For n < 341550071728321, the test is deterministic;
|
|
for larger n, the number of bases tested is given by the option variable
|
|
primep_number_of_tests, which is used by Maxima in primep. The bound for deterministic
|
|
test is also the same as in primep. */
|
|
|
|
miller_rabin(n) := block(
|
|
[v: [2, 3, 5, 7, 11, 13, 17], s, d, q: true, a],
|
|
if n < 19 then member(n, v) else (
|
|
[s, d]: facpow(n - 1, 2),
|
|
if n < 341550071728321 then ( /* see http://oeis.org/A014233 */
|
|
for a in v do (
|
|
if not sppp(n, a, s, d) then return(q: false)
|
|
),
|
|
q
|
|
) else (
|
|
thru primep_number_of_tests do (
|
|
a: 2 + random(n - 3),
|
|
if not sppp(n, a, s, d) then return(q: false)
|
|
),
|
|
q
|
|
)
|
|
)
|
|
)$
|