C#: Implementation of the Sieve of Atkin

Question

I was wondering if someone here have a good implementation of the Sieve of Atkin that they would like to share.

I am trying to implement it, but can't quite wrap my head around it. Here is what I have so far.

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>();
    }

    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])
                for (ulong k = n*n; k <= limit; k *= k)
                    isPrime[k] = false;

        primes.Add(2);
        primes.Add(3);
        for (ulong n = 5; n <= limit; n++)
            if (isPrime[n])
                primes.Add(n);
    }


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


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


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

I have pretty much just tried to "translate" the pseudocode listed at Wikipedia, but it isn't working correctly. So either I have misunderstood something or just done something wrong. Or most likely both...

Have a list of the first 500 primes which I use as a test and my implementation fails at number 40(or 41?).

Values differ at index [40]
Expected: 179
But was: 175

Are you able to find my mistake, do you have an implementation laying around that you could share, or both?

---

The exact test I am using looks like this:

public abstract class AtkinTests
{
    [Test]
    public void GetEnumerator_FirstFiveHundredNumbers_AreCorrect()
    {
        var sequence = new Atkin(2000000);
        var actual = sequence.Take(500).ToArray();
        var expected = First500;

        CollectionAssert.AreEqual(expected, actual);
    }

    private static readonly ulong[] First500 = new ulong[]
        {
            2, 3, 5, 7, 11, 13, 17, ...
        };
}

Answer 1

This code:

for (ulong k = n*n; k <= limit; k *= k)
  isPrime[k] = false;

doesn't seem to be a faithful translation of this pseudocode:

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

Your code looks like it will run for n * n, n ^ 4, n ^ 8, etc. i.e. squaring each time instead of adding n-squared each time. Try this:

ulong nSquared = n * n;
for (ulong k = nSquared; k <= limit; k += nSquared)
  isPrime[k] = false;

Answer 2

Here's another implementation. It uses BitArray to save memory. The Parallel.For needs .NET Framework 4.

static List<int> FindPrimesBySieveOfAtkins(int max)
{
    var isPrime = new BitArray((int)max+1, false);
    var sqrt = (int)Math.Sqrt(max);

    Parallel.For(1, sqrt, x =>
    {
        var xx = x * x;
        for (int y = 1; y <= sqrt; y++)
        {
            var yy = y * y;
            var n = 4 * xx + yy;
            if (n <= max && (n % 12 == 1 || n % 12 == 5))
                isPrime[n] = !isPrime[n];

            n = 3 * xx + yy;
            if (n <= max && n % 12 == 7)
                isPrime[n] = !isPrime[n];

            n = 3 * xx - yy;
            if (x > y && n <= max && n % 12 == 11)
                isPrime[n] = !isPrime[n];
        }
    });

    var primes = new List<int>() { 2, 3 };
    for (int n = 5; n <= sqrt; n++)
    {
        if (isPrime[n])
        {
            primes.Add(n);
            int nn = n * n;
            for (int k = nn; k <= max; k += nn)
                isPrime[k] = false;
        }
    }

    for (int n = sqrt + 1; n <= max; n++)
        if (isPrime[n])
            primes.Add(n);

    return primes;
}