Wednesday, November 20, 2013

The Sieve of Atkin in Common Lisp

So since I've been learning a bit of common lisp lately, I've decided to implement one of my favorite go-to algorithms for generating a list or primes: the sieve of Atkin. I won't go to in depth on the underlying mathematics of the algorithm (Wikipedia does a pretty good job as it is), but it basically sieves primes based on the solutions to a few specific quadratic forms. As of right now, the code doesn't work as well as I'd like. The speed seems fine, but generating a list of primes over 5000 gives me a stack overflow (whereas the C++ version I implemented earlier can generate primes under 1,000,000+). I'm thinking the recursive functions are the main culprits, even though I tried to make them tail-recursive in order to avoid this very issue . . . Anyway, if I find a way to optimize the code any better I'll post it here later.

;;;; Generates a list of primes using the "Sieve of Atkin."
;;;; More on this algorithm can be found at:
;;;; http://en.wikipedia.org/wiki/Sieve_of_Atkin

(setf *limit* 5000)    ; limit for maximum prime in list


;;; generate list of "nil" of size lim
(defun null-list (lim) 
  (nullness lim '()))

(defun nullness (lim lst)
  (if (eql lim 0)
    lst
    (nullness (- lim 1) (cons nil lst))))


;;; sieve the primes
(defun sieve (lim)
  (setf is_prime (null-list lim))
  (do ((x 1 (+ x 1)))
      ((> x (+ (round (sqrt lim)) 1)) is_prime)
    (do ((y 1 (+ y 1)))
        ((> y (+ (round (sqrt lim)) 1)) is_prime)

      (setf n (+ (* 4 x x) (* y y)))
      (if (and (<= n lim) (or (eql (mod n 12) 1) (eql (mod n 12) 5)))
        (setf (nth n is_prime) (not (nth n is_prime))))

      (setf n (+ (* 3 x x) (* y y)))
      (if (and (<= n lim) (eql (mod n 12) 7))
        (setf (nth n is_prime) (not (nth n is_prime))))

      (setf n (- (* 3 x x) (* y y)))
      (if (and (> x y) (<= n lim) (eql (mod n 12) 11))
        (setf (nth n is_prime) (not (nth n is_prime)))))))


;;; set the square multiples of primes to false
(defun square-sieve (lst lim)
  (do ((n 0 (+ n 1)))
      ((> n (+ (round (sqrt lim))) 1) lst)
    (if (eql (nth n lst) 't)
      (do ((k (* n n) (+ k (* n n))))
          ((> k (- lim 1)) lst)
        (setf (nth k lst) nil)))))


;;; print the primes
(defun print-primes (lst lim)
  (prime-print-print lst 0 '()))

(defun prime-print-print (lst n primes)
  (if (eql n (list-length lst))
    (cons 2 (cons 3 (reverse primes)))
    (if (eql (nth n lst) 't)
      (prime-print-print lst (+ n 1) (setf primes (cons n primes)))
      (prime-print-print lst (+ n 1) primes))))


;;; print the primes for the above limit
(print (print-primes (square-sieve (sieve *limit*) *limit*) *limit*))

No comments:

Post a Comment