The Sieve of Atkin in C#

I have previously written about the Sieve of Eratosthenes, which is an algorithm for finding primes. This algorithm worked very well for most of the prime related Euler Problems. However, for one of them it just didn’t do it. Well, it did it, but it did it kind of slow. The problem was to calculate the sum of all the primes below two million, and with that algorithm this took close to 2 seconds on my machine. Not too long you might think, but compared to my other solutions it is ages and ages. So I started to see if I could find a more efficient algorithm to use for that problem.

I quickly found one called the Sieve of Atkin. The Atkin algorithm works similarly to the original Eratosthenes one, but has a more fancy way of sieving out numbers which means it can work a lot faster. However, this also means that it needs to work towards an upper limit and that it have to find all the primes up to that limit in advance. In other words it cannot work incrementally and lazily like my Eratosthenes implementation.

Anyways, without further ado, here is my C# implementation of the Sieve of Atkin:

1
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
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
public class Atkin : IEnumerable<ulong>
{
  private readonly List<ulong> primes;
  private readonly ulong limit;

  public Atkin(ulong limit)
  {
      this.limit = limit;
      primes = new List<ulong>();
  }

  public IEnumerator<ulong> GetEnumerator()
  {
    if (!primes.Any())
        FindPrimes();

      foreach (var p in primes)
          yield return p;
  }

  IEnumerator IEnumerable.GetEnumerator()
  {
      return GetEnumerator();
  }

  private void FindPrimes()
  {
      var isPrime = new bool[limit + 1];
      var sqrt = Math.Sqrt(limit);

      for (ulong x = 1; x <= sqrt; x++)
          for (ulong y = 1; y <= sqrt; y++)
          {
              var n = 4 * x * x + y * y;
              if (n <= limit && (n % 12 == 1 || n % 12 == 5))
                  isPrime[n] ^= true;

              n = 3 * x * x + y * y;
              if (n <= limit && n % 12 == 7)
                  isPrime[n] ^= true;

              n = 3 * x * x - y * y;
              if (x > y && n <= limit && n % 12 == 11)
                  isPrime[n] ^= true;
          }

      for (ulong n = 5; n <= sqrt; n++)
          if (isPrime[n])
          {
              var s = n * n;
              for (ulong k = s; k <= limit; k += s)
                  isPrime[k] = false;
          }
               
      primes.Add(2);
      primes.Add(3);
      for (ulong n = 5; n <= limit; n+=2)
          if (isPrime[n])
              primes.Add(n);
  }
}

I am not going to explain my code. Instead I will just tell you to look at the algorithm at Wikipedia and look at my code while reading the pseudocode provided there. Why? Because it is pretty much a direct copy of it, just translated into C# :)

As a usage example, we can sum up all the primes below a 1000 like so:

var sum = new Atkin(1000)
  .Aggregate((sum, x) => sum + x);

Fantastic huh? The even more fantastic part is that when I used this class instead of the Eratosthenes class to solve the problem I mentioned in the beginning, my solution went from taking almost 2 seconds down to around 160 milliseconds :D

