这里提供了一种C＃中无界素数筛选的解决方案，可使用埃拉托色尼筛法（SoE）或阿特金筛法（SoA）算法实现。然而，我认为优化的SoA解决方案极其复杂，而真正的SoE在不增加太多复杂性的情况下提供了相同的性能。因此，这可能只是一个部分答案，它展示了如何实现更好的SoA算法并使用SoE实现无界序列，但只是暗示了如何组合这些内容以编写一个相当有效的SoA。
请注意，如果只讨论这些想法最快的实现，则可以跳转到本答案底部。
首先，我们应该评论一下这个练习的要点，即产生一个无界质数序列，以便使用IEnumerable扩展方法（例如Take()，TakeWhile()，Where()，Count()等），因为由于添加了方法调用级别，这些方法会减少一些性能，每次调用和返回枚举一个值至少增加28个机器周期，并添加几个方法调用的级别；尽管如此，即使在对枚举结果使用更好的命令式过滤技术以获得更好的速度时，拥有（有效）无限序列仍然很有用。
请注意，在以下更简单的示例中，我将范围限制为无符号32位数字（uint），因为超过该范围基本的SoE或SoA实现并不真正适用于效率，需要切换到“桶”或其他形式的增量筛法以提高效率。但是，对于完全优化的分段筛法（最快的实现），范围增加到64位范围，尽管如写成一篇文章，人们可能不希望筛选超过约100万亿（10的14次方），因为完整范围的运行时间需要数百年。
由于问题选择SoA的原因可能是错误的，首先将试除法（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秒内计算出32位数字范围内4294967291以下的203,280,221个质数（需要256兆字节的内存空间）。现在，为了与真正的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的模型，并根据
维基百科文章中建议使用单独的循环来调整“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的顶部数字常量，并消除primesSoE或primesSoA方法中的参数。 这在某种程度上是低效的，因为即使处理的实际素数值很小，筛选仍然在整个数字范围内完成，这使得确定小范围的质数需要比应有的时间更长。 还有其他理由采用分段版本的质数筛子，而不是这种极端的内存使用：首先，使用主要在CPU L1或L2数据缓存中的数据的算法由于更有效的内存访问而处理速度更快，其次，分段允许将工作轻松地分割成可以使用多核处理器后台攻略的页面，从而可以获得与使用核心数成比例的加速效果。
为了提高效率，我们应该选择一个CPU L1或L2缓存大小为至少16千字节的数组大小，以避免在筛选合数时发生缓存抖动，这意味着BitArray的大小可以是这个值的8倍（每字节8位），即128千位。 由于我们只需要筛选奇素数，因此这表示数字范围超过25万，这意味着分段版本仅使用八个段来筛选到Euler Problem 10所需的200万个数字。 可以保存最后一段的结果并继续进行，但这将排除将此代码适应多核情况，在多个线程上将整理工作委托给后台以充分利用多核处理器的可能性。 下面是（单线程）无限SoE的C＃代码：
static IEnumerable<uint> primesSoE() { yield return 2u; yield return 3u; yield return 5u;
const uint L1CACHEPOW = 14u + 3u, L1CACHESZ = (1u << (int)L1CACHEPOW);
const uint BUFSZ = L1CACHESZ / 15u * 15u;
var buf = new BitArray((int)BUFSZ);
const uint MAXNDX = (uint.MaxValue - 7u) / 2u;
var SQRTNDX = ((uint)Math.Sqrt(uint.MaxValue) - 7u) / 2u;
byte[] WHLPTRN = { 2, 1, 2, 1, 2, 3, 1, 3 };
byte[] WHLPOS = { 0, 2, 3, 5, 6, 8, 11, 12 };
byte[] WHLNDX = { 0, 0, 1, 2, 2, 3, 4, 4, 5, 5, 5, 6, 7, 7, 7,
0, 0, 1, 2, 2, 3, 4, 4, 5, 5, 5, 6, 7 };
byte[] WHLRNDUP = { 0, 2, 2, 3, 5, 5, 6, 8, 8, 11, 11, 11, 12, 15,
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;
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缓存时间时可以节省很多时间。许多复杂性在于每个基础素数每个段的新起始索引的计算，可以通过保存每个素数每个段的状态来减少这些复杂性，但以上代码可以轻松转换为在单独的后台线程上运行筛选过程以进一步提高4倍速度的机器(四核)。如果不使用BitArray类并通过位掩码寻址单个位位置将提供另外一倍左右的速度提升，并且进一步增加阶乘轮会使时间再减少约三分之二。更好地打包用阶乘轮分解消除的值的位数组将稍微提高大范围的效率，但也会在位操作的复杂性上增加一些。然而，所有这些SoE优化都无法达到Berstein SoA实现的极端复杂性，但其运行速度大约相同。
要将以上代码从SoE转换为SoA，我们只需要从有界情况转换到SoA筛选代码，但修改的是每个页面段的起始地址计算方式，类似于SoE情况下计算起始地址，但操作上更加复杂，因为二次方程的提前是通过数字的平方而不是通过简单的素数倍数进行的。我将把所需的适应性留给学生，尽管我真的看不到这个意义，因为经过合理优化的SoA并不比当前版本的SoE更快，并且复杂性要大得多;)。
编辑添加：
注意：以下代码已经被纠正，因为原始发布的代码虽然对于提供的质数序列是正确的，但它需要比必须剔除所有小于范围平方根的数字而只剔除到范围平方根的找到的基本质数慢约三分之一。
在IT技术方面，最有效和简单的优化措施是将每个段页的筛选操作转移到后台线程中。如果有足够的CPU核心，则在这种情况下，枚举素数的主要限制是枚举循环本身，此时，在上述参考机器上（使用C#），可以在约十秒钟内枚举出32位数字范围内的所有素数，而不需要其他优化措施。包括实现IEnumerable接口而不是使用内置迭代器的所有其他优化都可以将其减少到“32位数字范围内的203,280,221个素数大约需要六秒”（10亿个素数则需1.65秒），同样，大部分时间都花在枚举结果上。以下代码包含了许多这些修改，包括使用从Tasks使用的Dotnet Framework 4线程池中的后台任务，使用作为查找表实现的状态机来实现进一步的位打包以使素数枚举更快，并将算法重写为可枚举类，使用“自己编写”的迭代器以增加效率：
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;
const uint BUFSZ = CHNKSZ * PGSZ;
const uint BUFSZBTS = 15u * BUFSZ << 2;
static readonly byte[] WHLPTRN = { 2, 1, 2, 1, 2, 3, 1, 3 };
static readonly byte[] WHLPOS = { 0, 2, 3, 5, 6, 8, 11, 12 };
static readonly byte[] WHLNDX = { 0, 1, 1, 2, 3, 3, 4, 5, 5, 6, 6, 6, 7, 0, 0, 0 };
static readonly byte[] WHLRNDUP = { 0, 2, 2, 3, 5, 5, 6, 8, 8, 11, 11, 11, 12, 15, 15, 15,
17, 17, 18, 20, 20, 21, 23, 23, 26, 26, 26, 27, 30, 30, 30 };
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); }
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(); } }
请注意，对于小范围的质数（例如到一两百万），此代码并不比上一个版本更快，因为需要设置和初始化数组，但是对于大范围（四十亿及以上）来说，它要快得多。它比Atkin筛法的实现快约8倍，但对于大范围（四十亿及以上）来说，速度可提高20至50倍。代码中定义了常量，设置了基本缓存大小和每个线程要清除的数量（CHNKSZ），这些可能会略微调整性能。还可以尝试进行一些轻微的优化，例如进一步的位打包（如2、3、5、7轮而不仅仅是2、3、5轮），以减少大质数的执行时间，从而将打包减少约25%（可能将执行时间缩短到三分之二），并且通过轮阶乘预先清除页面段，直到17的因子为止，可以进一步减少约20%的时间，但这已经是纯C#代码所能做的极限，与C代码相比，C代码可以利用独特的本机代码优化。
无论如何，如果打算使用IEnumberable接口输出结果，就没有必要进行进一步的优化，因为大约三分之二的时间只用于枚举找到的质数，只有约三分之一的时间用于清除合数。更好的方法是编写类上的方法来实现所需的结果，例如CountTo、ElementAt、SumTo等，以完全避免枚举。
但我确实进行了额外的优化作为额外的答案，尽管增加了复杂性，但仍然显示出额外的好处，并进一步说明为什么不要使用SoA，因为完全优化的SoE实际上更好。
