Mark Engelberg <mark.engelb...@gmail.com> writes:
Hi Mark,
>> For number theory you often need things like
>>
>> (mod (expt n exp) m)
>
> Yes, and to make things like this fast, the trick is to perform the
> mod at each multiplication step, rather than waiting until the end.
So now I added this suggestion and the test is now really, really fast.
I tried some really big primes from the list of prime numbers on
wikipedia, and it always finishes instantly. And what's best: the
results seem to be correct, too. ;-)
Here's the code. Please feel free to comment on where I could do better
or use a more ideomatic style.
--8<---------------cut here---------------start------------->8---
(ns de.tsdh.clojure.math.primes
;; TODO: Only expt from math is needed, but how do I use :only with multiple
;; libs?
(:use [clojure.contrib math test-is]))
(def
#^{:doc "The chances that prime? returns true for a composite if *pseudo*
is true is (expt 4 (* -1 *pseudo-accuracy*))."}
*pseudo-accuracy* 10)
(defn factorize-out
"Factorizes out all x factors from n.
Examples:
(factorize-out 10 2) ==> 5, because 2^1 * 5 = 10
(factorize-out 90 3) ==> 10, because 3^2 * 10 = 90"
[n x]
(loop [acc n exp 0]
(if (= 0 (mod acc x))
(recur (/ acc x) (inc exp))
(hash-map :exponent exp :rest acc))))
(defn expt-mod
"Equivalent to (mod (expt n e) m), but faster.
http://en.wikipedia.org/wiki/Modular_exponentiation#An_efficient_method:_the_right-to-left_binary_algorithm";
[n e m]
(loop [r 1, b n, e e]
(if (= e 0)
r
(recur (if (odd? e)
(mod (* r b) m)
r)
(mod (expt b 2) m)
(bit-shift-right e 1)))))
(defn prime?
"Checks if n is a prime using the Miller-Rabin pseudo-primality test. Also
see *pseudo* and *pseudo-accuracy*."
[n]
(cond
(< n 2) false
(= n 2) true
(even? n) false
:else (let [m (factorize-out (dec n) 2)
d (:rest m)
s (:exponent m)]
(loop [k 1]
(if (> k *pseudo-accuracy*)
true
(let [a (+ 2 (rand-int (- n 4)))
x (expt-mod a d n)]
(if (or (= x 1) (= x (dec n)))
(recur (inc k))
(if (loop [r 1
x (expt-mod x 2 n)]
(cond
(or (= x 1) (> r (dec s))) false
(= x (dec n)) true
:else (recur (inc r) (mod (* x x) n))))
(recur (inc k))
false))))))))
(defn next-prime [n]
"Returns the next prime greater than n."
(cond
(< n 2) 2
:else (loop [n (inc n)]
(if (prime? n)
n
(recur (inc n))))))
;; test stuff
(def
#^{:private true}
first-1000-primes
'( 2 3 5 7 11 13 17 19 23 29
;; [Snipped, but the test works with all 1000 primes...]
7841 7853 7867 7873 7877 7879 7883 7901 7907 7919))
(deftest test-prime-fns
(loop [a (take 1000 (iterate next-prime 2))
b (take 1000 first-1000-primes)]
(when (and (empty? a) (empty? b))
(is (= (first a) (first b)))
(recur (rest a) (rest b)))))
--8<---------------cut here---------------end--------------->8---
Bye,
Tassilo
--~--~---------~--~----~------------~-------~--~----~
You received this message because you are subscribed to the Google Groups
"Clojure" group.
To post to this group, send email to clojure@googlegroups.com
To unsubscribe from this group, send email to
clojure+unsubscr...@googlegroups.com
For more options, visit this group at
http://groups.google.com/group/clojure?hl=en
-~----------~----~----~----~------~----~------~--~---
