C#:阿特金筛法的实现

Question

我想知道这里是否有人有他们愿意分享的阿特金筛法的良好实现。

我正在尝试实现它，但无法完全理解它。这是我目前所拥有的。

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

我几乎只是试图“翻译”在Wikipedia 列出的伪代码，但它无法正常工作。所以要么我误解了什么，要么只是做错了什么。或者很可能两者都...

列出前 500 个素数，我将其用作测试，但我的实现在第 40 个(或 41 个？)处失败。

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

您能发现我的错误吗？您是否有可以分享的实现方案，或者两者兼而有之？

---

我使用的确切测试如下所示:

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

这段代码:

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

似乎不是这个伪代码的忠实翻译:

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

您的代码看起来会运行 n * n、n ^ 4、n ^ 8 等，即每次平方而不是每次添加 n 平方。试试这个:

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