This is just a very straight forward implementation though, so if you have any suggestions to how it can be improved, made faster, cleaner, more efficient, et cetera, please let me know by leaving a comment!

  • Greg

    In the line

    for (ulong n = 5; n &lt;= limit; n++)

    could you use

    for (ulong n = 5; n &lt;= limit; n += 2)

    so you aren't checking to see if even numbers are primes?

    • http://www.geekality.net Torleif

      According to my unit test it seems to work at least, so thanks! :)

  • Mark Lawrence

    Ok I don’t know c# that well, but in your nested x/y for loops you compute both x * x and y * y three times, and 3 * x * x twice, could you compute these once only by using temporary variables? The third computation of n is done, then you check that x > y, why not perform this check first, then do the computation?

    • http://www.geekality.net Torleif

      Well, I suppose, but I have a feeling that allocating the variables might be slower than just doing the calculation… have you tried it out?

      • yofreke

        (my code uses longs)
        long xy = x * x + y * y;
        var n = 4 * xy;

        n = 3 * xy;
        I tried this and it took about 1/3 less of the time with my setup.

        • Ajay

          I still dnt understand. how could this code works??
          i mean take an example of nubmer 65..

          4*2*2 + 7*7 =65….. markdes as true in first loop….
          after then…. it would never be unmarked and would be in output inspite of being a composite….

          admin plzz solve this issue,,,,

          • http://www.geekality.net Torleif

            Like I said, I am not going to try to explain why the algorithm works. But it does, and there are tests in my code that proves it. Just try it out for yourself :)

            You can read about the algorithm at wikipedia,

  • Tom

    For large limit values, your algorithm will fail when it tries to initialize the isPrime array. Try passing ulong.MaxValue when you construct your Atkin object, and you’ll see what I mean.

    • http://www.geekality.net Torleif

      Ooh, good catch. Hadn’t tried it with that large values. Seems like it will fail for two reasons. First of all, the array can only contain int32.MaxValue items. Secondly, you’ll run out of memory :P

      Do you have any suggestions on how to fix it? I have a feeling this more a flaw of the algorithm itself than my implementation though… but yeah, if you have any smart solutions please let me know!

      • Ben

        The list of primes takes more memory stored as a list than as the bool array it’s calculated in.

        It also requires memory for both to change from one to the other.

        I replaced primes with the isPrime array (moved to class-level), changed the IEnumerator to this:

        public IEnumerator<ulong> GetEnumerator() {
                    if (!primeFlag) {
                        FindPrimes();
                        isPrime[2] = true;
                        isPrime[3] = true;
                        primeFlag = true;
                    }
                    for (ulong i = 2;i<limit;i++) {
                        if (isPrime[i]){
                            yield return i;
                        }
                    }
                }

        If memory still remains a problem, there are a number of memory optimisations that can be made, some at the expense of performance (for example, only representing odd numbers in isPrime).

  • Pingback: To find the list of primes till N. | I Live to Seek…()

  • Meet Shah

    Please atleast explain how 65 is not included in the primes even when it is first set true for being prime. 4*2*2 + 7*7 = 65 , 65%12 = 5

    • http://www.geekality.net Torleif

      It’s switched off again in the next step I believe. Try following the algorithm from Wikipedia yourself by hand and it should be clear what’s going on :)

    • Jim

      The ^= operator is an exclusive-or-equals operator. isPrime[n] ^= true returns true if isPrime[n] was false, returns false if isPrime[n] was true. Basically it toggles the boolean value in isPrime[n] to the opposite value from what it had. 65 has an even number of solutions, it toggles the boolean twice. Once for (2, 7) and once for (4, 1) since 4*2*2 + 7*7 = 65 and 4*4*4 + 1*1 = 65. The first one toggles isPrime[65] to true, the second one toggles isPrime[65] back to false. The Atkin sieve considers only those numbers for which an ODD number of solutions to the quadratic exist. 65 is an example of where more than one solution exists, but there are an EVEN number of them. Of those with an ODD number of solutions, only the ones that don’t have factors that are perfect squares (i.e. 9, 25, 49, …) are prime. The perfect squares, 4, 16, 36 don’t need to be checked as factors because no even numbers have solutions in the quadratics.

  • billa1972

    Fun post!

    I found BitArrays to be pretty efficient. Here’s a sieve algorithm that runs pretty quick

    public static IEnumerable<int> Seive(int candidate)
            {
         BitArray sieveContainer = new BitArray(candidate + 1, true); //assume all number are prime to start.
         int marker = 2; //start
         int factor = 2; //start.

         sieveContainer[0] = false;//0 is not prime
         sieveContainer[1] = false;//1 is not prime

         while ((marker * marker) <= sieveContainer.Length)
         {
              while ((factor += marker) <= candidate)
              {
         sieveContainer[factor] = false;
         };

         while (!sieveContainer.Get(++marker)) ;
         factor = marker;
                }

         //Return only the "true" indices
         return GetPrimes(sieveContainer);
            }