#! /usr/bin/env shell-sbcl

;; $Id: primes-sieve.lisp,v 1.2 2006/12/18 21:39:23 mca1001 Exp $

(defconstant +max-prime+ 400000000)
(declaim (optimize (speed 2)))
(deftype prime-t () `(integer 0 ,+max-prime+))
;(deftype prime-t () `(integer))
(deftype prime-flag-t () `(or prime-t null))

(defun is-prime (num bits)
  ;; assume the number is valid inside bits, and bits are initialise
  (declare (prime-t num)
	   (simple-bit-vector bits))
  (if (= 1 (bit bits num)) t))


(defmacro do-primes ((var start end bits &key stop-if) &body body)
  (declare (symbol var))
  (let ((start-var (gensym))
	(end-var (gensym))
	(bits-var (gensym)))
    `(let ((,start-var ,start)
	   (,end-var ,end)
	   (,bits-var ,bits))
      (declare (prime-t ,start-var ,end-var)
       (simple-bit-vector ,bits-var))
      (do ((,var ,start-var (1+ ,var)))
	  ((or ,stop-if
	       (>= ,var ,end-var)))
	(declare (prime-t ,var))
	(when (is-prime ,var ,bits-var)
	  ,@body)))))

(defun gen-primes (max)
  (declare (prime-t max))
  (let ((primes (make-array (1+ max)
			    :element-type 'bit :initial-element 1)))
    (loop for i from 2 to max
	  do
	  (if (is-prime i primes)
	      ;; tick off multiple of i
	      (do ((j (* i 2) (+ j i)))
		  ((>= j max))
		;; declaring j : could reach almost (* max 2)
		#+noisy
		(if (is-prime j primes)
		    (format t "tickoff ~s because ~s~%" j i))
		(setf (bit primes j) 0))))
    primes))

(defun gen-primes-broken (max)
  (declare (prime-t max))
  (let ((primes (make-array (1+ max)
			    :element-type 'bit :initial-element 1)))
    ;; tick off non-prime evens
    (do ((j 4 (+ j 2)))
	((>= j max))
      (declare (prime-t j))
      (setf (bit primes j) 0))
    (loop for i from 3 to max
	  do
	  (if (is-prime i primes)
	      ;; tick off powers of i
;;; need to unprime also all combinations of all primes to this point
;;; excepting where the product is clearly > max, e.g. two factors > sqrt(max)
	      (do ((j (* i i) (* j i)))
		  ((>= j max))
		(declare (prime-t j))
		(setf (bit primes j) 0))))
    primes))

(defun list-primes (start end bits)
  (let (rev-out)
    (do-primes (i start end bits)
      (setf rev-out (cons i rev-out)))
    (reverse rev-out)))

(let ((primes (time (gen-primes +max-prime+))))

  #+ doh
  (let ((l1 (list-primes 2 +max-prime+ primes))
	(l2 (list-primes 2 +max-prime+
			 (gen-primes2 +max-prime+))))
    (assert (eql l1 l2)
	    (l1 l2)
	    "l1 = ~s~%  l2 = ~s~%" l1 l2))

;(format t "list =~%~s~%" (list-primes 2 +max-prime+ primes))


  (let ((num 0))
    (declare (prime-t num))
(time
    (do-primes (i 2 +max-prime+ primes)
      (incf num))
)
    (format t "Of 1 .. ~s found ~s primes = ~f%~%"
	    +max-prime+ num (/ num +max-prime+ 0.01)))

  )