- The algorithm treats 2, 3 and 5 as special cases and just adds them to the set of primes to start with.
- Like Sieve of Eratosthenes, we start with a list of numbers we want to investigate. Suppose we want to find primes <=100, then we make a list for [5,1000] . As explained in (1), 2,3 and 5 are special cases and 4 is not a prime.
- The algorithm talks in terms of modulo-sixty remainders.
- All numbers with modulo-sixty remainder 1, 13, 17, 29, 37, 41, 49, or 53 have a modulo-four remainder of 1. These numbers are prime if and only if the number of solutions to 4
*x*^{2}+*y*^{2}=*n*is odd and the number is squarefree. A square free integer is one which is not divisible by any perfect square other than 1. - All numbers with modulo-sixty remainder 7, 19, 31, or 43 have a modulo-six remainder of 1. These numbers are prime if and only if the number of solutions to 3
*x*^{2}+*y*^{2}=*n*is odd and the number is squarefree. - All numbers with modulo-sixty remainder 11, 23, 47, or 59 have a modulo-twelve remainder of 11. These numbers are prime if and only if the number of solutions to 3
*x*^{2}−*y*^{2}=*n*is odd and the number is squarefree.

The points (4-6) are not intuitive and hence this algorithm is a little difficult to understand. But, if you are interested, check out the proofs on Atkin's paper. With these points in mind now, lets look at the following pseudocode from the Wikipedia article:

// arbitrary search limit

limit ← 1000000

// initialize the sieve

for i in [5, limit]: is_prime(i) ← false

// put in candidate primes:

// integers which have an odd number of

// representations by certain quadratic forms

for (x, y) in [1, √limit] × [1, √limit]:

n ← 4x²+y²

if (n ≤ limit) and (n mod 12 = 1 or n mod 12 = 5):

is_prime(n) ← ¬is_prime(n)

n ← 3x²+y²

if (n ≤ limit) and (n mod 12 = 7):

is_prime(n) ← ¬is_prime(n)

n ← 3x²-y²

if (x > y) and (n ≤ limit) and (n mod 12 = 11):

is_prime(n) ← ¬is_prime(n)

// eliminate composites by sieving

for n in [5, √limit]:

if is_prime(n):

// n is prime, omit multiples of its square; this is

// sufficient because composites which managed to get

// on the list cannot be square-free

is_prime(k) ← false, k ∈ {n², 2n², 3n², ..., limit}

print 2, 3

for n in [5, limit]:

if is_prime(n): print n

And finally lets convert this easy to understand algorithm into a python code:

[sourcecode language="Python"]

#!/usr/bin/python2.7 -tt

import sys

from math import sqrt, ceil, pow

def primes_below(limit): #sieve of Atkins

sqroot= int(ceil(sqrt(limit)))

is_prime=[False]*limit

primes = [2,3]

for x in xrange(sqroot+1):

for y in xrange(sqroot+1):

# n = 4*i^2 + j^2

n = 4*int(pow(x, 2)) + int(pow(y,2))

if n <= limit and (n % 12 == 1 or n % 12 == 5):

is_prime[n]= not is_prime[n]

# n = 3*i^2 + j^2

n = 3*int(pow(x, 2)) + int(pow(y,2))

if n <= limit and n % 12 == 7:

is_prime[n]= not is_prime[n]

# n = 3*i^2 - j^2

n = 3*int(pow(x, 2)) - int(pow(y,2))

if n <= limit and x > y and n % 12 == 11:

is_prime[n]= not is_prime[n]

for x in range(5, limit, 2):

if(is_prime[x] == True):

for y in xrange(x*x, limit, x):

is_prime[y]= False

primes.append(x);

print primes

def main():

primes_below(100)

if __name__=='__main__':

main()

[/sourcecode]

The C++ code for the algorithm can be found here. Enjoy!