gistfile1.clj

; The Sieve of Eratosthenes
;
; This is a simple example of how to calculate prime numbers.
;
; Typically the algorithm is described as making a list of numbers from 2 until the maximum number you
; want to test for primality(?).  You then take the first number that hasn't been crossed off (2 at this
; point), and cross off all of the numbers in the list that are multiples of it.  Then you repeat for the
; next number that hasn't been crossed off.  At the end all of the numbers that have not been crossed off
; are prime.
;
; Unforunately this isn't very good if you want to generate a large number of primes, because you have
; to pre-generate all of the numbers you want to test, then repeatedly iterate over that enormous list
; multiple times.
;
; What we need to do is do exactly what the name suggests: put each number through a sieve.  That sieve
; can vary as we find more prime numbers.
;
; We can do that with a priority queue that holds all of the primes we know about and their next
; composite number.  The queue is prioritised such that the head of the queue contains the prime with
; the lowest composite at the head.  Now all you need to do is check a number against the next
; composite at the head of the queue and take the appropriate action.
;
; Let's define some terms:
;   P = prime number
;   M = next composite for P
;   Q = the priority queue hold pairs (P M), and (P' M') is the head of this queue
;   N = the number we're testing for primality
;
; * If N<M' then N is prime: we create a pair of (N,2N) and push it into Q;
; * If N=M' then N is not prime: we increment M' by P' and push the pair back into Q;
; * If N>M' then we have hit a "dual case": we increment M' by P' and push the back into Q.
;
; A "dual case" is one where two primes P1 and P2 both have the same M, i.e. 2 & 3 both have 6.
(require '[clojure.data.priority-map :as priority-map])

; Notice that this is a `def` that is created by immediately calling a function that returns
; a lazy sequence.  This makes using the primes easier.
(defn primes-by-sieve-of-erathosthenes
  ((fn primes
    ; We know that 2 is the first prime and can start with N=3 and Q=[(2 4)]
    ([]            (cons 2 (primes 3 (priority-map/priority-map 2 4))))

    ; Subsequent primes come from the pair of N and Q
    ([n q] (let [[p m] (peek q)]
             (cond
               (< n m) (cons n (lazy-seq (primes (inc n) (assoc q n (+ n n)))))
               (= n m)                   (primes (inc n) (assoc q p (+ p m)))
               :else                     (primes n       (assoc q p (+ p m)))))))))  ; No change in N

; Cool facts:
(take 10 primes-by-sieve-of-erathosthenes) ; => (2 3 5 7 11 13 17 19 23 29)
(nth primes-by-sieve-of-erathosthenes 3000) ; => 27457