(defun sieve-odds (maximum)
"Sieve for odd numbers."
;; 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
;; also, maxi*2+1 is the highest odd number smaller than or equal to maximum
;; stop is the bitarray index upperbound of 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.
#|
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))))