Sieve Of Atkin16 Mar 2017
For some reason, everybody is obsessed with primes. Well, not everybody. Rather mathematicians and someone else. And an interesting problem with primes is how to generate them fast. Really fast. One of the proposals to reach this goal, which is pretty recent, is the Sieve of Atkin.
This sieve has in my opinion too much mathematics concepts for someone who doesn’t have a good mathematical knowledge, like me. That’s why I’m going to try to explain it with “easy” terms.
First we are going to take all the invertible numbers in Z60. Z60 represents the natural numbers in module 60. So Z60 is set of the numbers from 0 to 59.
If a number is >= 60 you can convert it to Z60 by doing the mod 60 operation. And the invertible numbers in Zn are those numbers that have a modular multiplicative inverse. Which means that if you multiply the number in Zn by another number in Zn you obtain 1 in module n. And which numbers have inverse? The numbers in Zn that are coprime with n.
Two numbers are coprime if their greatest common divisor is 1. For example:
5 is the highest number that both 25 and 60 are divisible by.
The only common divisor of 7 and 60 is 1, so they are coprime.
In the case of Z60, the invertible numbers are:
1 7 11 13 17 19 23 29 31 37 41 43 47 49 53 59
Binary Quadratic Form
Now, we are going to generate numbers with the form n = 60k + s where k is just a natural number and s an invertible number in Z60. If we want to generate all the primes from 60 to a limit, k will take values from 1 to ⌊limit/60⌋-1.
For example we will generate (60+1), (60+7) … (60+59), (120+1) …
We know that all the numbers with the previous form are primes candidates. Each one will be prime if it satisfies some rules which are based on binary quadratic forms. The sieve uses 3 algorithms (named 3.1, 3.2, and 3.3) depending on the module of the s used to generate the number.
There are three cases depending on the form of s. The number can be prime if the equation associated with its form has an odd number of solutions:
- s ∈ 1+4Z Which means that s mod 4 = 1. It will be prime if the equation 4x2 + y2 = n has an odd number of solutions (x,y) where x,y > 0.
- s ∈ 1+6Z Which means that s mod 6 = 1. It will be prime if the equation 3x2 + y2 = n has an odd number of solutions (x,y) where x,y > 0.
- s ∈ 11+12Z Which means that s mod 12 = 11. It will be prime if the equation 3x2 - y2 = n has an odd number of solutions (x,y) where x > y > 0.
For example, the number 67 (60 + 7):
s=7 so it has the form 1 + 6Z. The corresponding equation is 3x2 + y2 = n. The solutions that can be found are: x=1 y=8. So there are an odd number of solutions. Therefore 67 is still a prime candidate.
The hard part of this is how to find the solutions for those equations. The paper exposes efficient methods for this task but I’m not going to explain them because are quite complex and they are not important to understand the philosophy of this algorithm.
Finally, if after the previous algorithm the number is still a prime candidate we have to check if the number is square free. If it is, the number is prime!
By the way, square free means that if you factorize the number you won’t find any factor repeated. For example:
12 is not square free because it’s equal to 2·2·3, and we can see that the factor 2 is repeated.
How do we know if a number is square free? Well, actually we don’t know any algorithm that check it in polynomial time, so you can check if the number is divisible by the square of any of the primes that you have already calculated.
I’m sure that if you are interested in this algorithm is because of its performance. Unfortunately, I can’t swear that this sieve is faster than the Sieve of Eratosthenes. I don’t pretend to analyze deeply the reasons, although if you are interested in knowing the details I recommend you to have a look at this answer in StackOverflow.
Remember that performance is not only about the number of operations that you have to do to obtain the result; is also about cache misses, instruction level parallelism, branch predictions and much more.
I have made an implementation of this algorithm focused on making the code illustrative and simple, forgetting about performance related issues that would make it harder to understand.
Implementation in C++
Implementation in Maple
I want to thank Dr. G.M. Diaz-Toca her support, especially for understanding the original paper which requires a good mathematical background. I also want to thank her for her implementation of this sieve in Maple.
A.O.L. Atkin, D.J. Bernstein, Prime sieves using binary quadratic forms, Math. Comp. 73 (2004), 1023-1030.