Svi*_*ish 3 c# algorithm primes sieve sieve-of-atkin
我不知道这是否可能，但我只是想问一下。我的数学和算法技能有点让我失望：P
问题是我现在有这个类可以生成达到一定限制的素数：
public class Atkin : IEnumerable<ulong>
{
private readonly List<ulong> primes;
private readonly ulong limit;
public Atkin(ulong limit)
{
this.limit = limit;
primes = new List<ulong>();
}
private void FindPrimes()
{
var isPrime = new bool[limit + 1];
var sqrt = Math.Sqrt(limit);
for (ulong x = 1; x <= sqrt; x++)
for (ulong y = 1; y <= sqrt; y++)
{
var n = 4*x*x + y*y;
if (n <= limit && (n % 12 == 1 || n % 12 == 5))
isPrime[n] ^= true;
n = 3*x*x + y*y;
if (n <= limit && n % 12 == 7)
isPrime[n] ^= true;
n = 3*x*x - y*y;
if (x > y && n <= limit && n % 12 == 11)
isPrime[n] ^= true;
}
for (ulong n = 5; n <= sqrt; n++)
if (isPrime[n])
{
var s = n * n;
for (ulong k = s; k <= limit; k += s)
isPrime[k] = false;
}
primes.Add(2);
primes.Add(3);
for (ulong n = 5; n <= limit; n++)
if (isPrime[n])
primes.Add(n);
}
public IEnumerator<ulong> GetEnumerator()
{
if (!primes.Any())
FindPrimes();
foreach (var p in primes)
yield return p;
}
IEnumerator IEnumerable.GetEnumerator()
{
return GetEnumerator();
}
}
现在，我想要的是摆脱限制，以便序列永远不会停止（理论上）。
我想事情可能会是这样的：
最佳情况下，第二步应该只在旧限制和新限制之间起作用。换句话说，它不应该一次又一次地找到最低的素数。
有办法做到这一点吗？我的主要问题是我不太明白这个算法中的x例如y是什么。比如，我可以使用相同的算法类型，但将x和设置y为oldLimit（最初为 1）并运行到 吗newLimit？或者说这会如何运作？有聪明人对此有一些启发吗？
这样做的目的是让我不必设置该限制。这样我就可以使用 Linq 以及Take()我需要的任意数量的素数，而不用担心限制是否足够高等等。
下面是 C# 中无界素数筛选的解决方案，可以使用埃拉托斯特尼筛法 (SoE) 或阿特金筛法 (SoA) 算法实现；然而，我认为，与真正的 SoE 提供大约相同的性能而没有那么多复杂性相比，优化 SoA 解决方案的极端复杂性几乎不值得。因此，这可能只是部分答案，因为虽然它展示了如何实现更好的 SoA 算法并展示了如何使用 SoE 实现无界序列，但它仅暗示如何组合这些来编写相当高效的 SoA。
请注意，如果只需要讨论这些想法的最快实现，请跳至此答案的底部。
首先，我们应该评论一下这个练习的要点，即生成无界素数序列，以允许使用 IEnumerable 扩展方法，例如 Take()、TakeWhile()、Where()、Count() 等，因为这些方法给出了一些由于增加了方法调用级别，为每个调用增加了至少 28 个机器周期，并返回枚举一个值，并为每个函数添加了多个级别的方法调用，从而提高了性能；也就是说，即使对枚举结果使用更多命令式过滤技术以获得更快的速度，拥有（有效）无限序列仍然有用。
请注意，在下面的更简单的示例中，我将范围限制为无符号 32 位数字 (uint)，超过该范围基本 SoE 或 SoA 实现在效率方面并不真正合适，需要进行切换到“桶”或其他形式的增量筛以进行部分筛分以提高效率；然而，对于完全优化的页分段筛选（如最快的实现），范围增加到 64 位范围，尽管正如所写，运行时可能不希望筛选超过一百万亿（十的十四次方）完整范围需要长达数百年的时间。
由于问题选择 SoA 的原因可能是错误的，首先将 Trial Division (TD) 素数筛误认为是真正的 SoE，然后在 TD 筛上使用朴素的 SoA，让我们建立真正的有界 SoE，它可以在几个中实现lines 作为方法（可以根据问题的实现编码风格转换为类），如下所示：
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 位数字范围内的 203,280,221 个素数计算到 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 模型，并优化了Wikipedia 文章中的 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 实现实现此性能的复杂性。我的观点是“当人们非常努力地工作以获得与 SoE 相同的性能时，为什么选择 SoA？”。
现在，对于无界素数序列，最简单的方法是定义 Uint32.MaxValue 的 const top_number 并消除 primesSoE 或 primesSoA 方法中的参数。这有点低效，因为无论正在处理的实际素数值如何，剔除仍然在整个数字范围内完成，这使得确定小范围素数花费的时间比应有的时间长得多。除了这种情况和极端的内存使用之外，还有其他原因需要使用素数筛的分段版本：首先，使用主要位于 CPU L1 或 L2 数据缓存中的数据的算法由于内存访问效率更高而处理速度更快，并且其次，因为分段允许将作业轻松分割为多个页面，这些页面可以使用多核处理器在后台进行剔除，从而实现与所使用的核心数量成正比的加速。
为了提高效率，我们应该选择 CPU L1 或 L2 缓存大小的数组大小，对于大多数现代 CPU 而言，该数组大小至少为 16 KB，以避免在我们从数组中剔除素数组合时发生缓存抖动，这意味着 BitArray 可以具有大小是其大小的八倍（每字节 8 位）或 128 KB。由于我们只需要筛选奇数素数，这代表了超过 25 万的数字范围，这意味着分段版本将仅使用 8 个分段来筛选 200 万个，如欧拉问题 10 的要求。人们可以将结果保存在最后一段并从该点继续，但这将妨碍将此代码适应多核情况，在多核情况下，将剔除到多线程的后台以充分利用多核处理器。（单线程）无界 SoE 的 C# 代码如下：
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 转换为 SoA，我们只需要从有界情况更改为 SoA 剔除代码，但修改为每个页段重新计算起始地址，就像为 SoE 情况计算起始地址一样，但运算起来更加复杂，因为二次方程按数字的平方而不是简单的素数倍数前进。我将把所需的调整留给学生，尽管我并没有真正明白这一点，因为经过合理优化的 SoA 并不比当前版本的 SoE 更快，而且要复杂得多；）
编辑_添加：
注意：下面的代码已被更正，虽然原始发布的代码对于提供的素数序列是正确的，但它比所需的速度慢了一半，因为它剔除了范围平方根以下的所有数字，而不是仅剔除找到的碱基素数直到范围的平方根。
最有效和最简单的优化是将每个段页的剔除操作委托给后台线程，这样，在有足够的 CPU 核心的情况下，枚举素数的主要限制是枚举循环本身，在这种情况下，32 位数字中的所有素数在上述参考机器（C# 中）上可以在大约十秒内枚举范围，无需其他优化，所有其他优化包括编写 IEnumerable 接口实现而不是使用内置迭代器，将所有203,280,221 个素数的时间减少到大约六秒在 32 位数字范围内（1.65 秒到十亿），同样大部分时间都花在枚举结果上。以下代码进行了许多修改，包括使用任务使用的 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; //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(); } }
请注意，由于设置和初始化数组的开销，对于小范围的素数（最多一两百万），此代码并不比上一个版本快，但对于高达四十亿以上的较大范围来说，速度要快得多。它比阿特金筛问题的实施快大约 8 倍，但对于高达 40 亿以上的更大范围，速度会快 20 到 50 倍。代码中定义了常量，用于设置基本缓存大小以及每个线程要剔除的缓存数量 (CHNKSZ)，这可能会稍微调整性能。可以尝试一些进一步的优化，这可能会将大素数的执行时间减少多达两倍，例如针对 2,3,5,7 轮进行进一步的位打包，而不仅仅是 2,3,5 轮打包工作减少了约 25%（可能将执行时间减少到三分之二），并且通过轮阶乘预先剔除页面段高达 17 倍，进一步减少了约 20%，但这些只是程度与可以利用独特的本机代码优化的 C 代码相比，纯 C# 代码可以完成哪些工作。
无论如何，如果您打算按照问题的要求使用 IEnumberable 接口进行输出，则几乎不值得进一步优化，因为大约三分之二的时间用于枚举找到的素数，而只有大约三分之一的时间用于剔除合数。更好的方法是在类上编写方法来实现所需的结果，例如 CountTo、ElementAt、SumTo 等，从而完全避免枚举。
但我确实做了额外的优化作为额外的答案，尽管复杂性增加了，但它显示了额外的好处，并且进一步说明了为什么人们不想使用 SoA，因为完全优化的 SoE 实际上更好。