Author: aaa. Date 2018-03-22 08:15:16, views: 354, Raw
(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))))