sieve

(defun sieve-odds (maximum)
  "Prime numbers sieve for odd numbers.
   Returns a list with all primes that are less or equal to maximum."
  ;; maxi represents the number of odd numbers that are >= 3 and <= maximum
  ;; it's equal to (floor (1- maximum) 2)
  ;; because there are (/ EVENX 2)      odd numbers that are <= EVENX
  ;;   so there are (floor (1- EVENX) 2) odd numbers that are >= 3 and <= EVENX
  ;; also there are (floor (1+ ODDX) 2) odd numbers <= ODDX
  ;;   so there are (floor ODDX 2)=(floor (1- ODDX) 2) odd numbers that are >= 3 and <= ODDX
  ;; thus maxi is equal to (floor (1- maximum) 2) or equivalently (ash (1- maximum) -1)
  ;; also, maxi*2+1 is the highest odd number smaller than or equal to maximum
  ;; stop is the bitarray index upperbound for the floored square root of maximum
  ;; i is the bitarray index of the
  ;; (i+1)th odd natural number (let's call this number oddi).
  ;; i=1 maps to oddi=3, i=2 maps to 5,...
  (loop :with maxi = (ash (1- maximum) -1)
        :with stop = (ash (isqrt maximum) -1)
        :with sieve = (make-array (1+ maxi) :element-type 'bit :initial-element 0)
        :for i :from 1 :to maxi
        :for odd-number = (1+ (ash i 1))
        :when (zerop (sbit sieve i))
          ;; the first bit (sbit sieve 0) is initialized to 0 so 3 will always be collected
          :collect odd-number :into values
        :when (<= i stop)
          ;; next we're gonna sift out the multiples of oddi that are <= stop
          :do (loop :for j :from (* i (1+ i) 2) :to maxi :by odd-number
                    ;; j is the bitarray index of the (j-1)th odd number
                    ;; (let's call this number oddj)
                    ;; on every step j grows by oddj
                    ;; so that it can reach through all the bitarray indexes of oddi's multiples
                    ;; j starts from (* 2 (* i (1+ i))), which is (expt oddi 2)'s index
                    ;; because (expt oddi 2) is the smallest multiple of oddi which hasn't been
                    ;; yet sifted out: all multiples in the form oddi*oddk where k<i were
                    ;; sifted in a previous iteration(when i was equal to k).
                    ;; So because (expt oddi 2) is the smallest multiple to be sifted out
                    ;; the outer loop must end when i > stop (=the index of (isqrt maximum)).
                    #|
                    Here's how to figure out the bitarray index of (expt oddi 2)
                    array index                i:  1 2 3 4  5  6  7  8  9 10 11 12
                    the corresponding value oddi:  3 5 7 9 11 13 15 17 19 21 23 25
                    An index of 2 corresponds to the value 1+2*2=5
                    An index of 3 corresponds to the value 1+3*2=7
                    So to get from    index 2 (number 5) to index 12 (5^2)
                    means to get from     1+2*2=5        to     1+12*2=25=(1+2*2)^2

                    Thus the index of the square of oddi is the x in this equation:
                    (1+i*2)^2=1+x*2
                    1+4i+4i^2=1+2x
                    x=(1+4i+4i^2-1)/2
                    x=4(i+i^2)/2
                    x=2(i+i^2)
                    x=2(i(1+i))
                    |#
                    ;; Therefore j starts from this 2(i(1+i)), goes up through all the multiples
                    ;; and sifts them out.
                    ;; The loop stops when it reaches the highest multiple
                    ;; of oddi with the index smaller or equal to maxi
                    ;; (that's the highest multiple indexed by our bitarray)
                    :do (setf (sbit sieve j) 1))
        :finally (return (cons 2 values))))