Primality testing- II (Sieve of Atkins)

In the previous post on primality testing, we discussed successive division (brute force method) and sieve of Eratosthenes for checking if a number is prime and generating all primes below a given number respectively. Lets move to the Sieve of Atkins algorithm today which is a modified version of sieve of Eratosthenes. Sieve of Atkins is something I wouldn't suggest for a programming contest for it is difficult to understand and remember. To be frank, there are parts to it that I don't understand why they are right but it it works pretty fast when coupled with a good implementation. Note that if implemented badly, this could perform worse! The purpose behind putting this article on the blog is to spread the word about the existence of such an algorithm and that the code can be used when needed ;) Lets check out the algorithm now. Some excerpts first:

1. The algorithm treats 2, 3 and 5 as special cases and just adds them to the set of primes to start with.


2. 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.


3. The algorithm talks in terms of modulo-sixty remainders.


4. 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 4x² + y² = 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.


5. 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 3x² + y² = n is odd and the number is squarefree.


6. 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 3x² − y² = 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!