fork download
copy 
  1. ; modest prime-number library
  2.  
  3. (defn gcd "greatest common divisor" [a b]
  4.   (if (zero? b) a (gcd b (mod a b))))
  5.  
  6. (println "gcd of 8 and 12 is" (gcd 8 12))
  7.  
  8. (defn powerMod "modular exponentiation" [b e m]
  9.   (defn m* [p q] (mod (* p q) m))
  10.   (loop [b b, e e, x 1]
  11.     (if (zero? e) x
  12.       (if (even? e) (recur (m* b b) (/ e 2) x)
  13.         (recur (m* b b) (quot e 2) (m* b x))))))
  14.  
  15. (def b 2988348162058574136915891421498819466320163312926952423791023078876139)
  16. (def e 2351399303373464486466122544523690094744975233415544072992656881240319)
  17. (def m 10000000000000000000000000000000000000000)
  18. (println "b ^ e (mod m) =" (powerMod b e m))
  19.  
  20. (defn primes "primes less than n" [n]
  21.   (let [sieve (boolean-array n true)]
  22.     (loop [p 2, ps (list)]
  23.       (cond (= n p) (reverse ps)
  24.             (aget sieve p)
  25.               (do (doseq [i (range (* p p) n p)]
  26.                     (aset sieve i false))
  27.                   (recur (+ p 1) (cons p ps)))
  28.             :else (recur (+ p 1) ps)))))
  29.  
  30. (println "primes less than 30 are" (primes 30))
  31.  
  32. (defn prime? "miller-rabin primality test" [n]
  33.   (defn witness? [n a]
  34.     (loop [d (- n 1), s 0]
  35.       (if (even? d) (recur (/ d 2) (+ s 1))
  36.         (let [t (powerMod a d n)]
  37.           (if (= t 1) false ; probably prime
  38.             (loop [t t, s s]
  39.               (if (zero? s) true ; composite
  40.                 (if (= t (- n 1)) false ; probably prime
  41.                   (recur (powerMod t 2 n) (- s 1))))))))))
  42.   (if (= n 2) true
  43.     (loop [i 25] ; arbitrary number of witness trials
  44.       (if (zero? i) true ; probably prime
  45.         (let [a (+ 2 (rand-int (min 100000 (- n 2))))]
  46.           (if (witness? n a) false ; composite
  47.             (recur (- i 1))))))))
  48.  
  49. (println "2 is" (if (prime? 2) "prime" "composite"))
  50. (println "87 is" (if (prime? 87) "prime" "composite"))
  51. (println "2^89-1 is" (if (prime? 618970019642690137449562111) "prime" "composite"))
  52.  
  53. (defn factors "pollard-rho factorization" [n]
  54.   (defn rho [n c]
  55.     (defn f [x] (mod (+ (* x x) c) n))
  56.     (loop [h 5, t 2, d 1]
  57.       (if (< 1 d) d
  58.         (recur (f (f h)) (f t) (gcd (- t h) n)))))
  59.   (if (<= -1 n 1) (list n)
  60.     (if (neg? n) (cons -1 (factors (- n)))
  61.       (loop [n n, c 1, fs (list)]
  62.         (if (= n 1) fs (if (prime? n) (sort < (cons n fs))
  63.           (if (even? n) (recur (/ n 2) c (cons 2 fs))
  64.             (let [d (rho n c)]
  65.               (if (prime? d)
  66.                   (recur (/ n d) (+ c 1) (cons d fs))
  67.                   (recur n (+ c 1) fs))))))))))
  68.  
  69. (println "factors of 13290059 are" (factors 13290059))
Success #stdin #stdout 1.54s 389120KB
comments (?)
stdin
copy 
Standard input is empty
stdout
copy 
gcd of 8 and 12 is 4
b ^ e (mod m) = 1527229998585248450016808958343740453059N
primes less than 30 are (2 3 5 7 11 13 17 19 23 29)
2 is prime
87 is composite
2^89-1 is prime
factors of 13290059 are (3119 4261)
https://ideone.com/3zRNil
language:
Clojure (clojure 1.10.0)
created:
10 years ago
visibility:
public

Share or Embed source code

Discover > Sphere Engine API

The brand new service which powers Ideone!

Discover > IDE Widget

Widget for compiling and running the source code in a web browser!