sieve.clj

(ns sieve.core
  (:use clojure.data.priority-map)
  (:gen-class))

(defn doto-last [coll f & args]
  (let [idx (dec (count coll))
        args (conj args (coll idx))]
    (assoc coll idx (apply f args))))

(defn prime-gaps [gk-1]
  "return the next constealltion of gaps between primes given the previous
  reference: http://staff.washington.edu/fbholt/papers/primegaps.pdf"
  (let [k (inc (first gk-1))]
    (loop [in (drop 2 (flatten (repeat k gk-1)))
           cum-diff (map #(* %1 k) gk-1)
           sum (nth gk-1 1)
           out [(+ (first gk-1) sum)]]
      (if (empty? in)
        (sequence out)
        (let [n (first in)
              t (first cum-diff)
              in (rest in)]
          (if (= sum t)
            (recur in (next cum-diff) n (doto-last out + n))
            (recur in cum-diff (+ sum n) (conj out n))))))))

(defn prime-cand
  ([]
   "quick and dirty wheel for testing, an infinite wheel using G(7)"
   (let [g3 '(4 2)
         w (cycle (prime-gaps (prime-gaps g3)))] ;; g7
     (prime-cand (+ 1 (first w)) (rest w))))
  ([n w]
   (cons n (lazy-seq (prime-cand (+ n (first w))
                                 (rest w))))))

(defn adjust [comp-map n]
  "enforce invariant: n <= all enties in the map, allows wheel optimization"
  (let [k (key (peek comp-map))
        v (val (peek comp-map))]
    (if (<= n v)
      comp-map
      (recur (assoc comp-map k (+ k v)) n))))

(defn sieve [candidates comp-map]
   (let [n (first candidates)
         comp-map (adjust comp-map n)
         v (val (peek comp-map))]
     (if (= n v)
       (recur (rest candidates) comp-map)
       (cons n (lazy-seq (sieve (rest candidates)
                                (assoc comp-map n (* n n))))))))

(defn primes
  "yeilds an unbounded sequence of primes with a priotiry map (with wheel optimization)
  reference: https://www.cs.hmc.edu/~oneill/papers/Sieve-JFP.pdf"
  ([] (conj (primes (prime-cand)) 7 5 3 2))
  ([cand]
   (let [n (first cand)]
     (cons n (lazy-seq (sieve (rest cand)
                              (priority-map n (* n n))))))))

(defn -main [& args]
  (println (take 5000 (primes))))