在我对这个问题的原始回答中,为了帮助更好地理解阿特金筛算法,我扩展了维基百科阿特金筛 (SoA) 伪代码,以纠正给出的简化“模 12”算法中的一些缺陷,与完整的“模 60”算法。同样,原始的 Wikipedia 伪代码算法导致实现效率低下;具体来说,“x”和“y”变量的范围都超出了必要的范围,模 12 的简化会导致额外的 20% 冗余工作。但是,正如该答案中所讨论的那样,该伪代码仍然效率不高,因为模测试仍在最里面的二次解循环中,因此对于那些未通过模测试的解决方案,这些解循环中几乎有一半是无效的;
Atkin 和 Bernstein 在他们的“primegen”程序中使用的避免循环冗余的算法可以通过基于我之前回答中的“命中”表的以下观察结果更好地理解:
“命中”表的对称性取决于“y”是奇数(垂直对称)还是偶数(水平对称)。
因此,如果我们希望模式的前进沿水平方向进行,如图所示,我们可以翻转“3x+”和“3x-;odd x -> even”的“x”和“y”循环的顺序y" 案例表示奇数的 'x' 案例。
现在很容易消除两个翻转的“3x”表的零情况,它们在外循环中的条件水平只有五次重复,因为“y mod 3 等于 0”的情况加上中间的三个额外情况“轴”列,其中“y mod 5 等于 0”。
对于最后一个“3x-;偶数 x -> 奇数 y”的情况,只需过滤那些“(y + 2) mod 3 等于 0”的情况加上额外的两个情况,就可以很容易地消除大部分零最后一行“(x mod 60)等于59”,其中“(y + 2)mod 5等于0”(对称垂直轴列的重复消除,始终消除)。尽管这些“y 测试”似乎必须发生在内部循环中,因为垂直对称轴的每侧只有五个“命中”,但这些可以由“x”循环内联编码。
最难消除的是模式更复杂的“4x”表:有五个不同的水平模式,每个都可以识别,因为“第一列的(x mod 60)”是不同的,其中一个领先上面第一个表中的零行模式的模数为 5,另一个模数为 25;现在可以通过将案例“内联”为每个模式的垂直对称轴两侧的偏移量来确定五种不同的模式。
对于第一个表,到下一行的间距正好是 4 × ((x + 1)² - x²) 或 8x + 4,对于后两个(交换的)新表,是 (y + 2)² - y² 或 4y + 4 , 最后一个表为 3 × ((x + 2)² - x²) 或 12x + 12。
对于所有(现在)垂直对称的表格,中心轴偏移量可以从第一列的偏移量确定为 (y + 16)² - y² 或 32y + 256(对于长行)和 3 × ((x + 4) ² - x²) 或 24x + 48 用于短行。
以类似的方式,垂直轴两侧的对称偏移可以很容易地从垂直轴的偏移中计算出来,例如通过对 (y + 2)² - y² 或 4 + 4x 与 (y - 2) ² - y² 或 4 - 4y 其中 y 是长排案例的中心轴偏移量,正负两个位置;短行“x”案例的情况类似,但只是将“3x”“x”值的整体规模乘以 3。
上述由水平变量条件确定的过滤器似乎打破了我们消除内部循环中的条件代码的目标,但是当确定模式不是内部循环而是开始一系列新的内部循环时,该规则并不适用。模式中每个有效元素的阵列中水平方向的模式。这是有效的,因为我们可以在数学上证明下一个模式中的每个相同等效位置(具有相同的模数)按如下方式分开:对于长行情况,水平模式前进 30,模式由 (x + 30)² - x² 或 60x + 900 所以 60 × (x + 15); 对于短排情况,水平图案前进 10,图案由 3 × ((x + 10)² - x²) 分隔,因此 60 * (x + 5)。由于两者的模 60 为零,这就是为什么在每个位置都有相同模数的重复模式的确切解释。因此,一旦我们找到每个左上角位置的偏移量,我们就可以确定非零模式位置的起始位置的模数和偏移量,并使用从该点到筛阵列限制的模式关系滚动一个内部循环.
所有上述列和行增量都可以在不使用模加法的“除法”操作的情况下完成,其中算法只需对模数/商进行“溢出检查和调整”操作即可跟踪商和 60 的模数对在每次操作之后配对,但是执行“检查和调整”的条件代码可能比编写的方式在计算上更昂贵,为此,大多数现代优化编译器将使用溢出截断特性生成乘法运算,以替换计算上昂贵的除法运算用于除以小整数,这通常比运算组合加上“检查和调整”方法所需的条件代码或使用直接除法运算所需的时间更少。
上述消除组对于候选素数的压缩位表示非常有效,因为 16 个有效模值中的每一个都可以表示为 16 位字中的一位,并且字索引表示每个模 60 的索引车轮; 因此,我们可以通过仅推进整数个 16 位字来推进可被 60 整除的偶数增量。
对于此处的伪代码用作非页面分段算法,只需将模检查移出最内层循环就足够了,而无需消除结果中间循环的所有“命中未命中”,尽管它们仍然效率不是最高的,因为中间循环仍然处理“命中未命中”,这不会增加超过 1% 的运行时间;但是,对于分页实现,即使在外部循环中,“未命中”也会大大增加计算开销,因为这些循环在该页面的“x”值范围内每页运行一次,因此存在一个大约对数的因子与埃拉托色尼筛 (SoE) 相比,SoA 的筛分范围的平方根更多。
// arbitrary search limit
limit ← 1000000000
// the set of base wheel prime candidates
s ∈ {7,11,13,17,19,23,29,31, 37,41,43,47,49,53,59,61}
// Initialize the sieve as an array of wheels with
// enough wheels to include the representation for limit
for n ← 60 * w + x where w ∈ {0,...(limit - 7) ÷ 60}, x in s:
sieve(n) ← false // start with no primes.
// Put in candidate primes:
// integers which have an odd number of
// representations by certain quadratic forms.
while n ≤ limit when n ← 4x²+y² // all x's and only odd y's
where x ∈ {1,...}, y ∈ {1,31,...}: // skip by pattern width in y
for cb ≤ limit when midy ← y + i, cb ← n - y² + midy²
where i ∈ {0,2,...28}: // middle level loop to "hit" test modulos
cbd ← (cb - 7) ÷ 60, cm ← cb mod 60
if cm in {1,13,17,29,37,41,49,53}:
while c ≤ limit when c ← 60 × (cbd - midy² + ((midy + 15 × j) × j)²) + cm
where j {0,...}:
sieve(c) ← ¬sieve(c) // toggles the affected bit
while n ≤ limit when n ← 3x²+y² // only odd x's and even y's
where x ∈ {1,3,...}, y ∈ {2,32,...}: // skip by pattern width in y
for cb ≤ limit when midy ← y + i, cb ← n - y² + midy²
where i ∈ {0,2,...28}: // middle level loop to "hit" test modulos
cbd ← (cb - 7) ÷ 60, cm ← cb mod 60
if cm in {7,19,31,43}:
while c ≤ limit when c ← 60 × (cbd - midy² + ((midy + 15 × j) × j)²) + cm
where j {0,...}:
sieve(c) ← ¬sieve(c) // toggles the affected bit
while n ≤ limit when n ← 3x²-y² //even/odd & odd/even combos
where x ∈ {2,3,...}, y ∈ {x-1,x-31,...,1}: // skip by pattern width in y
for midy ≥ 1 and cb ≤ limit when midy ← y - i, cb ← n + y² - midy²
where i ∈ {0,2,...28}: // middle level loop to "hit" test modulos
cbd ← (cb - 7) ÷ 60, cm ← cb mod 60
if cm in {11,23,47,59}:
while yn ≥ 1 and c ≤ limit when yn = midy - 30 * j and
c ← 60 × (cbd + midy² - ((midy - 15 × j) × j)²) + cm
where j {0,...}:
sieve(c) ← ¬sieve(c) // toggles the affected bit
// Eliminate prime squares by sieving, only for those occurrences on the wheel
for sqr ≤ limit where sqr ← n², n in {7,11,13,17,19,23,29,31, 37,41,43,47,49,53,59,61,...}:
if is_prime(n):
// n is prime, omit multiples of its square; this is
// sufficient because composites which managed to get
// on the list cannot be square-free
if c0 ≤ limit where c0 ← sqr × (m mod 60) where m in s:
c0d ← (c0 - 7) ÷ 60, cm ← (c0 - 7) mod 60 + 7 // means cm in s
while c ≤ limit for c ← 60 × (c0d + sqr × j) + cm
where j ∈ {0,...}:
sieve(c) ← false // cull composites of this prime square
output 2, 3, 5
for n ≤ limit when n ← 60 × k + x where k ∈ {0..}, x in s:
if sieve(n): output n
特别注意,最里面的二次循环非常简单,但每个循环以线性增加的变量递增,而使用 Eratosthenes 级数的素数平方自由消除以固定量“sqr”递增。这些简单的循环将在现代 CPU 上快速执行,如果内存访问速度足够快,例如小筛子范围或使用页面分段版本(其中页面适合 L2/L1 CPU 缓存),则每个循环大约需要 4 个 CPU 时钟周期具有大约此范围的内存访问时间。
当阵列大到许多兆字节时,平均内存访问时间可能平均在 20 到 100 个 CPU 时钟周期之间(取决于 RAM 的速度、CPU 的内存访问效率和使用效率)中间缓存),这将是执行速度的主要限制因素,而不是循环的紧密性,而不是算法所需的切换/剔除操作的数量。这些特殊的算法,包括 SoA 和 SoE,对大型缓冲区的内存接口造成了特殊的压力,因为它们跳过缓冲区,每次增加整个轮周乘以乘数,并重复不同的偏移量,而不是处理所有“依次撞击”车轮;这是由于将“命中”测试移动到中间循环,这节省了“命中”测试,但结果是这个更大的缓冲区“跳跃”,SoA的平均“跳跃”高于SoE。因此,该算法在数字到 10 亿的范围内的执行时间估计约为每次操作 60 个 CPU 时钟周期(由于较小的平均“跳跃”,SoE 将减少一个因子)乘以切换/剔除操作的数量(对于 Atkin 筛来说约为 260,000)除以 CPU 时钟速度。因此,3.5 GHz CPU 将在大约 4 秒内在十亿范围内执行该算法,SoA 的非分页实现和这种风格的 SoE 的执行速度之间的主要差异是大型内存访问效率的差异。缓冲器,
鉴于对于大缓冲区,内存访问时间是限制速度因素,因此找到具有最少操作数的算法变得非常重要。Atkin 和 Bernstein 在他们的“primegen”示例程序中实现的 SoA 将“命中”操作的常数因子作为与筛选范围的比率,并且通过从我的第一个答案中实现一个简单版本的程序,我们可以确定这个因子大约是筛选的主要范围的 0.26:这意味着对于十亿的筛选范围,大约有 260,000 次操作。
对于使用 210 元素轮子的高轮子分解实现的 SoE,每次旋转 48 个素数候选“命中”(素数为 2;3;5;7 的小轮子)结合额外素数 11 的预剔除;13;17;19,筛选范围 n 的剔除操作数约为 n 倍 (ln((ln n) / 2) + M - 1/(ln n) - 1/2 - 1/3 - 1 /5 - 1/7 - 1/11 - 1/13 - 1/17 - 1/19) * 48 / 210 其中“*”表示乘以,“/”表示除以,M 是 Meissel-Mertens 常数 M ≈ 0.2614972128... 校正由 (ln ln n) 估计产生的操作数,倒数“ln n”项是从基本素数的平方开始的剔除调整(小于范围的平方根) n); 对于 10 亿的 n,与 n 相比的因子约为 0。
这与 SoA 的操作数量大致相同,并与 Atkin 和 Bernstein 比较测试中使用的简单 2;3;5 轮形成对比,后者的系数为 0.404805008 或约 400,000 次操作,筛分范围为 10 亿。对于更高分解的 SoE,这个更好的比率意味着这个完全可实现的实现的常数因子与 SoA 的常数因子大致相同,除了由于循环复杂性导致的常数因子的任何差异。
EDIT_ADD1: 为了比较,以下是高度分解和预先剔除的 SoE 的等效伪代码,其编写方式与上述 SoA 算法相同:
// arbitrary search limit
limit ← 1000000000
// the set of base wheel primes
r ∈ {23,29,31,37,41,43,47,53, 59,61,67,71,73,79,83, //positions + 19
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}
// an array of length 11 times 13 times 17 times 19 = 46189 wheels initialized
// so that it doesn't contain multiples of the large wheel primes
for n where n ← 210 × w + x where w ∈ {0,...46189}, x in r: // already
if (n mod cp) not equal to 0 where cp ∈ {11,13,17,19}: // no 2,3,5,7
mstr(n) ← true else mstr(n) ← false // factors
// Initialize the sieve as an array of the smaller wheels with
// enough wheels to include the representation for limit
for n where n ← 210 × w + x, w ∈ {0,...(limit - 19) ÷ 210}, x in r:
sieve(n) ← mstr(n mod (210 × 46189)) // init pre-culled primes.
// Eliminate composites by sieving, only for those occurrences on the
// wheel using wheel factorization version of the Sieve of Eratosthenes
for n² ≤ limit when n ← 210 × k + x where k ∈ {0..}, x in r
if sieve(n):
// n is prime, cull its multiples
s ← n² - n × (x - 23) // zero'th modulo cull start position
while c0 ≤ limit when c0 ← s + n × m where m in r:
c0d ← (c0 - 23) ÷ 210, cm ← (c0 - 23) mod 210 + 23 //means cm in r
while c ≤ limit for c ← 210 × (c0d + n × j) + cm
where j ∈ {0,...}:
sieve(c) ← false // cull composites of this prime
output 2, 3, 5, 7, 11, 13, 17, 19,
for n ≤ limit when n ← 210 × k + x where k ∈ {0..}, x in r:
if sieve(n): output n
请注意,尽管此代码被高度分解,但它的复杂程度要低得多. 还要注意必须为 SoA 代码运行多少外部循环,其中外部循环的总数与 SoA 的筛分范围的平方根几乎相同,其中筛分范围除以SoE 的该范围外环的平方根的自然对数 - 对于 10 亿筛分范围的 SoA,SoE 的因子大约是 SoA 的十倍。这是 SoA 所做的权衡:更多的外部循环,每个内部循环的操作更少,缓冲区“跳跃”更大 操作之间;然而,这也是使 SoA 每次操作效率较低的原因,尤其是对于较大的筛网范围。 END_EDIT_ADD1
至于更大范围的代码复杂性,这些总是作为页面分段筛来实现,以避免将整个筛阵列放在 RAM 内存中(内存限制会限制范围)并利用 CPU 缓存局部性以减少如上所述的平均内存访问时间。对于 SoE,这可以根据保存的基本素数表示的副本轻松实现,其中大约有范围的平方根除以该范围的平方根数的对数;对于 SoA,需要参考相同的基本素数表示以进行素数平方自由消除,但还需要使用序列偏移值加上保存的“x”或“y”值(或保存两个“x”的当前值'和'y' ) 以便在下一页偏移量处恢复,保存的序列总数约为范围的平方根(不除以根 n 的自然对数),或者新的 'x'/'y' 对偏移量可以是以相当大的计算开销为每个新页面的开始计算。第一种方法大大增加了内存使用开销,特别是对于需要为每个线程保留一份副本的多线程实现;对于这些更大的范围,与 SoE 相比,第二个不使用更多内存但使用更多计算能力。可以以相当大的计算开销为每个新页面的开始计算对偏移量。第一种方法大大增加了内存使用开销,特别是对于需要为每个线程保留一份副本的多线程实现;对于这些更大的范围,与 SoE 相比,第二个不使用更多内存但使用更多计算能力。可以以相当大的计算开销为每个新页面的开始计算对偏移量。第一种方法大大增加了内存使用开销,特别是对于需要为每个线程保留一份副本的多线程实现;对于这些更大的范围,与 SoE 相比,第二个不使用更多内存但使用更多计算能力。
为了充分利用现代多核 CPU,页面分段算法也很重要,它可以通过使用的内核数量的比率减少给定任务所需的时钟时间,特别是对于每个页面处理都可以进行的页面分割算法分配给单独的线程。
与 Atkin 和 Bernstein 的 SoA 的“primegen”示例实现的大约 4.5 个 CPU 时钟周期相比,“primesieve”SoE 实现管理大约 2.5 个 CPU 时钟周期的平均操作时间的方式是通过极端内联循环展开页面分割;这种技术也可以应用于 SoA,但是这种特殊的伪代码不适合这样做,并且需要使用在单个循环(每个二次方程循环)中处理所有模数的“命中率”优化的替代技术;但是,这是一个复杂的算法,似乎超出了这个问题的范围。可以说,与将这种额外优化添加到 SoA 所需的相比,复杂性的跳跃到目前为止微不足道,
因此,虽然像在原始维基百科伪代码中那样实现一个幼稚的 SoA 算法,或者甚至在这两个答案中使用更复杂的 SoA 伪代码是相当容易的,但编写一些甚至可以与 Bernstein 的“primegen”竞争的东西是非常复杂的" (SoA) 更不用说更快的 "primesieve" (SoE),而编写与 "primegen" 性能接近的代码很容易,而编写与 "primesieve" 竞争的东西只需要多一点工作使用稍作修改的 SoE 算法来支持页面分割。
EDIT_ADD2: 为了将这个理论付诸实践,我通常会编写一个示例程序来演示这些原理;然而,在这种情况下,我们在“primegen”源代码中拥有我们真正需要的一切,这是 SoA 的发明者的参考实现。在我的机器(i7-2600K,3.5 GHz)上,缓冲区大小从 8192 调整为 32768 以反映我的 L1 缓存大小,这个“primespeed”在大约 0.365 秒内将素数筛选到十亿。事实上,Atkin 和 Bernstein 在与他们的 Eratosthenes 筛算法相比时有点作弊,因为他们只为其分配了大约 4 Kilobyte 缓冲区。将该缓冲区大小调整到大致相同的大小(“#define B32 1001”到“#define B32 8008”)会导致“eratspeed”的执行时间大约为 0.42 秒,达到大约 10 亿,但仍然稍微慢一些。然而,正如上面所讨论的,它正在做更多的工作,因为它做了大约 4.048 亿次剔除操作,而 SoA 只做了大约 0 次。260 亿次组合切换/剔除操作。使用上述 SoE 伪代码中的算法来改变这一点,具有最大车轮分解意味着 SoE 实际上会比 SoA 执行的操作略少,约为 2.505 亿,这意味着执行时间将减少到大约三分之二或更多或大约 0.265 到 0.3 秒,使其比“primegen”更快。
同时,如果将“primespeed”与“eratspeed”代码进行比较,就会发现分页的 SoE 远没有 SoA 代码复杂,正如上面的伪代码所示。
当然,人们可能会反驳说,SoA 将筛分范围操作的 0.26 倍与人们想要的筛分范围保持在一个比值,而 SoE 比率根据上面给出的等式攀升——筛分范围增加了大约 20%甚至达到一千亿。然而,对两种算法进行轻微调整以增加缓冲区大小并将其“子切片”以在 L1 缓存大小的切片中处理它,SoE 将按预期执行,而 SoA 的开销由于中间的额外比率而继续上升与 SoE 相比,SoA 的 /inner 循环。虽然 SoE 的剔除操作增加了大约 20%,但 SoA 的二次处理循环增加了大约 20%,如果有的话,净增益非常小。如果对 SoA 进行了这些优化,那么它可能只是在一些非常大的筛分范围(例如 10 到 19 次方左右)上赶上最大轮式分解 SoE,但我们正在谈论的是在台式计算机上进行数百年的处理,即使使用多处理也是如此。
SoA可能适用于从非常大的范围偏移(例如 10 次方的 20 次方)中筛选子范围(即单个页面段),但在处理过程中仍然存在额外内存或计算的成本使用 SoE 不需要的所有“x”和“y”值。 END_EDIT_ADD2
就像我的第一个答案一样,比较表明,SoA 在几乎所有方面都输给了高度轮式分解的 SoE,如下所示:
SoA 代码要复杂得多,并且更难编写和维护,因此更难实现页面分段,这对于更大的范围非常重要,以便利用最大的 CPU 能力和内存访问速度。
由于稍微复杂的内部循环代码和更多数量的外部循环,SoA 可能会在 CPU 时钟周期上损失至少 20%。
在 1 到 40 亿的相当大的范围内,SoA 与高度分解的 SoE(比上述讨论的更高的分解)相比具有稍多的操作。
随着范围变大,SoA 通过具有 log (log n) 的改进因子在非常大的范围内获胜,但该因素仅允许 SoA 在非常高的数字范围(可能大约10 的 19 次方)如果有的话,那么只有在相对复杂性保持不变的情况下,这是不可能的。
所以我再说一遍:“当埃拉托色尼筛在几乎所有方面都更好,特别是不太复杂,同时具有较少的恒定计算因子开销时,为什么要使用阿特金筛呢?" 基于我在这里的结论,我认为对于任何合理的筛分范围,SoA 都是劣等筛分,并且可能永远不会编写阿特金筛分的页面分段版本,但已经用各种计算机语言为筛分埃拉托色尼. 虽然阿特金筛在智力和数学上很有趣,但它并不实用,而且阿特金和伯恩斯坦的比较在过度限制他们在比较中使用的埃拉托色尼筛的实现的轮因式分解方面存在缺陷,以显示优势大约 30%,这实际上应该是大约 60%,除了阿特金筛子的持续实施开销。