; modest prime-number library
(defn gcd "greatest common divisor" [a b]
(if (zero? b) a (gcd b (mod a b))))
(println "gcd of 8 and 12 is" (gcd 8 12))
(defn powerMod "modular exponentiation" [b e m]
(defn m* [p q] (mod (* p q) m))
(loop [b b, e e, x 1]
(if (zero? e) x
(if (even? e) (recur (m* b b) (/ e 2) x)
(recur (m* b b) (quot e 2) (m* b x))))))
(def b 2988348162058574136915891421498819466320163312926952423791023078876139)
(def e 2351399303373464486466122544523690094744975233415544072992656881240319)
(def m 10000000000000000000000000000000000000000)
(println "b ^ e (mod m) =" (powerMod b e m))
(defn primes "primes less than n" [n]
(let [sieve (boolean-array n true)]
(loop [p 2, ps (list)]
(cond (= n p) (reverse ps)
(aget sieve p)
(do (doseq [i (range (* p p) n p)]
(aset sieve i false))
(recur (+ p 1) (cons p ps)))
:else (recur (+ p 1) ps)))))
(println "primes less than 30 are" (primes 30))
(defn prime? "miller-rabin primality test" [n]
(defn witness? [n a]
(loop [d (- n 1), s 0]
(if (even? d) (recur (/ d 2) (+ s 1))
(let [t (powerMod a d n)]
(if (= t 1) false ; probably prime
(loop [t t, s s]
(if (zero? s) true ; composite
(if (= t (- n 1)) false ; probably prime
(recur (powerMod t 2 n) (- s 1))))))))))
(if (= n 2) true
(loop [i 25] ; arbitrary number of witness trials
(if (zero? i) true ; probably prime
(let [a (+ 2 (rand-int (min 100000 (- n 2))))]
(if (witness? n a) false ; composite
(recur (- i 1))))))))
(println "2 is" (if (prime? 2) "prime" "composite"))
(println "87 is" (if (prime? 87) "prime" "composite"))
(println "2^89-1 is" (if (prime? 618970019642690137449562111) "prime" "composite"))
(defn factors "pollard-rho factorization" [n]
(defn rho [n c]
(defn f [x] (mod (+ (* x x) c) n))
(loop [h 5, t 2, d 1]
(if (< 1 d) d
(recur (f (f h)) (f t) (gcd (- t h) n)))))
(if (<= -1 n 1) (list n)
(if (neg? n) (cons -1 (factors (- n)))
(loop [n n, c 1, fs (list)]
(if (= n 1) fs (if (prime? n) (sort < (cons n fs))
(if (even? n) (recur (/ n 2) c (cons 2 fs))
(let [d (rho n c)]
(if (prime? d)
(recur (/ n d) (+ c 1) (cons d fs))
(recur n (+ c 1) fs))))))))))
(println "factors of 13290059 are" (factors 13290059))