为什么我对Atkins筛子的天真实现排除了5？

Question

我编写了一个非常天真的基于维基百科效率低下但清晰的伪码的Atkin筛子实现。我最初是用MATLAB编写算法的，它省略了5作为素数。我还用Python编写了算法，结果也是一样的。

从技术上讲，我知道为什么要排除5；在n = 4*x^2 + y^2的步骤中，当x == 1和y == 1时，n == 5和y==1。这只发生一次，所以5从素数翻转到非素数，从不翻转。

为什么我的算法与维基百科上的算法不匹配？虽然我做了一些肤浅的调整(例如每次迭代只计算x^2一次，在第一个方程中使用它时存储mod(n，12)的值，等等)。他们不应该改变算法的逻辑。

我从几个 讨论 相关中读到了Atkin的筛子，但是我不知道有什么不同在我的实现中造成了问题。

Python代码：

def atkin1(limit):
    res = [0] * (limit + 1)
    res[2] = 1
    res[3] = 1
    res[5] = 1

    limitSqrt = int(math.sqrt(limit))
    for x in range(1, limitSqrt+1):
        for y in range(1, limitSqrt+1):
            x2 = x**2
            y2 = y**2
            n = 4*x2 + y2
            if n == 5:
                print('debug1')
            nMod12 = n % 12
            if n <= limit and (nMod12 == 1 or nMod12 == 5):
                res[n] ^= 1

            n = 3*x2 + y2
            if n == 5:
                print('debug2')
            if n <= limit and (n % 12 == 7):
                res[n] ^= 1

            if x > y:
                n = 3*x2 - y2
                if n == 5:
                    print('debug3')
                if n <= limit and (n % 12 == 11):
                    res[n] ^= 1

    ndx = 5
    while ndx <= limitSqrt:
        m = 1
        if res[ndx]:
            ndx2 = ndx**2
            ndx2Mult =m * ndx2
            while ndx2Mult < limit:
                res[ndx2Mult] = 0
                m += 1
                ndx2Mult = m * ndx2
        ndx += 1

    return res

MATLAB代码

function p = atkin1(limit)

% 1. Create a results list, filled with 2, 3, and 5
res = [0, 1, 1, 0, 1];

% 2. Create a sieve list with an entry for each positive integer; all entries of
% this list should initially be marked nonprime (composite).
res = [res, zeros(1, limit-5)];

% 3. For each entry number n in the sieve list, with modulo-sixty remainder r:

limitSqrt = floor(sqrt(limit));
for x=1:limitSqrt
    for y=1:limitSqrt
        x2 = x^2;       y2 = y^2;

        % If r is 1, 13, 17, 29, 37, 41, 49, or 53, flip the entry for each
        % possible solution to 4x^2 + y^2 = n.
        n = 4*x2 + y2;
        nMod12 = mod(n, 12); 
        if n <= limit && (nMod12 == 1 || nMod12 == 5)
            res(1, n) = ~res(1, n);
        end

        % If r is 7, 19, 31, or 43, flip the entry for each possible solution
        % to 3x^2 + y^2 = n.
        n = 3*x2 + y2;
        if n <= limit && mod(n, 12) == 7
            res(1, n) = ~res(1, n);
        end

        % If r is 11, 23, 47, or 59, flip the entry for each possible solution
        % to 3x^2 - y^2 = n when x > y.
        if x > y
            n = 3*x2 - y2;
            if n <= limit && mod(n, 12) == 11
                res(1, n) = ~res(1, n);
            end
        end

        % If r is something else, ignore it completely.
    end
end

   % 4. Start with the lowest number in the sieve list.
ndx = 5;
while ndx < limitSqrt
    m = 1;
    if res(ndx)
        % 5. Take the next number in the sieve list still marked prime.
        % 6. Include the number in the results list.
        % 7. Square the number and mark all multiples of that square as nonprime.
        ndx2 = ndx^2;
        ndx2Mult = m * ndx2;
        while ndx2Mult < limit
            res(ndx2Mult) = 0;
            m = m + 1;
            ndx2Mult = m * ndx2;
        end
    end

    % 8. Repeat steps five through eight.
    ndx = ndx + 1;
end

p = find(res); % Find the indexes of nonzerogo
end

维基百科伪码

// arbitrary search limit
limit ← 1000000         

// initialize the sieve
is_prime(i) ← false, ∀ i ∈ [5, limit] 

// 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

Answer 1

在维基百科对算法的文字描述中，“结果列表”和“筛子列表”是两件不同的事情。您的Matlab代码只有一个向量用于这两个。但是5的筛子列表中的初始值应该是“非素数”。

Answer 2

在不看代码的情况下，我在您的描述中读到了以下内容：

so 5 is flipped from prime to nonprime and never flipped back

我想这就是问题所在，它应该被初始化为false (因此是非素数)，因为只有当它保持为素数时才会切换。