/* 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
      )
   )
)$