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!

Primality testing- I (brute force, Sieve of Eratosthenes)

If you are into algorithms, you would realize that Primality testing is a very important topic. In simple terms, Primality testing means finding if a number is prime or not. Lets discuss a few tricks for the same. There is actually a lot of number theory based on primality testing and it is impossible to cover all of them here. Nevertheless lets discuss some of them. This post starts with a simple brute-force algorithm and introduces Sieve of Eratosthenes. With these simple things covered, we will try and discuss more of primality testing in the following posts. Lets get started with the brute force approach!

NAIVE ALGORITHM:

Lets say the given number is n. The simplest(brute-force) way to find if a number is prime or not is to check if the number is divisible by any number in the range 2 to n-1. This follows from the definition of a prime number and should be easy to understand. At this moment, lets revisit a basic school grade mathematics property which says:

If a number, say n, is divisible by an integer k, then either k or n/k has to be <= sqrt(n)

This means you dont need to check if n is divisible by any of the numbers in [0, n-1], rather you can just check for [0, sqrt(n)]. Well, Give yourself a pat on the back because that seems like a good enough optimization. If you were making 10000 checks in the former case, you have reduced it to just 100 iterations with the addition of this property.  Lets have a look at a small python code to illustrate this (source: github):

[sourcecode language="python"]

def isPrime(n):
maxPossibleFactor = int( floor( sqrt(n) ) )
possibleFactors = range( 2, maxPossibleFactor+1 )
divisibleFactors = [ n % i == 0 for i in possibleFactors ]
return not any( divisibleFactors )

[/sourcecode]

That was easy. Right?

SIEVE OF ERAOSTHENES:

Lets move ahead and look at a slightly different question: Identify all the prime numbers below 1000. A naive way to do this is to call the above python function for [0-1000].  Sounds like bad... Can we do better? Yes, we can. Turns out there is a well known algorithm for doing this- Sieve of Eratosthenes. Lets see how this algorithm works (The explanation below is  copied from here):

Make a list of all the integers less than or equal to n (and greater than one). Strike out the multiples of all primes less than or equal to the square root of n, then the numbers that are left are the primes.

For example, to find all the primes less than or equal to 30, first list the numbers from 2 to 30.

2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30

The first number 2 is prime, so keep it (we will color it green) and cross out its multiples (we will color them red), so the red numbers are not prime.

2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30

The first number left (still black) is 3, so it is the first odd prime. Keep it and cross out all of its multiples. We know that all multiples less than 9 (i.e. 6) will already have been crossed out, so we can start crossing out at 3²=9.

2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30

Now the first number left (still black) is 5, the second odd prime. So keep it also and cross out all of its multiples (all multiples less than 5²=25 have already been crossed out, and in fact 25 is the only multiple not yet crossed out).

2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30

The next number left, 7, is larger than the square root of 30, so there are no multiples of 7 to cross off that haven't already been crossed off (14 and 28 by 2, and 21 by 3), and therefore the sieve is complete. Therefore all of the numbers left are primes: {2, 3, 5, 7, 11, 13, 17, 19, 23, 29}. Notice we just found these primes without dividing.

Now, lets look at a fun animation to illustrate the algorithm (source: wikipedia):


With the algorithm clear now, lets look at a python code (source: stackoverflow).

[sourcecode language="python"]
def primes_sieve2(limit):
a = [True] * limit # Initialize the primality list
a[0] = a[1] = False

for (i, isprime) in enumerate(a):
if isprime:
for n in xrange(i*i, limit, i): # Mark factors non-prime
a[n] = False
[/sourcecode]

If you have a trouble understanding some part of  this post or any other related question, please drop a comment. And by the time, we have the next post ready, you may try some questions related to prime numbers at http://projecteuler.net/problems. Happy coding :)