Daniel Lyons <[email protected]> writes: Hi Daniel,
> In your particular case I think you might want to take a look at the > source code in the contrib project which implements a lazy list of > primes: > > http://code.google.com/p/clojure-contrib/source/browse/trunk/src/clojure/contrib/lazy_seqs.clj I came up with another approach on lazily calculating primes using the Miller Rabin algorithm. Calculating the primes < 100.000 akes about 12 seconds on my 2.2 GHz dual core. Here it is: --8<---------------cut here---------------start------------->8--- (ns de.tsdh.math.primes (:use [clojure.contrib [math :only [expt]] [test-is :only [deftest is]] [def :only [defvar defvar-]]])) (defvar *pseudo-accuracy* 10 "The chances that prime? returns true for a composite is (expt 4 (* -1 *pseudo-accuracy*)).") (defn factorize-out "Factorizes out all x factors from n. Examples: (factorize-out 10 2) ==> {:exponent 1, :rest 5}, because 2^1 * 5 = 10 (factorize-out 90 3) ==> {:exponent 2, :rest 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-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 "Returns the next prime greater than n." [n] (cond (< n 2) 2 :else (loop [n (inc n)] (if (prime? n) n (recur (inc n)))))) (defn primes "Returns a lazy sequence of prime numbers" [] (iterate next-prime 2)) --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 [email protected] Note that posts from new members are moderated - please be patient with your first post. To unsubscribe from this group, send email to [email protected] For more options, visit this group at http://groups.google.com/group/clojure?hl=en -~----------~----~----~----~------~----~------~--~---
