;;;; 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*))
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.
Subscribe to:
Post Comments (Atom)
No comments:
Post a Comment