我写了一个非常天真的Atna Sieve实现,基于
Wikipedia’s inefficient but clear pseudocode.我最初在MATLAB中编写了算法,它省略了5作为素数.我也用Python编写了算法,结果相同.
从技术上讲,我知道为什么要排除5;在n = 4 * x ^ 2 y ^ 2的步骤中,当x == 1且y == 1时n == 5.这仅发生一次,因此5从素数转换为非素数并且从不翻转.
为什么我的算法与维基百科上的算法不匹配？虽然我做了一些表面调整(例如,在每次迭代中只计算一次x ^ 2,在第一个等式中使用时存储mod(n,12)的值等),但它们不应该改变逻辑.算法.
我在several discussions related阅读了Atkin的Sieve,但我不知道在我的实现中产生问题的区别是什么.
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];
% 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);
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);
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
解决方法
版权声明：本文内容由互联网用户自发贡献，该文观点与技术仅代表作者本人。本站仅提供信息存储空间服务，不拥有所有权，不承担相关法律责任。如发现本站有涉嫌侵权/违法违规的内容， 请发送邮件至 dio@foxmail.com 举报，一经查实，本站将立刻删除。
