我基于Wikipedia 低效但清晰的伪代码编写了一个极其幼稚的 Atkin 筛子实现。我最初在 MATLAB 中编写了算法，它省略了 5 作为质数。我还用 Python 编写了算法，结果相同。
从技术上讲，我知道为什么5 被排除在外；在步骤 where n = 4*x^2 + y^2，n == 5 当 x == 1 和 y == 1 时。这只发生一次，所以 5 从素数翻转到非素数并且永远不会翻转回来。
为什么我的算法与 Wikipedia 上的算法不匹配？虽然我做了一些表面上的调整（例如，每次迭代只计算一次 x^2，在第一个方程中使用它时存储 mod(n, 12) 的值等），但它们不应该改变逻辑算法。
我阅读了几个与阿特金筛子 相关的讨论 ，但我不知道哪些差异在我的实现中造成了问题。
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