sieve_comp.clj

(ns sieve-comp)

(defmacro measure
  [& body]
  `(do
     (System/gc)
     (let [start# (System/nanoTime)
           result# (do ~@body)]
       [result# (- (System/nanoTime) start#)])))

(defn relate
  ([measured-a measured-b] (relate true measured-a measured-b))
  ([strict [result-a elapsed-a] [result-b elapsed-b]]
   (when (or (not strict) (= result-a result-b))
     (float (/ elapsed-a elapsed-b)))))

(defn order-of-growth
  [sizes times]
  (->> [sizes times]
       (map #(partition 2 1 %))
       (apply map (fn [[n1 n2] [t1 t2]]
                    (/ (Math/log (/ t2 t1))
                       (Math/log (/ n2 n1)))))))

(defn measure-growth
  [f sizes]
  (->> sizes
       (map #(-> % f last measure second))
       (order-of-growth sizes)))

(def wheel-2-3-5-7
  (cycle [2 4 2 4 6 2 6 4 2 4 6 6 2 6 4 2 6 4 6 8 4 2 4 2 4 8 6 4 6 2 4 6 2 6 6
          4 2 4 6 2 6 4 2 4 2 10 2 10]))

(defn spin
  [[x & xs] n]
  (cons n (lazy-seq (spin xs (+ n x)))))

(defn skip-2-3-5-7 [] (concat [2 3 5 7] (spin wheel-2-3-5-7 11)))


(defn sieve
  "int -> int[]
  Returns a seq of primes up to 'limit inclusive."
  ([limit] (sieve limit (skip-2-3-5-7)))
  ([limit candidates]
   (letfn [(mark-multiples! [unmarked p]
             (->> (range (* p p) (inc limit) p)
                  (reduce #(assoc! %1 %2 nil) unmarked)))]
     (loop [[x & xs] candidates
            primes (->> (range 2 (inc limit))
                        (reduce conj! (transient [nil nil])))]
       (if (> (* x x) limit)
         (remove nil? (persistent! primes))
         (recur xs (if (primes x) (mark-multiples! primes x) primes)))))))


(defn sieve-array
  ([limit] (sieve-array limit (skip-2-3-5-7)))
  ([limit candidates]
   (let [primes (long-array (concat [0 0] (range 2 (inc limit))))]
     (doseq [x candidates
             :while (<= (* x x) limit)
             :when (pos? (aget primes x))
             i (range (* x x) (inc limit) x)]
       (aset primes i 0))
     (-> (areduce primes x xs (transient [])
                  (if (pos? (aget primes x)) (conj! xs x) xs))
         persistent!))))

(defn sieve-packed
  [limit]
  (let [last-bit (long (quot (- limit 3) 2))
        last-candidate-bit (long (quot (- (Math/sqrt limit) 3) 2))
        size (long (inc (quot limit 64)))
        composites (long-array size 0)]  ;; zero means prime

    (loop [bit 0]
      (when (<= bit last-candidate-bit)
        (let [w (bit-shift-right bit 6)
              offset (bit-and bit 63)
              word (aget composites w)
              p (+ 3 bit bit)]
          (when (zero? (bit-and word (bit-shift-left 1 offset)))
            (loop [m (+ bit (* p (inc bit)))]
              (when (<= m last-bit)
                (let [m-bit (bit-shift-right m 6)
                      m-offset (bit-and m 63)
                      m-word (aget composites m-bit)]
                  (aset composites m-bit
                        ^long (bit-or m-word (bit-shift-left 1 m-offset))))
                (recur (long (+ p m)))))))
        (recur (inc bit))))

    (letfn [(produce [^long bit]
              (when (<= bit last-bit)
                (lazy-seq
                  (let [w (bit-shift-right bit 6)
                        offset (bit-and bit 63)
                        word (aget composites w)]
                    (if (zero? (bit-and word (bit-shift-left 1 offset)))
                      (cons (+ 3 bit bit) (produce (inc bit)))
                      (produce (inc bit)))))))]
      (cons 2 (produce 0)))))

(defn sieve-set
  [limit]
  (loop [[x & xs] (iterate inc 2)
         composites (transient #{})]
    (if (> (* x x) limit)
      (remove composites (range 2 (inc limit)))
      (recur xs (if (composites x)
                  composites
                  (->> (range (* x x) (inc limit) x)
                       (reduce conj! composites)))))))

(defn sieve-bit-packed
  "from http://rosettacode.org/wiki/Sieve_of_Eratosthenes#Clojure
  with slightly cleaned up names."
  [n]
  (let [root (long (Math/sqrt n))
        root-index (long (/ (- root 3) 2))
        index (long (/ (- n 3) 2))
        composites (long-array (inc (/ index 64)))
        prime? #(zero? (bit-and (aget composites (bit-shift-right % 6))
                                (bit-shift-left 1 (bit-and % 63))))
        cull-primes (fn [i]
                      (let [p (long (+ i i 3))]
                        (loop [i (bit-shift-right (- (* p p) 3) 1)]
                          (if (<= i index)
                            (let [w (bit-shift-right i 6)]
                              (aset composites w
                                    (bit-or (aget composites w)
                                            (bit-shift-left 1 (bit-and i 63))))
                              (recur (+ i p)))))))
        cull (fn []
               (loop [i 0]
                 (when (<= i root-index)
                   (when (prime? i)
                     (cull-primes i))
                   (recur (inc i)))))]
    (letfn [(next-prime [i]
              (if (<= i index)
                (cons (+ i i 3) (lazy-seq (next-prime (loop [i (inc i)]
                                                        (if (or (> i index)
                                                                (prime? i))
                                                          i
                                                          (recur (inc i)))))))))]
      (when (>= n 2)
        (cons 2 (when (>= n 3)
                  (cull)
                  (lazy-seq (next-prime 0))))))))