Aaron Mugatroyd的最后一个答案提供了Sieve of Atkin（SoA）的Python源代码翻译，但它可以在以下几个方面进行改进，因为它错过了一些重要的优化：
- 他的答案没有使用完整的模60原始Atkin和Bernstein版本的筛子，而是使用了维基百科文章中略微改进的伪代码的变体，因此使用了大约0.36倍的数字筛选范围组合切换/削减操作;下面的我的代码使用相当高效的非页面段伪代码，如我在评论Sieve of Atkin时的回答中所述，它使用数字范围的约0.26倍来减少工作量，使工作量减少了约2/7。
- 他的代码通过仅具有奇数表示来减小缓冲区大小，而我的代码进一步位打包以消除任何可被三和五整除的数字表示以及暗示为“仅奇数”的可被二整除的数字表示;这将内存要求进一步减少了近一半（到8/15），并帮助更好地利用CPU缓存，从而降低平均内存访问时间，进一步增加速度。
- 我的代码使用快速查找表（LUT）pop count技术计算素数的数量，需要几乎没有时间来计算，而他使用的位逐位技术大约需要一秒钟;但是，在此示例代码中，即使这段小时间也被从计时部分中删除。
- 最后，我的代码对位操作进行了优化，以获得最少的内部循环代码。例如，它不使用连续的右移一来生成奇数表示索引，并且实际上几乎不使用位移位通过将所有内部循环编写为常量模数（等于位位置）操作。同样，Aaron的翻译代码在操作方面相当低效，例如在素数平方自由削减中，它将素数的平方添加到索引中，然后检查奇数结果，而不仅仅是两倍的平方，因此不需要进行检查;然后，在内部循环中执行cull操作之前，将数字向右移动一位（除以二），就像在所有循环中一样，使得甚至检查都变得多余。这种低效的代码在使用这种“大筛选缓冲区数组”技术的大范围内执行时间不会产生太大影响，因为每个操作的大部分时间用于RAM内存访问（对于十亿范围而言约为37个CPU时钟周期或更多），但对于适合CPU缓存的较小范围，执行时间将比必要的慢得多;换句话说，它设置了每个操作的太高的最低限制。
代码如下：
using System;
using System.Collections;
using System.Collections.Generic;
using System.Linq;
namespace NonPagedSoA {
class SoA : IEnumerable<ulong> {
private ushort[] buf = null;
private long cnt = 0;
private long opcnt = 0;
private static byte[] modPRMS = { 7, 11, 13, 17, 19, 23, 29, 31, 37, 41, 43, 47, 49, 53, 59, 61 };
private static ushort[] modLUT;
private static byte[] cntLUT;
static SoA() {
modLUT = new ushort[60];
for (int i = 0, m = 0; i < modLUT.Length; ++i) {
if ((i & 1) != 0 || (i + 7) % 3 == 0 || (i + 7) % 5 == 0) modLUT[i] = 0;
else modLUT[i] = (ushort)(1 << (m++));
}
cntLUT = new byte[65536];
for (int i = 0; i < cntLUT.Length; ++i) {
var c = 0;
for (int j = i; j > 0; j >>= 1) c += j & 1;
cntLUT[i] = (byte)c;
}
}
public SoA(ulong range) {
this.opcnt = 0;
if (range < 7) {
if (range > 1) {
cnt = 1;
if (range > 2) this.cnt += (long)(range - 1) / 2;
}
this.buf = new ushort[0];
}
else {
this.cnt = 3;
var nrng = range - 7; var lmtw = nrng / 60;
this.buf = new ushort[lmtw + 1];
ulong n = 6;
for (uint x = 1, y = 3; n <= nrng; n += (x << 3) + 4, ++x, y = 1) {
var cb = n; if (x <= 1) n -= 8;
for (uint i = 0; i < 15 && cb <= range; cb += (y << 2) + 4, y += 2, ++i) {
var cbd = cb / 60; var cm = modLUT[cb % 60];
if (cm != 0)
for (uint c = (uint)cbd, my = y + 15; c < buf.Length; c += my, my += 30) {
buf[c] ^= cm;
}
}
}
n = 0;
for (uint x = 1, y = 2; n <= nrng; n += ((x + x + x) << 2) + 12, x += 2, y = 2) {
var cb = n;
for (var i = 0; i < 15 && cb <= range; cb += (y << 2) + 4, y += 2, ++i) {
var cbd = cb / 60; var cm = modLUT[cb % 60];
if (cm != 0)
for (uint c = (uint)cbd, my = y + 15; c < buf.Length; c += my, my += 30) {
buf[c] ^= cm;
}
}
}
n = 4;
for (uint x = 2, y = x - 1; n <= nrng; n += (ulong)(x << 2) + 4, y = x, ++x) {
var cb = n; int i = 0;
for ( ; y > 1 && i < 15 && cb <= nrng; cb += (ulong)(y << 2) - 4, y -= 2, ++i) {
var cbd = cb / 60; var cm = modLUT[cb % 60];
if (cm != 0) {
uint c = (uint)cbd, my = y;
for ( ; my >= 30 && c < buf.Length; c += my - 15, my -= 30) {
buf[c] ^= cm;
}
if (my > 0 && c < buf.Length) { buf[c] ^= cm; }
}
}
if (y == 1 && i < 15) {
var cbd = cb / 60; var cm = modLUT[cb % 60];
if ((cm & 0x4822) != 0 && cbd < (ulong)buf.Length) { buf[cbd] ^= cm; }
}
}
for (uint i = 0, w = 0, pd = 0, pn = 0, msk = 1; w < this.buf.Length ; ++i) {
uint p = pd + modPRMS[pn];
ulong sqr = (ulong)p * (ulong)p;
if (sqr > range) break;
if ((this.buf[w] & msk) != 0) {
ulong s = sqr - 7;
for (int j = 0; s <= nrng && j < modPRMS.Length; s = sqr * modPRMS[j] - 7, ++j) {
var cd = s / 60; var cm = (ushort)(modLUT[s % 60] ^ 0xFFFF);
for (ulong c = cd; c < (ulong)this.buf.Length; c += sqr) {
this.buf[c] &= cm;
}
}
}
if (msk >= 0x8000) { msk = 1; pn = 0; ++w; pd += 60; }
else { msk <<= 1; ++pn; }
}
var ndx = nrng % 60;
for (; modLUT[ndx] == 0; --ndx) ;
this.buf[lmtw] &= (ushort)((modLUT[ndx] << 1) - 1);
}
}
public long Count {
get {
long cnt = this.cnt;
for (int i = 0; i < this.buf.Length; ++i) cnt += cntLUT[this.buf[i]];
return cnt;
}
}
public long Ops {
get {
return this.opcnt;
}
}
public IEnumerator<ulong> GetEnumerator() {
yield return 2; yield return 3; yield return 5;
ulong pd = 0;
for (uint i = 0, w = 0, pn = 0, msk = 1; w < this.buf.Length; ++i) {
if ((this.buf[w] & msk) != 0)
yield return pd + modPRMS[pn];
if (msk >= 0x8000) { msk = 1; pn = 0; ++w; pd += 60; }
else { msk <<= 1; ++pn; }
}
}
IEnumerator IEnumerable.GetEnumerator() {
return this.GetEnumerator();
}
}
class Program {
static void Main(string[] args) {
Console.WriteLine("This program calculates primes by a simple full version of the Sieve of Atkin.\r\n");
const ulong n = 1000000000;
var elpsd = -DateTime.Now.Ticks;
var gen = new SoA(n);
elpsd += DateTime.Now.Ticks;
Console.WriteLine("{0} primes found to {1} using {2} operations in {3} milliseconds.", gen.Count, n, gen.Ops, elpsd / 10000);
Console.Write("\r\nPress any key to exit:");
Console.ReadKey(true);
Console.WriteLine();
}
}
}
这段代码的速度比Aaron的代码快了大约两倍（在i7-2700K（3.5 GHz）上以64位或32位模式运行，缓冲区大小约为16.5兆字节，并且在筛选范围为十亿时进行了约0.258十亿次的组合开关/质数平方自由的剔除操作（取消对“++this.opcnt”语句的注释可查看），相比之下，如果不计算时间，他的代码需要5.4 / 6.2秒（32位/ 64位），使用了近似0.359十亿次的组合开关/剔除操作，用于筛选最多达到十亿的数），但这并不意味着Atkin筛法比Eratosthenes筛法更快，因为如果将上述SoA实现中使用的技术应用到使用最大轮因子分解的SoE中，SoE的速度将与此相同。
分析：虽然完全优化过的SoE和SoA的操作数量在筛选范围为十亿时大致相同，但是对于这些未分页的实现，主要瓶颈是在筛选缓冲区大小超过CPU缓存大小（我的i7具有32 KiloBytes L1高速缓存、256 Kilobytes L2高速缓存和8 Megabytes L3高速缓存，分别需要1、4和20个时钟周期的时间访问）后的内存访问。此时，内存访问可能超过100个时钟周期。
当将算法适应于分页分段时，两种方法的内存访问速度都有约8倍的改进。然而，随着筛选范围变得非常大，SoE仍然比SoA更优越，因为由于Culling扫描中快速增加到许多倍页面缓冲区大小的概率极小的数质数自由部分的实现难度，所以难以实现该算法。另外，也更严重的是，对于每个x值的新起始点计算与每个页面缓冲区最低表示的y值相比，Paged SoA的效率损失很大，随着筛选范围的增长。
EDIT_ADD: 在Aaron Murgatroyd使用的仅限于奇数的SoE中，对于一亿个筛子范围，使用了大约10.26亿次筛选操作，所以比SoA多进行了四倍的操作，因此应该运行得慢四倍。但是，即使在这里实现的SoA有一个更复杂的内部循环，尤其是由于仅用于奇数的SoE筛除中具有较高比例的步幅比SoA短得多，因此尽管筛子缓冲区远远超过CPU缓存大小（更好地使用缓存关联性），但天真的仅限于奇数的SoE具有更好的平均内存访问时间。这就解释了为什么上面的SoA只比仅限于奇数的SoE快两倍，尽管理论上看起来只执行了四分之一的工作。
如果使用与上述SoA相同的常数取模内部循环的类似算法，并实现相同的2/3/5轮因数分解，则SoE将减少约0.405亿个操作，因此仅比SoA多进行50％的操作，并且理论上略慢于SoA，但可能由于筛除步幅仍然比平均值小一些而以与SoA几乎相同的速度运行，使用这种“天真的”大型内存缓冲区。将轮因数分解增加到2/3/5/7轮意味着对于10亿个筛子范围，SoE的筛选操作减少到约0.314，可能使该版本的SoE在此算法中以相同的速度运行。
可以通过预筛选2/3/5/7/11/13/17/19素数因子的筛子数组（复制模式）来进一步使用轮因数分解，几乎不会增加执行时间成本，从而将总筛选操作数量降至10亿个筛子范围下的大约0.251亿，并且即使对于这些大型内存缓冲区版本，SoE仍然比上述优化版本的SoA具有更少的代码复杂度，速度更快或者与之相当。
因此，可以看出，与天真的或仅限于奇数的或2/3/5轮因数分解版本相比，SoE的操作次数可以大大减少，使操作次数与SoA相同，同时由于较简单的内部循环和更高效的内存访问，每个操作的时间实际上可能更短。END_EDIT_ADD
EDIT_ADD2: 我在这里添加了一个使用类似于上面SoA伪代码中链接答案所示的常数模/bit位置技术的SoE代码。 尽管具有高形式分解和预处理，使得实际被削减的操作总数实际上少于两十亿个素数范围内的组合切换/削减操作，但该代码比上述SoA要简单得多。 以下是该代码：
EDIT_FINAL 下面是已校正的代码及相关注释 END_EDIT_FINAL
using System;
using System.Collections;
using System.Collections.Generic;
namespace NonPagedSoE {
class SoE : IEnumerable<ulong> {
private ushort[] buf = null;
private long cnt = 0;
private long opcnt = 0;
private static byte[] basePRMS = { 2, 3, 5, 7, 11, 13, 17, 19 };
private static byte[] modPRMS = { 23, 29, 31, 37, 41, 43, 47, 53, 59, 61, 67, 71, 73, 79, 83, 89,
97, 101, 103, 107, 109, 113, 121, 127, 131, 137, 139, 143, 149, 151, 157, 163,
167, 169, 173, 179, 181 ,187 ,191 ,193, 197, 199, 209, 211, 221, 223, 227, 229 };
private static byte[] gapsPRMS = { 6, 2, 6, 4, 2, 4, 6, 6, 2, 6, 4, 2, 6, 4, 6, 8,
4, 2, 4, 2, 4, 8, 6, 4, 6, 2, 4, 6, 2, 6, 6, 4,
2, 4, 6, 2, 6, 4, 2, 4, 2, 10, 2, 10, 2, 4, 2, 4 };
private static ulong[] modLUT;
private static byte[] cntLUT;
static SoE() {
modLUT = new ulong[210];
for (int i = 0, m = 0; i < modLUT.Length; ++i) {
if ((i & 1) != 0 || (i + 23) % 3 == 0 || (i + 23) % 5 == 0 || (i + 23) % 7 == 0) modLUT[i] = 0;
else modLUT[i] = 1UL << (m++);
}
cntLUT = new byte[65536];
for (int i = 0; i < cntLUT.Length; ++i) {
var c = 0;
for (int j = i ^ 0xFFFF; j > 0; j >>= 1) c += j & 1;
cntLUT[i] = (byte)c;
}
}
public SoE(ulong range) {
this.opcnt = 0;
if (range < 23) {
if (range > 1) {
for (int i = 0; i < modPRMS.Length; ++i) if (modPRMS[i] <= range) this.cnt++; else break;
}
this.buf = new ushort[0];
}
else {
this.cnt = 8;
var nrng = range - 23; var lmtw = nrng / 210; var lmtwt3 = lmtw * 3;
this.buf = new ushort[lmtwt3 + 3];
for (int pn = modPRMS.Length - 4; pn < modPRMS.Length; ++pn) {
uint p = modPRMS[pn] - 210u; ulong pt3 = p * 3;
ulong s = p * p - 23;
ulong xrng = Math.Min(9699709, nrng);
ulong nwrds = (ulong)Math.Min(138567, this.buf.Length);
for (int j = 0; s <= xrng && j < modPRMS.Length; s += p * gapsPRMS[(pn + j++) % 48]) {
var sm = modLUT[s % 210];
var si = (sm < (1UL << 16)) ? 0UL : ((sm < (1UL << 32)) ? 1UL : 2UL);
var cd = s / 210 * 3 + si; var cm = (ushort)(sm >> (int)(si << 4));
for (ulong c = cd; c < nwrds; c += pt3) {
this.buf[c] |= cm;
}
}
}
for (long i = 138567; i < this.buf.Length; i += 138567)
if (i + 138567 <= this.buf.Length)
Array.Copy(this.buf, 0, this.buf, i, 138567);
else Array.Copy(this.buf, 0, this.buf, i, this.buf.Length - i);
for (uint i = 0, w = 0, wi = 0, pd = 0, pn = 0, msk = 1; w < this.buf.Length; ++i) {
uint p = pd + modPRMS[pn];
ulong sqr = (ulong)p * (ulong)p;
if (sqr > range) break;
if ((this.buf[w] & msk) == 0) {
ulong s = sqr - 23; ulong pt3 = p * 3;
for (int j = 0; s <= nrng && j < modPRMS.Length; s += p * gapsPRMS[(pn + j++) % 48]) {
var sm = modLUT[s % 210];
var si = (sm < (1UL << 16)) ? 0UL : ((sm < (1UL << 32)) ? 1UL : 2UL);
var cd = s / 210 * 3 + si; var cm = (ushort)(sm >> (int)(si << 4));
for (ulong c = cd; c < (ulong)this.buf.Length; c += pt3) {
this.buf[c] |= cm;
}
}
}
++pn;
if (msk >= 0x8000) { msk = 1; ++w; ++wi; if (wi == 3) { wi = 0; pn = 0; pd += 210; } }
else msk <<= 1;
}
var ndx = nrng % 210;
for (; modLUT[ndx] == 0; --ndx) ;
var cmsk = (~(modLUT[ndx] - 1)) << 1;
this.buf[lmtwt3++] |= (ushort)cmsk;
this.buf[lmtwt3++] |= (ushort)(cmsk >> 16);
this.buf[lmtwt3] |= (ushort)(cmsk >> 32);
}
}
public long Count {
get {
long cnt = this.cnt;
for (int i = 0; i < this.buf.Length; ++i) cnt += cntLUT[this.buf[i]];
return cnt;
}
}
public long Ops {
get {
return this.opcnt;
}
}
public IEnumerator<ulong> GetEnumerator() {
yield return 2; yield return 3; yield return 5; yield return 7;
yield return 11; yield return 13; yield return 17; yield return 19;
ulong pd = 0;
for (uint i = 0, w = 0, wi = 0, pn = 0, msk = 1; w < this.buf.Length; ++i) {
if ((this.buf[w] & msk) == 0)
yield return pd + modPRMS[pn];
++pn;
if (msk >= 0x8000) { msk = 1; ++w; ++wi; if (wi == 3) { wi = 0; pn = 0; pd += 210; } }
else msk <<= 1;
}
}
IEnumerator IEnumerable.GetEnumerator() {
return this.GetEnumerator();
}
}
class Program {
static void Main(string[] args) {
Console.WriteLine("This program calculates primes by a simple maximually wheel factorized version of the Sieve of Eratosthenes.\r\n");
const ulong n = 1000000000;
var elpsd = -DateTime.Now.Ticks;
var gen = new SoE(n);
elpsd += DateTime.Now.Ticks;
Console.WriteLine("{0} primes found to {1} using {2} operations in {3} milliseconds.", gen.Count, n, gen.Ops, elpsd / 10000);
Console.Write("\r\nPress any key to exit:");
Console.ReadKey(true);
Console.WriteLine();
}
}
}
这段代码实际上比上面的结构体数组运行速度快几个百分点，因为操作略少，对于十亿级别的大数组来说，主要瓶颈是内存访问时间，约为40到100个CPU时钟周期，具体取决于CPU和内存规格。这意味着代码优化（除了减少总操作数之外）是无效的，因为大部分时间都花在等待内存访问上。总之，使用巨大的内存缓冲区不是筛选大范围的最有效方法，使用页面分割可以将SoE的因子提高多达8倍，同时使用相同的最大轮因子分解（也为多处理铺平道路）。
在实现页面分割和多处理方面，与SoE相比，当范围远高于40亿时，SoA确实存在缺陷，因为由于素数平方自由处理和计算更多的页面起始地址相关的页面处理开销因素，SoA的渐近复杂度降低所带来的任何收益都会迅速消耗殆尽；或者，通过以巨大的内存消耗为代价，在RAM内存中存储标记，并进一步降低访问这些标记存储结构的效率来克服这一问题。
简而言之，与完全轮因子分解的SoE相比，SoA并不是一个实用的筛法，因为随着渐近复杂度的提高，它的性能越来越接近于完全优化的SoE，但由于实际实现的细节，如相对内存访问时间和页面分割复杂性以及一般更为复杂和难写，它开始失去效率。在我看来，它更像是一个有趣的智力概念和思维锻炼，而不是实用的筛法。
总有一天，我会将这些技术适应于多线程页面分割的埃拉托色尼筛法，使其在C#中的速度与Atkin和Bernstein的“primegen” SoA实现相当，并在超过40亿的大范围内甚至单线程时也能获得额外的速度提升，而在我的i7上进行多线程处理时可获得多达4倍的速度提升（包括超线程的八个核心）。
