; 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))