

 This page uses content from Wikipedia. The original article was at Miller–Rabin primality test. The list of authors can be seen in the page history. As with Rosetta Code, the text of Wikipedia is available under the GNU FDL. (See links for details on variance)



The Miller–Rabin primality test or Rabin–Miller primality test is a primality test: an algorithm which determines whether a given number is prime or not.

The algorithm, as modified by Michael O. Rabin to avoid the generalized Riemann hypothesis, is a probabilistic algorithm.

The pseudocode, from Wikipedia is:

Input: n > 2, an odd integer to be tested for primality;
       k, a parameter that determines the accuracy of the test
Output: composite if n is composite, otherwise probably prime
write n − 1 as 2s·d with d odd by factoring powers of 2 from n − 1
LOOP: repeat k times:
   pick a randomly in the range [2, n − 1]
   x ← ad mod n
   if x = 1 or x = n − 1 then do next LOOP
   repeat s − 1 times:
      x ← x2 mod n
      if x = 1 then return composite
      if x = n − 1 then do next LOOP
   return composite
return probably prime

 The nature of the test involves big numbers, so the use of "big numbers" libraries (or similar features of the language of your choice) are suggested, but not mandatory.
 Deterministic variants of the test exist and can be implemented as extra (not mandatory to complete the task)


(ns test-p.core
  (:require [clojure.math.numeric-tower :as math])
  (:require [clojure.set :as set]))
 
(def WITNESSLOOP "witness")
(def COMPOSITE "composite")
 
(defn m* [p q m]
  " Computes (p*q) mod m "
  (mod (*' p q) m))
 
(defn power
  "modular exponentiation (i.e. b^e mod m"
  [b e m]
  (loop [b b, e e, x 1]
    (if (zero? e)
      x
      (if (even? e) (recur (m* b b m) (quot e 2) x)
                    (recur (m* b b m) (quot e 2) (m* b x m))))))
 
; Sequence of random numbers to use in the test
(defn rand-num [n]
  " random number between 2 and n-2 "
  (bigint (math/floor (+' 2 (*' (- n 4) (rand))))))
 
; Unique set of random numbers
(defn unique-random-numbers [n k]
  " k unique random numbers between 2 and n-2 "
  (loop [a-set #{}]
    (cond
      (>= (count a-set) k) a-set
      :else (recur (conj a-set (rand-num n))))))
 
(defn find-d-s [n]
  " write n − 1 as 2s·d with d odd "
  (loop [d (dec n), s 0]
    (if (odd? d)
      [d s]
      (recur (quot d 2) (inc s)))))
 
(defn random-test
  ([n] (random-test n (min 1000 (bigint (/ n 2)))))
  ([n k]
  " Random version of primality test"
  (let [[d s] (find-d-s n)
        ; Individual Primality Test
        single-test (fn [a s]
                      (let [z (power a d n)]
                       (if (some #{z} [1 (dec n)])
                         WITNESSLOOP
                         (loop [x (power z 2 n), r s]
                           (cond
                             (= x 1) COMPOSITE
                             (= x (dec n)) WITNESSLOOP
                             (= r 0) COMPOSITE
                             :else (recur (power x 2 n) (dec r)))))))]
    ; Apply Test
    ;(not-any? #(= COMPOSITE (local-test % s))
    ;          (unique-random-numbers n k))))
    (not-any? #(= COMPOSITE (single-test % s)) (unique-random-numbers n k)))))
 
;; Testing
(println "Primes beteen 900-1000:")
(doseq [q (range 900 1000)
        :when (random-test q)]
  (print " " q))
(println)
(println "Is Prime?" 4547337172376300111955330758342147474062293202868155909489 (random-test 4547337172376300111955330758342147474062293202868155909489))
(println "Is Prime?" 4547337172376300111955330758342147474062293202868155909393 (random-test 4547337172376300111955330758342147474062293202868155909393))
(println "Is Prime?" 643808006803554439230129854961492699151386107534013432918073439524138264842370630061369715394739134090922937332590384720397133335969549256322620979036686633213903952966175107096769180017646161851573147596390153
         (random-test 643808006803554439230129854961492699151386107534013432918073439524138264842370630061369715394739134090922937332590384720397133335969549256322620979036686633213903952966175107096769180017646161851573147596390153))
 
(println "Is Prime?" 743808006803554439230129854961492699151386107534013432918073439524138264842370630061369715394739134090922937332590384720397133335969549256322620979036686633213903952966175107096769180017646161851573147596390153
         (random-test 743808006803554439230129854961492699151386107534013432918073439524138264842370630061369715394739134090922937332590384720397133335969549256322620979036686633213903952966175107096769180017646161851573147596390153))
 