Miller-Rabin primality test (Clojure)

From LiteratePrograms

Other implementations: C | C, GMP | Clojure | Groovy | Java | Python | Ruby | Scala

The Miller-Rabin primality test is a simple probabilistic algorithm for determining whether a number is prime or composite that is easy to implement. It proves compositeness of a number using the following formulas:

Suppose 0 < a < n is coprime to n (this is easy to test using the GCD). Write the number n−1 as , where d is odd. Then, provided that all of the following formulas hold, n is composite:

for all

If a is chosen uniformly at random and n is prime, these formulas hold with probability 1/4. Thus, repeating the test for k random choices of a gives a probability of 1 − 1 / 4^k that the number is prime. Moreover, Gerhard Jaeschke showed that any 32-bit number can be deterministically tested for primality by trying only a=2, 7, and 61.

(defn factorize-out [n x]
  (loop [acc n exp 0]
    (if (zero? (rem acc x))
      (recur (/ acc x) (inc exp))
      (hash-map :exponent exp :rest acc))))
(use '[clojure.contrib.math :only (expt)])
(defn expt-rem [n e m]
  (loop [r 1, b n, e e]
                (if (zero? e) r
                        (recur
                                (if (odd? e) (rem (* r b) m) r)
                                (rem (expt b 2) m)
                                (bit-shift-right e 1)))))
(defn miller-rabin? [accuracy num]
        (cond
                (even? num) (= num 2)
                (< num 2) 'false?
                :else
                (let
                        [m (factorize-out (dec num) 2)
                        d (:rest m)
                        s (:exponent m)
                        witness? (fn [r x]
                                (cond
                                        (or (= x 1)(> r (dec s))) 'false
                                        (= x (dec num)) 'true
                                        :else (recur (inc r)(rem (* x x) num))))
                        investigate? (fn [k]
                                (if (> k accuracy) 'true
                                        (let
                                                [a (+ 2 (rand-int (- num 4)))
                                                x (expt-rem a d num)]
                                                (if (or (= x 1)(= x (dec num))(witness? 1 (expt-rem x 2 num)))
                                                        (recur (inc k))
                                                        'false))))]
                        (investigate? 1))))
(def *prime-accuracy* 10)
(def isPrime? (partial miller-rabin? *prime-accuracy*))

Download code