Sieve Of Atkin
February 12, 2010
Our solution follows the one in Wikipedia:
(define (atkin limit)
(let ((v (make-vector (+ limit 1) #f)))
(define (flip k)
(vector-set! v k
(not (vector-ref v k))))
(do ((x 1 (+ x 1))) ((> (* x x) limit))
(do ((y 1 (+ y 1))) ((> (* y y) limit))
(let* ((n (+ (* 4 x x) (* y y))) (nmod12 (modulo n 12)))
(if (and (<= n limit) (member nmod12 '(1 5))) (flip n)))
(let* ((n (+ (* 3 x x) (* y y))) (nmod12 (modulo n 12)))
(if (and (<= n limit) (= nmod12 7)) (flip n)))
(let* ((n (- (* 3 x x) (* y y))) (nmod12 (modulo n 12)))
(if (and (> x y) (<= n limit) (= nmod12 11)) (flip n)))))
(do ((n 5 (+ n 1))) ((> (* n n) limit))
(if (vector-ref v n)
(do ((k (* n n) (+ k (* n n)))) ((> k limit))
(vector-set! v k #f))))
(let ((ps '(3 2)))
(do ((n 5 (+ n 1))) ((> n limit) (reverse ps))
(if (vector-ref v n) (set! ps (cons n ps)))))))
This implementation isn’t particularly fast. Dan Bernstein, co-creator of the sieve, has a screamingly fast implementation, much faster than any Sieve of Eratosthenes, at http://cr.yp.to/primegen.html. You can run our implementation at http://programmingpraxis.codepad.org/WAgoUQGS.
Pages: 1 2
The equations that are like: x² + y² = n
should be: x² + y² = k
Thanks. Fixed.
According to the Wikipedia page the equation x^2 + y^2 = k should read 4*x^2 + y^2 = k. Your code is also using the factor 4.
Fixed again. Thanks.
Tightened up the loop limits a little bit. It’s about 50% slower than a plain sieve on getting the primes less than a million (0.96 sec vs. 0.63 sec).
from math import sqrt limit = 1000000 def atkins(limit=1000000): '''use sieve of atkins to find primes <= limit.''' primes = [False] * limit sqrt_limit = int( sqrt( limit ) ) x_limit = int( sqrt( ( limit + 1 ) / 2 ) ) for x in xrange( 1, x_limit ): xx = x*x for y in xrange( 1, sqrt_limit ): yy = y*y n = 3*xx + yy if n <= limit and n%12 == 7: primes[n] = not primes[n] n += xx if n <= limit and n%12 in (1,5): primes[n] = not primes[n] if x > y: n -= xx + 2*yy if n <= limit and n%12 == 11: primes[n] = not primes[n] for n in xrange(5,limit): if primes[n]: for k in range(n*n, limit, n*n): primes[k] = False return [2,3] + filter(primes.__getitem__, xrange(5, limit))escucho maravillas de la criba de atkin pero todavia no he visto un codigo que me funcione yo hice otro tipo de criba qe trabaja diferebte y quiero comparar yo saco hasta 600 millones la saco en 14 seg en visual c pero no se si esta criba es mejor o peor, trabajo con una lenovo n200