我不知道这是否可行，但我只是想问一下。我的数学和算法能力使我有些失望：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();
}
}
现在，我想摆脱限制，以便该序列永远不会停止（理论上）。
我认为它可能会像这样：
从一些琐碎的限制开始
找到所有素数到极限
产生所有新发现的素数
增加限制（通过将旧限制加倍或加倍或类似的方法）
转到步骤2
最理想的是，第二步只需在旧限制和新限制之间进行即可。换句话说，它不必一次又一次找到最低的素数。
有办法吗？我的主要问题是我不太了解此算法中的
x和y。就像，我可以使用相同的算法，但是将x和y设置为oldLimit（最初为1）并运行到newLimit吗？或如何运作？有什么聪明的想法可以阐明吗？这样做的目的是让我不必设置该限制。这样我就可以使用Linq并仅使用
Take()无论我需要多少素数，而不必担心该限制是否足够高，等等。
最佳答案
这是C＃中无界素数筛选的一种解决方案，可以使用Eratosthenes筛（SoE）或Atkin筛（SoA）算法实现。但是，我认为，给出的优化SoA解决方案的极端复杂性几乎不值得，而真正的SoE却可以提供相同的性能而又没有那么复杂。因此，这也许只是部分答案，尽管它展示了如何实现更好的SoA算法，并展示了如何使用SoE实现无限序列，但仅暗示了如何组合这些元素以编写合理有效的SoA。
请注意，如果只需要讨论这些想法的最快实现，请跳至该答案的底部。
首先，我们应该对本练习的观点进行评论，以产生无穷大的质数序列，以允许使用IEnumerable扩展方法，例如Take（），TakeWhile（），Where（），Count（）等，因为这些方法会浪费一些由于增加了方法调用级别，因此性能提高，每个调用至少增加28个机器周期，并返回枚举一个值，并且每个函数添加多个级别的方法调用；也就是说，即使对枚举的结果使用更多的命令式过滤技术以提高速度，具有（有效）无限的序列仍然有用。
请注意，在下面的简单示例中，我将范围限制为无符号32位数字（uint），因为超出了范围，基本的SoE或SoA实现对于效率并不合适，因此需要切换使用“桶”或其他形式的增量筛进行部分筛分以提高效率；但是，对于最快实现中的完全优化的页面分段筛，该范围会增加到64位范围，尽管如写的那样，它可能不希望在运行时筛分超过约100万亿（十至十四次幂）整个过程可能需要数百年的时间。
由于问题可能是出于错误的原因选择了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秒内，在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模型，并根据该文章中的建议对各个二次解使用单独的循环对the Wikipedia article中的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，我们应该选择CPU L1或L2高速缓存大小的数组大小至少为16 KB，以避免在我们从数组中剔除素数的组合时避免高速缓存发生抖动，这意味着BitArray可以具有一个大小的八倍（每字节8位）或128 Kilobits。由于我们只需要筛选奇数质数，这代表了超过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快，而且要复杂得多，但我并没有真正明白这一点；）
EDIT_ADD：
注意：以下代码已得到纠正，因为原始发布的代码对于提供的素数序列是正确的，但它比范围平方根以下的所有数字都被剔除的速度要慢一半找到的底数会填满范围的平方根。
最有效和最简单的优化是将每个段页面的剔除操作委托给后台线程，这样，如果有足够的CPU核心，枚举素数的主要限制是枚举循环本身，在这种情况下，所有32位数字中的素数可以在上述参考计算机（在C＃中）上在大约十秒钟内枚举范围，而无需其他优化，而所有其他优化包括编写IEnumerable接口实现，而不是使用内置的迭代器将所有203,280,221个素数减少到大约六秒钟在32位数字范围（1.65秒到十亿）中，大部分时间只是在枚举结果。下面的代码进行了许多修改，包括使用Tasks使用的Dotnet Framework 4 ThreadPool中的后台任务，使用实现为查找表的状态机来实现进一步的位打包以使质数枚举更快，并重写算法作为使用“自己动手”迭代器的可枚举类，以提高效率：
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(); } }
请注意，由于设置和初始化数组的开销，此代码并不能快于上一个素数范围的最新版本（最多一两百万），而对于不超过40亿以上的较大范围，此代码却要快得多。它比问题的阿特金筛子实施速度快约8倍，但对于更大范围（高达40亿多）的速度，则要快20到50倍。在代码中定义了常量，用于设置基本缓存大小以及每个线程要删除的常量数（CHNKSZ），这可能会稍微影响性能。可以尝试进行一些细微的进一步优化，这些优化可能将大质数的执行时间减少多达两倍，例如对2,3,5,7轮进行进一步的位打包，而不是仅对2,3,5轮进行位打包。减少了约25％的打包（也许将执行时间减少到三分之二），并通过轮乘分解将页面段预选为因子17，从而进一步减少了约20％，但是这些程度与可以利用独特的本机代码优化的C代码相比，纯C＃代码可以完成的工作。
无论如何，如果一个问题打算使用IEnumberable接口进行输出，则几乎不值得进行进一步的优化，因为问题需要这样做，因为大约三分之二的时间仅用于枚举找到的素数，而只有大约三分之一的时间被用于剔除复合数字。更好的方法是在类上编写实现所需结果的方法，例如CountTo，ElementAt，SumTo等，以便完全避免枚举。
但是，我确实做了additional answer的额外优化，尽管增加了复杂性，但仍然显示了其他好处，并且进一步表明了为什么人们不希望使用SoA，因为经过完全优化的SoE实际上更好。
关于c# - C#:如何使Atkin筛网增量，我们在Stack Overflow上找到一个类似的问题： https://stackoverflow.com/questions/1569393/