Generación de números primos rápidos en Clojure
lisp primes (12)
He estado trabajando en resolver problemas de Project Euler en Clojure para mejorar, y ya he encontrado la generación de números primos un par de veces. Mi problema es que solo lleva demasiado tiempo. Esperaba que alguien pudiera ayudarme a encontrar una forma eficiente de hacerlo de una manera Clojure-y.
Cuando lo hice, lo forcé brutalmente. Eso fue fácil de hacer. Pero el cálculo de los números primos 10001 tomó 2 minutos de esta manera en un Xeon 2.33GHz, demasiado largo para las reglas, y demasiado largo en general. Aquí estaba el algoritmo:
(defn next-prime-slow
"Find the next prime number, checking against our already existing list"
([sofar guess]
(if (not-any? #(zero? (mod guess %)) sofar)
guess ; Then we have a prime
(recur sofar (+ guess 2))))) ; Try again
(defn find-primes-slow
"Finds prime numbers, slowly"
([]
(find-primes-slow 10001 [2 3])) ; How many we need, initial prime seeds
([needed sofar]
(if (<= needed (count sofar))
sofar ; Found enough, we''re done
(recur needed (concat sofar [(next-prime-slow sofar (last sofar))])))))
Al reemplazar next-prime-slow con una rutina más nueva que tomó en cuenta algunas reglas adicionales (como la propiedad 6n +/- 1) pude acelerar las cosas hasta aproximadamente 70 segundos.
Luego intenté hacer un colador de Eratóstenes en puro Clojure. No creo haber sacado todos los errores, pero me rendí porque era demasiado lento (incluso peor que el anterior, creo).
(defn clean-sieve
"Clean the sieve of what we know isn''t prime based"
[seeds-left sieve]
(if (zero? (count seeds-left))
sieve ; Nothing left to filter the list against
(recur
(rest seeds-left) ; The numbers we haven''t checked against
(filter #(> (mod % (first seeds-left)) 0) sieve)))) ; Filter out multiples
(defn self-clean-sieve ; This seems to be REALLY slow
"Remove the stuff in the sieve that isn''t prime based on it''s self"
([sieve]
(self-clean-sieve (rest sieve) (take 1 sieve)))
([sieve clean]
(if (zero? (count sieve))
clean
(let [cleaned (filter #(> (mod % (last clean)) 0) sieve)]
(recur (rest cleaned) (into clean [(first cleaned)]))))))
(defn find-primes
"Finds prime numbers, hopefully faster"
([]
(find-primes 10001 [2]))
([needed seeds]
(if (>= (count seeds) needed)
seeds ; We have enough
(recur ; Recalculate
needed
(into
seeds ; Stuff we''ve already found
(let [start (last seeds)
end-range (+ start 150000)] ; NOTE HERE
(reverse
(self-clean-sieve
(clean-sieve seeds (range (inc start) end-range))))))))))
Esto es malo. También causa desbordamientos de pila si el número 150000 es más pequeño. Esto a pesar del hecho de que estoy usando recurre. Esa puede ser mi culpa.
Luego intenté con un tamiz, usando métodos Java en Java ArrayList. Eso llevó bastante tiempo y memoria.
Mi último intento es un tamiz usando un hash-map Clojure, insertando todos los números en el tamiz y luego disociando los números que no son primos. Al final, toma la lista de claves, que son los números primos que encontró. Toma entre 10 y 12 segundos encontrar 10000 números primos. No estoy seguro de que esté completamente depurado todavía. También es recursivo (usando recur y loop), ya que intento ser Lispy.
Entonces, con este tipo de problemas, el problema 10 (resume todos los números primos por debajo de 2000000) me está matando. A mi código más rápido se le ocurrió la respuesta correcta, pero tardó 105 segundos en hacerlo, y necesitaba bastante memoria (le di 512 MB solo para no tener que preocuparme por ello). Mis otros algoritmos tardan tanto que siempre los detuve primero.
Podría usar un tamiz para calcular muchos primos en Java o C bastante rápido y sin usar tanta memoria. Sé que debo estar perdiendo algo en mi estilo Clojure / Lisp que está causando el problema.
¿Hay algo que estoy haciendo realmente mal? ¿Clojure es un poco lento con grandes secuencias? Al leer algunas de las discusiones de Euler sobre el proyecto, las personas calcularon los primeros 10000 primos en otros Lisps en menos de 100 milisegundos. Me doy cuenta de que la JVM puede ralentizar las cosas y Clojure es relativamente joven, pero no esperaría una diferencia de 100x.
¿Alguien puede iluminarme de una manera rápida para calcular los números primos en Clojure?
Aquí hay una implementación simple y agradable:
http://clj-me.blogspot.com/2008/06/primes.html
... pero está escrito para alguna versión previa a la versión 1.0 de Clojure. Consulte lazy_seqs en Clojure Contrib para obtener uno que funcione con la versión actual del idioma.
Empecé a usar Clojure, así que no sé si es bueno, pero esta es mi solución:
(defn divides? [x i]
(zero? (mod x i)))
(defn factors [x]
(flatten (map #(list % (/ x %)) (filter #(divides? x %) (range 1 (inc (Math/floor (Math/sqrt x))))))))
(defn prime? [x]
(empty? (filter #(and divides? (not= x %) (not= 1 %)) (factors x))))
(def primes
(filter prime? (range 2 java.lang.Integer/MAX_VALUE)))
(defn sum-of-primes-below [n]
(reduce + (take-while #(< % n) primes)))
Aquí hay un simple tamiz en Scheme:
http://telegraphics.com.au/svn/puzzles/trunk/programming-in-scheme/primes-up-to.scm
Aquí hay una carrera para números primos hasta 10,000:
#;1> (include "primes-up-to.scm")
; including primes-up-to.scm ...
#;2> ,t (primes-up-to 10000)
0.238s CPU time, 0.062s GC time (major), 180013 mutations, 130/4758 GCs (major/minor)
(2 3 5 7 11 13...
(defn sieve
[[p & rst]]
;; make sure the stack size is sufficiently large!
(lazy-seq (cons p (sieve (remove #(= 0 (mod % p)) rst)))))
(def primes (sieve (iterate inc 2)))
con un tamaño de pila de 10M, obtengo la primo 1001 en ~ 33 segundos en un macbook de 2.1Gz.
Vea el último ejemplo aquí: http://clojuredocs.org/clojure_core/clojure.core/lazy-seq
;; An example combining lazy sequences with higher order functions
;; Generate prime numbers using Eratosthenes Sieve
;; See http://en.wikipedia.org/wiki/Sieve_of_Eratosthenes
;; Note that the starting set of sieved numbers should be
;; the set of integers starting with 2 i.e., (iterate inc 2)
(defn sieve [s]
(cons (first s)
(lazy-seq (sieve (filter #(not= 0 (mod % (first s)))
(rest s))))))
user=> (take 20 (sieve (iterate inc 2)))
(2 3 5 7 11 13 17 19 23 29 31 37 41 43 47 53 59 61 67 71)
Así que acabo de comenzar con Clojure, y sí, esto surge mucho en Project Euler, ¿no? Escribí un algoritmo de división de prueba bastante rápido, pero realmente no escala demasiado antes de que cada carrera de divisiones se vuelva prohibitivamente lenta.
Así que comencé de nuevo, esta vez usando el método de tamiz:
(defn clense
"Walks through the sieve and nils out multiples of step"
[primes step i]
(if (<= i (count primes))
(recur
(assoc! primes i nil)
step
(+ i step))
primes))
(defn sieve-step
"Only works if i is >= 3"
[primes i]
(if (< i (count primes))
(recur
(if (nil? (primes i)) primes (clense primes (* 2 i) (* i i)))
(+ 2 i))
primes))
(defn prime-sieve
"Returns a lazy list of all primes smaller than x"
[x]
(drop 2
(filter (complement nil?)
(persistent! (sieve-step
(clense (transient (vec (range x))) 2 4) 3)))))
Uso y velocidad:
user=> (time (do (prime-sieve 1E6) nil))
"Elapsed time: 930.881 msecs
Estoy bastante contento con la velocidad: se está quedando sin un REPL ejecutándose en un MBP 2009. Es casi todo rápido porque evito por completo el Clojure idiomático y en vez me vuelvo como un mono. También es 4 veces más rápido porque estoy usando un vector transitorio para trabajar en el tamiz en lugar de permanecer completamente inmutable.
Editar: Después de un par de sugerencias / correcciones de errores de Will Ness ahora se ejecuta mucho más rápido.
Usando una matriz de Java
(defmacro loopwhile [init-symbol init whilep step & body]
`(loop [~init-symbol ~init]
(when ~whilep ~@body (recur (+ ~init-symbol ~step)))))
(defn primesUnderb [limit]
(let [p (boolean-array limit true)]
(loopwhile i 2 (< i (Math/sqrt limit)) 1
(when (aget p i)
(loopwhile j (* i 2) (< j limit) i (aset p j false))))
(filter #(aget p %) (range 2 limit))))
Uso y velocidad:
user=> (time (def p (primesUnderb 1e6)))
"Elapsed time: 104.065891 msecs"
Muy tarde para la fiesta, pero voy a dar un ejemplo, usando Java BitSets:
(defn sieve [n]
"Returns a BitSet with bits set for each prime up to n"
(let [bs (new java.util.BitSet n)]
(.flip bs 2 n)
(doseq [i (range 4 n 2)] (.clear bs i))
(doseq [p (range 3 (Math/sqrt n))]
(if (.get bs p)
(doseq [q (range (* p p) n (* 2 p))] (.clear bs q))))
bs))
Ejecutando esto en una Macbook Pro 2014 (2.3GHz Core i7), obtengo:
user=> (time (do (sieve 1e6) nil))
"Elapsed time: 64.936 msecs"
Me doy cuenta de que esta es una pregunta muy antigua, pero recientemente terminé buscando lo mismo y los enlaces aquí no fueron lo que estoy buscando (restringido a tipos funcionales tanto como sea posible, generando perezosamente cada cada primo que quiero) .
Me encontré con una buena implementación de F # , por lo que todos los créditos son suyos. Simplemente lo porté a Clojure:
(defn gen-primes "Generates an infinite, lazy sequence of prime numbers"
[]
(let [reinsert (fn [table x prime]
(update-in table [(+ prime x)] conj prime))]
(defn primes-step [table d]
(if-let [factors (get table d)]
(recur (reduce #(reinsert %1 d %2) (dissoc table d) factors)
(inc d))
(lazy-seq (cons d (primes-step (assoc table (* d d) (list d))
(inc d))))))
(primes-step {} 2)))
El uso es simple
(->>
(gen-primes)
(filter #(< % 2000000)
; do what you need with the stuff
)
Aquí hay otro enfoque que celebra Clojure''s Java interop . Esto lleva 374 ms en un Core 2 Duo de 2,4 Ghz (ejecutando un solo hilo). BigInteger#isProbablePrime implementación eficiente de Miller-Rabin en BigInteger#isProbablePrime tratara con la verificación de primalidad.
(def certainty 5)
(defn prime? [n]
(.isProbablePrime (BigInteger/valueOf n) certainty))
(concat [2] (take 10001
(filter prime?
(take-nth 2
(range 1 Integer/MAX_VALUE)))))
La certeza Miller-Rabin de 5 probablemente no sea muy buena para números mucho más grandes que esto. Esa certeza es igual a 96.875% que es primo ( 1 - .5^certainty )
Después de llegar a este hilo y buscar una alternativa más rápida que las que ya están aquí, me sorprende que nadie esté relacionado con el siguiente artículo de Christophe Grand :
(defn primes3 [max]
(let [enqueue (fn [sieve n factor]
(let [m (+ n (+ factor factor))]
(if (sieve m)
(recur sieve m factor)
(assoc sieve m factor))))
next-sieve (fn [sieve candidate]
(if-let [factor (sieve candidate)]
(-> sieve
(dissoc candidate)
(enqueue candidate factor))
(enqueue sieve candidate candidate)))]
(cons 2 (vals (reduce next-sieve {} (range 3 max 2))))))
Además de una versión perezosa:
(defn lazy-primes3 []
(letfn [(enqueue [sieve n step]
(let [m (+ n step)]
(if (sieve m)
(recur sieve m step)
(assoc sieve m step))))
(next-sieve [sieve candidate]
(if-let [step (sieve candidate)]
(-> sieve
(dissoc candidate)
(enqueue candidate step))
(enqueue sieve candidate (+ candidate candidate))))
(next-primes [sieve candidate]
(if (sieve candidate)
(recur (next-sieve sieve candidate) (+ candidate 2))
(cons candidate
(lazy-seq (next-primes (next-sieve sieve candidate)
(+ candidate 2))))))]
(cons 2 (lazy-seq (next-primes {} 3)))))
Basado en el comentario de Will, aquí está mi opinión sobre postponed-primes :
(defn postponed-primes-recursive
([]
(concat (list 2 3 5 7)
(lazy-seq (postponed-primes-recursive
{}
3
9
(rest (rest (postponed-primes-recursive)))
9))))
([D p q ps c]
(letfn [(add-composites
[D x s]
(loop [a x]
(if (contains? D a)
(recur (+ a s))
(persistent! (assoc! (transient D) a s)))))]
(loop [D D
p p
q q
ps ps
c c]
(if (not (contains? D c))
(if (< c q)
(cons c (lazy-seq (postponed-primes-recursive D p q ps (+ 2 c))))
(recur (add-composites D
(+ c (* 2 p))
(* 2 p))
(first ps)
(* (first ps) (first ps))
(rest ps)
(+ c 2)))
(let [s (get D c)]
(recur (add-composites
(persistent! (dissoc! (transient D) c))
(+ c s)
s)
p
q
ps
(+ c 2))))))))
Presentación inicial para comparación:
Aquí está mi intento de transferir este generador de números primos de Python a Clojure. El siguiente devuelve una secuencia infinitamente floja.
(defn primes
[]
(letfn [(prime-help
[foo bar]
(loop [D foo
q bar]
(if (nil? (get D q))
(cons q (lazy-seq
(prime-help
(persistent! (assoc! (transient D) (* q q) (list q)))
(inc q))))
(let [factors-of-q (get D q)
key-val (interleave
(map #(+ % q) factors-of-q)
(map #(cons % (get D (+ % q) (list)))
factors-of-q))]
(recur (persistent!
(dissoc!
(apply assoc! (transient D) key-val)
q))
(inc q))))))]
(prime-help {} 2)))
Uso:
user=> (first (primes))
2
user=> (second (primes))
3
user=> (nth (primes) 100)
547
user=> (take 5 (primes))
(2 3 5 7 11)
user=> (time (nth (primes) 10000))
"Elapsed time: 409.052221 msecs"
104743
editar:
Comparación de rendimiento, donde postponed-primes utiliza una cola de primos vista hasta ahora en lugar de una llamada recursiva a postponed-primes :
user=> (def counts (list 200000 400000 600000 800000))
#''user/counts
user=> (map #(time (nth (postponed-primes) %)) counts)
("Elapsed time: 1822.882 msecs"
"Elapsed time: 3985.299 msecs"
"Elapsed time: 6916.98 msecs"
"Elapsed time: 8710.791 msecs"
2750161 5800139 8960467 12195263)
user=> (map #(time (nth (postponed-primes-recursive) %)) counts)
("Elapsed time: 1776.843 msecs"
"Elapsed time: 3874.125 msecs"
"Elapsed time: 6092.79 msecs"
"Elapsed time: 8453.017 msecs"
2750161 5800139 8960467 12195263)