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.

About these ads

Pages: 1 2

6 Responses to “Sieve Of Atkin”

  1. Jon said

    The equations that are like: x² + y² = n
    should be: x² + y² = k

  2. programmingpraxis said

    Thanks. Fixed.

  3. Paul Hofstra said

    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.

  4. programmingpraxis said

    Fixed again. Thanks.

  5. Mike said

    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))
    
  6. ivan said

    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

Leave a Reply

Fill in your details below or click an icon to log in:

WordPress.com Logo

You are commenting using your WordPress.com account. Log Out / Change )

Twitter picture

You are commenting using your Twitter account. Log Out / Change )

Facebook photo

You are commenting using your Facebook account. Log Out / Change )

Google+ photo

You are commenting using your Google+ account. Log Out / Change )

Connecting to %s

Follow

Get every new post delivered to your Inbox.

Join 611 other followers

%d bloggers like this: