这是 C# 中无界素数筛选的解决方案，可以使用埃拉托色尼筛 (SoE) 或阿特金筛 (SoA) 算法实现；然而，我坚持认为，与真正的 SoE 提供大致相同的性能而没有那么多复杂性相比，给定的优化 SoA 解决方案的极端复杂性是不值得的。因此，这可能只是部分答案，虽然它展示了如何实现更好的 SoA 算法并展示了如何使用 SoE 实现无界序列，但它仅暗示了如何将这些组合起来以编写一个相当有效的 SoA。
请注意，如果只需要讨论这些想法的最快实现，请跳至此答案的底部。
首先，我们应该评论这个练习的要点，即生成一个无界的素数序列，以允许使用 IEnumerable 扩展方法，例如 Take()、TakeWhile()、Where()、Count() 等，因为这些方法放弃了一些由于增加了方法调用的级别，每次调用至少增加了 28 个机器周期并返回以枚举一个值，并为每个函数添加了多个级别的方法调用；也就是说，即使对枚举结果使用更多命令式过滤技术以提高速度，拥有（有效）无限序列仍然有用。
请注意，在以下更简单的示例中，我将范围限制为无符号 32 位数字 (uint)用于部分筛分以提高效率的“桶”或其他形式的增量筛；然而，对于在最快的实现中完全优化的页面分段筛选器，范围增加到 64 位范围，尽管在编写时可能不希望筛选超过大约 100 万亿（10 次方）作为运行时间整个范围可能需要数百年的时间。
由于问题选择 SoA 的原因可能是错误的，首先将 Trial Division (TD) 素筛误认为是真正的 SoE，然后在 TD 筛上使用朴素 SoA，让我们建立真正的有界 SoE，它可以在几个行作为方法（可以根据问题的实现编码风格转换为类）如下：
static IEnumerable<uint> primesSoE(uint top_number) {
if (top_number < 2u) yield break;
yield return 2u; if (top_number < 3u) yield break;
var BFLMT = (top_number - 3u) / 2u;
var SQRTLMT = ((uint)(Math.Sqrt((double)top_number)) - 3u) / 2u;
var buf = new BitArray((int)BFLMT + 1,true);
for (var i = 0u; i <= BFLMT; ++i) if (buf[(int)i]) {
var p = 3u + i + i; if (i <= SQRTLMT) {
for (var j = (p * p - 3u) / 2u; j <= BFLMT; j += p)
buf[(int)j] = false; } yield return p; } }
这可以在 i7-2700K (3.5 GHz) 上在大约 16 毫秒内计算出 200 万个素数，在同一台机器上大约 67 秒内将 203,280,221 个素数在 32 位数字范围内计算到 4,294,967,291（假设备用 256 MB RAM记忆）。
现在，使用上面的版本与 SoA 进行比较是不公平的，因为真正的 SoA 会自动忽略 2、3 和 5 的倍数，因此引入车轮分解对 SoE 执行相同操作会产生以下代码：
static IEnumerable<uint> primesSoE(uint top_number) {
if (top_number < 2u) yield break;
yield return 2u; if (top_number < 3u) yield break;
yield return 3u; if (top_number < 5u) yield break;
yield return 5u; if (top_number < 7u) yield break;
var BFLMT = (top_number - 7u) / 2u;
var SQRTLMT = ((uint)(Math.Sqrt((double)top_number)) - 7u) / 2u;
var buf = new BitArray((int)BFLMT + 1,true);
byte[] WHLPTRN = { 2, 1, 2, 1, 2, 3, 1, 3 };
for (uint i = 0u, w = 0u; i <= BFLMT; i += WHLPTRN[w], w = (w < 7u) ? ++w : 0u)
if (buf[(int)i]) { var p = 7u + i + i; if (i <= SQRTLMT) {
var pX2 = p + p; uint[] pa = { p, pX2, pX2 + p };
for (uint j = (p * p - 7u) / 2u, m = w; j <= BFLMT;
j += pa[WHLPTRN[m] - 1u], m = (m < 7u) ? ++m : 0u)
buf[(int)j] = false; } yield return p; } }
上面的代码在大约 10 毫秒内计算出 200 万个素数，在与上述相同的机器上大约 40 秒内计算出 32 位数字范围内的素数。
接下来，让我们确定我们可能在此处编写的 SoA 版本与根据上述代码的 SoE 相比，就执行速度而言，实际上是否具有任何优势。下面的代码遵循上面SoE的模型，优化了维基百科文章中的SoA伪代码至于使用单独的循环来调整“x”和“y”变量的范围，如那篇文章中所建议的那样，只为奇数解求解二次方程（和无平方消除），结合“3* x^2" 二次以在同一遍中求解正负第二项，并消除计算上昂贵的模运算，生成的代码比那里发布的伪代码的原始版本快三倍以上，并且在这里的问题中使用。有界 SoA 的代码如下：
static IEnumerable<uint> primesSoA(uint top_number) {
if (top_number < 2u) yield break;
yield return 2u; if (top_number < 3u) yield break;
yield return 3u; if (top_number < 5u) yield break;
var BFLMT = (top_number - 3u) / 2u; var buf = new BitArray((int)BFLMT + 1, false);
var SQRT = (uint)(Math.Sqrt((double)top_number)); var SQRTLMT = (SQRT - 3u) / 2u;
for (uint x = 1u, s = 1u, mdx12 = 5u, dmdx12 = 0u; s <= BFLMT; ++x, s += ((x << 1) - 1u) << 1) {
for (uint y = 1u, n = s, md12 = mdx12, dmd12 = 8u; n <= BFLMT; y += 2, n += (y - 1u) << 1) {
if ((md12 == 1) || (md12 == 5)) buf[(int)n] = buf[(int)n] ^ true;
md12 += dmd12; if (md12 >= 12) md12 -= 12; dmd12 += 8u; if (dmd12 >= 12u) dmd12 -= 12u; }
mdx12 += dmdx12; if (mdx12 >= 12u) mdx12 -= 12u; dmdx12 += 8u; if (dmdx12 >= 12u) dmdx12 -= 12u; }
for (uint x = 1u, s = 0u, mdx12 = 3u, dmdx12 = 8u; s <= BFLMT; ++x, s += x << 1) {
int y = 1 - (int)x, n = (int)s, md12 = (int)mdx12, dmd12 = ((-y - 1) << 2) % 12;
for (; (y < 0) && (uint)n <= BFLMT; y += 2, n += (-y + 1) << 1) {
if (md12 == 11) buf[(int)n] = buf[(int)n] ^ true;
md12 += dmd12; if (md12 >= 12) md12 -= 12; dmd12 += 4; if (dmd12 >= 12) dmd12 -= 12; }
if (y < 1) { y = 2; n += 2; md12 += 4; dmd12 = 0; } else { n += 1; md12 += 2; dmd12 = 8; }
if (md12 >= 12) md12 -= 12; for (; (uint)n <= BFLMT; y += 2, n += (y - 1) << 1) {
if (md12 == 7) buf[(int)n] = buf[(int)n] ^ true;
md12 += dmd12; if (md12 >= 12) md12 -= 12; dmd12 += 8; if (dmd12 >= 12) dmd12 -= 12; }
mdx12 += dmdx12; if (mdx12 >= 12) mdx12 -= 12; dmdx12 += 4; if (dmdx12 >= 12) dmdx12 -= 12; }
for (var i = 0u; i<=BFLMT; ++i) if (buf[(int)i]) { var p = 3u+i+i; if (i<=SQRTLMT) { var sqr = p*p;
for (var j = (sqr - 3ul) / 2ul; j <= BFLMT; j += sqr) buf[(int)j] = false; } yield return p; } }
由于未完全优化的操作数量，这仍然是发布的车轮分解 SoE 算法的两倍多。可以对 SoA 代码进行进一步优化，如对原始（非伪代码）算法使用模 60 运算，并使用位打包仅处理不包括 3 和 5 的倍数，以使代码在性能上与 SoE 相当接近甚至稍微超过它，但我们越来越接近 Berstein 实现的复杂性以实现这种性能。我的观点是“为什么是 SoA，当一个人非常努力地获得与 SoE 相同的性能时？”。
现在对于无界素数序列，最简单的方法是定义 Uint32.MaxValue 的 const top_number 并消除 primesSoE 或 primesSoA 方法中的参数。这有点低效，因为无论正在处理的实际素数值如何，剔除仍然在整个数字范围内完成，这使得确定小范围素数的时间比它应该花费的时间长得多。除了这个和极端内存使用之外，还有其他原因需要使用分段版本的素数筛：首先，使用主要在 CPU L1 或 L2 数据缓存中的数据的算法处理速度更快，因为内存访问效率更高，
为了提高效率，我们应该选择 CPU L1 或 L2 缓存大小的数组大小，对于大多数现代 CPU 来说，它至少为 16 KB，以避免缓存抖动，因为我们从数组中剔除素数的组合，这意味着 BitArray 可以有一个8 倍的大小（每字节 8 位）或 128 千位。由于我们只需要筛选奇数素数，这代表了超过 25 万的数字范围，这意味着分段版本将仅使用 8 段来筛选欧拉问题 10 所要求的 200 万。人们可以将结果从最后一段并从那一点继续，但这会排除将此代码调整为多核情况，在这种情况下，将剔除多个线程的后台以充分利用多核处理器。
static IEnumerable<uint> primesSoE() { yield return 2u; yield return 3u; yield return 5u;
const uint L1CACHEPOW = 14u + 3u, L1CACHESZ = (1u << (int)L1CACHEPOW); //for 16K in bits...
const uint BUFSZ = L1CACHESZ / 15u * 15u; //an even number of wheel rotations
var buf = new BitArray((int)BUFSZ);
const uint MAXNDX = (uint.MaxValue - 7u) / 2u; //need maximum for number range
var SQRTNDX = ((uint)Math.Sqrt(uint.MaxValue) - 7u) / 2u;
byte[] WHLPTRN = { 2, 1, 2, 1, 2, 3, 1, 3 }; //the 2,3,5 factorial wheel, (sum) 15 elements long
byte[] WHLPOS = { 0, 2, 3, 5, 6, 8, 11, 12 }; //get wheel position from index
byte[] WHLNDX = { 0, 0, 1, 2, 2, 3, 4, 4, 5, 5, 5, 6, 7, 7, 7, //get index from position
0, 0, 1, 2, 2, 3, 4, 4, 5, 5, 5, 6, 7 }; //allow for overflow
byte[] WHLRNDUP = { 0, 2, 2, 3, 5, 5, 6, 8, 8, 11, 11, 11, 12, 15, //allow for overflow...
15, 15, 17, 17, 18, 20, 20, 21, 23, 23, 26, 26, 26, 27 };
uint BPLMT = (ushort.MaxValue - 7u) / 2u; var bpbuf = new BitArray((int)BPLMT + 1, true);
for (var i = 0; i <= 124; ++i) if (bpbuf[i]) { var p = 7 + i + i; //initialize baseprimes array
for (var j = (p * p - 7) / 2; j <= BPLMT; j += p) bpbuf[j] = false; } var pa = new uint[3];
for (uint i = 0u, w = 0, si = 0; i <= MAXNDX;
i += WHLPTRN[w], si += WHLPTRN[w], si = (si >= BUFSZ) ? 0u : si, w = (w < 7u) ? ++w : 0u) {
if (si == 0) { buf.SetAll(true);
for (uint j = 0u, bw = 0u; j <= BPLMT; j += WHLPTRN[bw], bw = (bw < 7u) ? ++bw : 0u)
if (bpbuf[(int)j]) { var p = 7u+j+j; var pX2=p+p; var k = p * (j + 3u) + j;
if (k >= i + BUFSZ) break; pa[0] = p; pa[1] = pX2; pa[2] = pX2 + p; var sw = bw; if (k < i) {
k = (i - k) % (15u * p); if (k != 0) { var os = WHLPOS[bw]; sw = os + ((k + p - 1u) / p);
sw = WHLRNDUP[sw]; k = (sw - os) * p - k; sw = WHLNDX[sw]; } } else k -= i;
for (; k<BUFSZ; k+=pa[WHLPTRN[sw]-1], sw=(sw<7u) ? ++sw : 0u) buf[(int)k]=false; } }
if (buf[(int)si]) yield return 7u + i + i; } }
上面的代码大约需要 16 毫秒来筛选最多 200 万个素数，大约需要 30 秒来筛选完整的 32 位数字范围。此代码比使用阶乘轮的类似版本更快，因为即使代码在计算方面更复杂，但在不破坏 CPU 缓存的情况下节省了大量时间。大部分复杂性在于计算每个段的每个基本素数的新起始索引，这可以通过保存每个段的每个素数的状态来减少，但是上面的代码可以很容易地转换，以便在对于具有四核的机器，单独的后台线程进一步加快了四倍的速度，甚至更多地具有八核。不使用 BitArray 类并通过位掩码寻址各个位位置将提供大约两倍的加速，并且阶乘轮的进一步增加将提供大约三分之二的时间减少。将位数组更好地包装在由轮因式分解消除的值的消除索引中将略微提高大范围的效率，但也会在一定程度上增加位操作的复杂性。然而，所有这些 SoE 优化都不会接近 Berstein SoA 实现的极端复杂性，但运行速度大致相同。将位数组更好地包装在由轮因式分解消除的值的消除索引中将略微提高大范围的效率，但也会在一定程度上增加位操作的复杂性。然而，所有这些 SoE 优化都不会接近 Berstein SoA 实现的极端复杂性，但运行速度大致相同。将位数组更好地包装在由轮因式分解消除的值的消除索引中将略微提高大范围的效率，但也会在一定程度上增加位操作的复杂性。然而，所有这些 SoE 优化都不会接近 Berstein SoA 实现的极端复杂性，但运行速度大致相同。
要将上面的代码从 SoE 转换为 SoA，我们只需要从有界情况更改为 SoA 剔除代码，但要修改为每个页段重新计算起始地址，就像为 SoE 情况计算起始地址一样，但在操作上甚至更复杂一些，因为二次方程是通过数字的平方而不是简单的素数倍数来推进的。我将把所需的调整留给学生，尽管考虑到经过合理优化的 SoA 并不比当前版本的 SoE 快并且要复杂得多，我并没有真正理解这一点；）
编辑_添加：
注意：下面的代码已经更正，虽然原始发布的代码对于提供的素数序列是正确的，但它比它需要的速度慢了一半，因为它剔除了范围平方根以下的所有数字，而不仅仅是找到的基础素数最多为范围的平方根。
最有效和最简单的优化是将每个段页面的剔除操作委托给后台线程，以便在有足够 CPU 内核的情况下，枚举素数的主要限制是枚举循环本身，在这种情况下，所有素数都是 32 位数range 可以在上述参考机器（在 C# 中）上在大约 10 秒内枚举，无需其他优化，所有其他优化包括编写 IEnumerable 接口实现而不是使用内置迭代器，将所有 203,280,221 个素数减少到大约 6 秒在 32 位数字范围内（1.65 秒到 10 亿），同样大部分时间都花在列举结果上。以下代码进行了许多修改，包括使用 Tasks 使用的 Dotnet Framework 4 ThreadPool 中的后台任务，使用实现为查找表的状态机来实现进一步的位打包以使素数枚举更快，以及重写算法作为一个使用“roll-your-own”迭代器来提高效率的可枚举类：
class fastprimesSoE : IEnumerable<uint>, IEnumerable {
struct procspc { public Task tsk; public uint[] buf; }
struct wst { public byte msk; public byte mlt; public byte xtr; public byte nxt; }
static readonly uint NUMPROCS = (uint)Environment.ProcessorCount + 1u; const uint CHNKSZ = 1u;
const uint L1CACHEPOW = 14u, L1CACHESZ = (1u << (int)L1CACHEPOW), PGSZ = L1CACHESZ >> 2; //for 16K in bytes...
const uint BUFSZ = CHNKSZ * PGSZ; //number of uints even number of caches in chunk
const uint BUFSZBTS = 15u * BUFSZ << 2; //even in wheel rotations and uints (and chunks)
static readonly byte[] WHLPTRN = { 2, 1, 2, 1, 2, 3, 1, 3 }; //the 2,3,5 factorial wheel, (sum) 15 elements long
static readonly byte[] WHLPOS = { 0, 2, 3, 5, 6, 8, 11, 12 }; //get wheel position from index
static readonly byte[] WHLNDX = { 0, 1, 1, 2, 3, 3, 4, 5, 5, 6, 6, 6, 7, 0, 0, 0 }; //get index from position
static readonly byte[] WHLRNDUP = { 0, 2, 2, 3, 5, 5, 6, 8, 8, 11, 11, 11, 12, 15, 15, 15, //allow for overflow...
17, 17, 18, 20, 20, 21, 23, 23, 26, 26, 26, 27, 30, 30, 30 }; //round multiplier up
const uint BPLMT = (ushort.MaxValue - 7u) / 2u; const uint BPSZ = BPLMT / 60u + 1u;
static readonly uint[] bpbuf = new uint[BPSZ]; static readonly wst[] WHLST = new wst[64];
static void cullpg(uint i, uint[] b, int strt, int cnt) { Array.Clear(b, strt, cnt);
for (uint j = 0u, wp = 0, bw = 0x31321212u, bi = 0u, v = 0xc0881000u, bm = 1u; j <= BPLMT;
j += bw & 0xF, wp = wp < 12 ? wp += bw & 0xF : 0, bw = (bw > 3u) ? bw >>= 4 : 0x31321212u) {
var p = 7u + j + j; var k = p * (j + 3u) + j; if (k >= (i + (uint)cnt * 60u)) break;
if ((v & bm) == 0u) { if (k < i) { k = (i - k) % (15u * p); if (k != 0) {
var sw = wp + ((k + p - 1u) / p); sw = WHLRNDUP[sw]; k = (sw - wp) * p - k; } }
else k -= i; var pd = p / 15;
for (uint l = k / 15u + (uint)strt * 4u, lw = ((uint)WHLNDX[wp] << 3) + WHLNDX[k % 15u];
l < (uint)(strt + cnt) * 4u; ) { var st = WHLST[lw];
b[l >> 2] |= (uint)st.msk << (int)((l & 3) << 3); l += st.mlt * pd + st.xtr; lw = st.nxt; } }
if ((bm <<= 1) == 0u) { v = bpbuf[++bi]; bm = 1u; } } }
static fastprimesSoE() {
for (var x = 0; x < 8; ++x) { var p = 7 + 2 * WHLPOS[x]; var pr = p % 15;
for (int y = 0, i = ((p * p - 7) / 2); y < 8; ++y) { var m = WHLPTRN[(x + y) % 8]; i %= 15;
var n = WHLNDX[i]; i += m * pr; WHLST[x * 8 + n] = new wst { msk = (byte)(1 << n), mlt = m,
xtr = (byte)(i / 15),
nxt = (byte)(8 * x + WHLNDX[i % 15]) }; }
} cullpg(0u, bpbuf, 0, bpbuf.Length); } //init baseprimes
class nmrtr : IEnumerator<uint>, IEnumerator, IDisposable {
procspc[] ps = new procspc[NUMPROCS]; uint[] buf;
Task dlycullpg(uint i, uint[] buf) { return Task.Factory.StartNew(() => {
for (var c = 0u; c < CHNKSZ; ++c) cullpg(i + c * PGSZ * 60, buf, (int)(c * PGSZ), (int)PGSZ); }); }
public nmrtr() {
for (var i = 0u; i < NUMPROCS; ++i) ps[i] = new procspc { buf = new uint[BUFSZ] };
for (var i = 1u; i < NUMPROCS; ++i) { ps[i].tsk = dlycullpg((i - 1u) * BUFSZBTS, ps[i].buf); } buf = ps[0].buf; }
public uint Current { get { return this._curr; } } object IEnumerator.Current { get { return this._curr; } }
uint _curr; int b = -4; uint i = 0, w = 0; uint v, msk = 0;
public bool MoveNext() {
if (b < 0) if (b == -1) { _curr = 7; b += (int)BUFSZ; }
else { if (b++ == -4) this._curr = 2u; else this._curr = 7u + ((uint)b << 1); return true; }
do { i += w & 0xF; if ((w >>= 4) == 0) w = 0x31321212u; if ((this.msk <<= 1) == 0) {
if (++b >= BUFSZ) { b = 0; for (var prc = 0; prc < NUMPROCS - 1; ++prc) ps[prc] = ps[prc + 1];
ps[NUMPROCS - 1u].buf = buf; var low = i + (NUMPROCS - 1u) * BUFSZBTS;
ps[NUMPROCS - 1u].tsk = dlycullpg(i + (NUMPROCS - 1u) * BUFSZBTS, buf);
ps[0].tsk.Wait(); buf = ps[0].buf; } v = buf[b]; this.msk = 1; } }
while ((v & msk) != 0u); if (_curr > (_curr = 7u + i + i)) return false; else return true; }
public void Reset() { throw new Exception("Primes enumeration reset not implemented!!!"); }
public void Dispose() { }
}
public IEnumerator<uint> GetEnumerator() { return new nmrtr(); }
IEnumerator IEnumerable.GetEnumerator() { return new nmrtr(); } }
请注意，由于设置和初始化数组的开销，此代码对于小范围的素数（最多 1 或 200 万）并不比上一个版本快，但对于高达 40 亿以上的更大范围则要快得多。它比问题的阿特金筛子的实施快了大约 8 倍，但对于高达 40 亿以上的更大范围，它的速度要快 20 到 50 倍。在设置基本缓存大小的代码中定义了常量，以及每个线程（CHNKSZ）要剔除多少个常量，这可能会稍微调整性能。可以尝试一些轻微的进一步优化，这可能会将大素数的执行时间减少多达两倍，例如对 2,3,5,7 轮进行进一步的位打包，而不仅仅是 2,3，
无论如何，如果一个人打算使用 IEnumberable 接口作为问题所需的输出，则几乎不值得进一步优化，因为大约三分之二的时间仅用于枚举找到的素数，而只有大约三分之一的时间用于剔除合数。更好的方法是在类上编写方法来实现所需的结果，例如 CountTo、ElementAt、SumTo 等，从而完全避免枚举。
但是我确实做了额外的优化作为一个额外的答案，尽管有额外的复杂性，但它显示了额外的好处，并且进一步说明了为什么一个人不想使用 SoA，因为完全优化的 SoE 实际上更好。