20

我一直在尝试学习生成素数的算法,我在维基百科上遇到了阿特金筛子。我几乎了解算法的所有部分,除了少数。以下是问题:

  1. 下面的三个二次方程是如何形成的?4x^2+y^2、3x^2+y^2 和 3x^2-y2
  2. 维基百科中的算法讨论了模 60,但我不明白在下面的伪代码中如何/在哪里使用它。
  3. 这些提示 1、5、7 和 11 是如何找到的?

以下是来自维基百科的伪代码供参考:

// arbitrary search limit                                                   
limit ← 1000000                                                             

// initialize the sieve                                                     
for i in [5, limit]: is_prime(i) ← false                                    

// put in candidate primes:                                                 
// integers which have an odd number of                                     
// representations by certain quadratic forms                               
for (x, y) in [1, √limit] × [1, √limit]:                                    
    n ← 4x²+y²                                                              
    if (n ≤ limit) and (n mod 12 = 1 or n mod 12 = 5):                      
        is_prime(n) ← ¬is_prime(n)                                          
    n ← 3x²+y²                                                              
    if (n ≤ limit) and (n mod 12 = 7):                                      
        is_prime(n) ← ¬is_prime(n)                                          
    n ← 3x²-y²                                                              
    if (x > y) and (n ≤ limit) and (n mod 12 = 11):                         
        is_prime(n) ← ¬is_prime(n)                                          

// eliminate composites by sieving                                          
for n in [5, √limit]:                                                       
    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                                
        is_prime(k) ← false, k ∈ {n², 2n², 3n², ..., limit}                 

print 2, 3                                                                  
for n in [5, limit]:                                                        
    if is_prime(n): print n  
4

2 回答 2

39

您引用的 Wikipedia 文章中的阿特金筛法伪代码包含您的问题的答案或关于 Will Ness 提供链接的文章的讨论,尽管您可能无法将这些信息放在一起。简答如下:

  1. 这三个方程来自阿特金的数学证明,即所有素数都将作为这三个方程中的一个或多个的解出现,对于变量“x”和“y”的所有有效值具有适当的模条件,事实上他进一步证明了真正的素数将是那些对这三个方程有奇数解的数字(当从 False 的初始条件切换奇数次时保留为 True),每个方程的模条件不包括那些可以被平方整除的数字找到小于或等于测试数平方根的素数。

  2. 真正的阿特金筛法基于一组模 60 条件;此伪代码表示一种简化,其中使用一组模 12 条件 (5 × 12 = 60) 的每个方程的条件范围较小 - 但是,由于在此引入了冗余,因此需要完成 20% 的额外工作新的一组条件。这也是这个简化的伪代码需要从 5 而不是正常的 7 开始其素数扫描的原因,并且要从基本素数 5 开始进行基本素数平方免费消除,但会增加执行时间的成本,因为5 其他处理不当。这种简化的原因可能是牺牲了一些额外的代码开销,以减轻必须检查值的代码复杂性,

  3. “提醒”(应该拼写为“余数” )由您引用的伪代码中的“mod”运算符找到,它们代表可能使用的任何计算机语言中的 运算符,通常由计算机语言中的“%”符号表示,例如如 C、Java、C#、F# 等。该运算符产生整数余数 在被除数(第一个数字,在“mod”之前)除以除数(第二个数字,在“mod”符号之后)的整数除法之后。余数只是 4 而不是在完整模 60 算法中使用的 16 的原因是由于简化为模 12 算法。您会注意到,在模 12 条件下,“4x”条件通过了 25,这通常会被模 60 条件消除,因此许多包含 25 作为因子的复合材料需要通过额外的 5 平方素数自由检查来消除. 类似地,55 通过“3x+”检查,35 通过“3x-”检查,而对于完整的模 60 算法,它们不会。

正如 Wikipedia 文章“谈话”部分中所讨论的和上面所暗示的,这个引用的伪代码永远不会比仅基本赔率的埃拉托色尼筛 (SoE) 实现快得多,更不用说由于效率低下而使用相同程度的轮子分解的实现了:变量 'x' 和 'y' 在很多情况下(在文章中提到)不需要在整个范围内覆盖范围的平方根,正确使用模 60 轮可以恢复使用模 12 简化(如上所述),并且在二次方程的解中存在模格模式,因此使用(计算速度较慢的)模运算的条件不需要通过使用自动跳过那些不满足那些模条件的解的算法来测试晶格图案(仅在完整的 Atkin 和 Bernstein 论文中非常模糊地提及)。

要回答一个您没有问但应该问的问题:“为什么要使用阿特金筛?” .

实施阿特金筛法 (SoA) 而不是埃拉托色尼筛法 (SoE) 的主要原因是 SOA 更快的是“互联网知识”。这种信念有两个原因,如下所示:

  1. 假设 SoA 更快,因为它的渐近计算复杂度比 SoE 小一个因子 log(log n),其中 n 是筛选的素数范围。实际上,从 2 的 32 次方(40 亿以上)到 2 的 64 次方(大约 2 后跟 19 个零)的范围,这个因数是 6 比 5 等于 1.2。这意味着,与 32 位数字范围相比,在筛选到 64 位数字范围时,真正的 SoE 应该是预期的 1.2 倍,仅按线性比率计算,而如果一切都是理想的,SoA 将具有线性关系. 您可以理解,首先,对于数字范围如此巨大的差异来说,这是一个非常小的因素,其次,只有当两种算法在某个合理的素数范围内具有相同或接近相同的性能时,这种优势才有效。

    “互联网知识”没有清楚地理解的是,这些数字仅适用于比较给定范围内的性能比率与 相同算法的另一个给定范围内的性能比,而不是作为不同算法之间的比较。因此,它无法证明 SoA 将比 SoE 更快,因为对于特定 SoE 算法的给定筛选范围,SoA 可能以更大的开销开始,如以下优化的 SoE 示例所示。

  2. 由于 Atkin 和 Bernstein 根据他们在 Wikipedia 文章中链接的论文进行并发表的计算比较,SoA 被认为更快。虽然这项工作是准确的,但它仅适用于他们对他们比较的 SoE 算法所做的人为限制:由于 SoA 算法基于模 60 分解,它代表两个 2,3,5 分解车轮旋转,因此他们还限制了 SoE算法到相同的车轮分解。这样做,SoE 在 10 亿个测试范围内执行了大约 424,000 次合数剔除操作,而经过测试的完全优化的 SoA 具有大约 326,000 次组合切换和无正方形剔除操作,每个操作与每个合数剔除操作花费的时间大致相同在 SoE 中,因为是用相同的风格编写的。对于这组特定的车轮分解条件,这正是阿特金斯和伯恩斯坦比较测试所显示的。

    但是,SoA 不会响应进一步级别的车轮分解,因为 2、3、5 级别已“融入”算法,而 SoE 确实响应进一步级别:使用 2、3、5、7 车轮因式分解导致大约相同数量的操作,这意味着每个操作的性能大致相同。甚至可以使用比 2、3、5、7 级别更高的轮子分解级别,以使 SoE 的操作数量比 SoA 减少约 16.7%,从而获得更好的性能。实现这些额外级别的车轮分解所需的优化实际上比实现原始优化 SoA 的代码复杂性要简单。实现这些可比较的页面分段实现的内存占用大致相同,即页面缓冲区的大小加上基本素数数组。

    此外,两者都将受益于以“机器状态查找”样式编写,这将有助于使用内联代码和极端循环展开进行更好的优化,但由于算法更简单,SoE 比 SoA 更适合这种样式。

因此,对于高达约 32 位数字范围的筛子范围,最大优化的 SoE 比 SoA 快约 20%(或更多轮因式分解);然而,SoA 具有这种渐近计算复杂性的优势,因此它会在某个点上赶上。该点大约在 log (log n) 与 log (log (2^32)) 或 5 的比率为 1.2 或大约 2 乘以 10 的 19 次方的范围内——一个非常大的数字。 如果在现代台式计算机上运行优化的 SoA 需要大约三分之一秒来计算 32 位数字范围内的素数,并且如果该实现是理想的,因为随着范围的增加线性时间增加,计算这个范围的素数大约需要 45 年。然而,在这些高范围内的素数分析通常是在小块中完成的,对于非常大的筛子,SoA 的使用与 SoE 相比在理论上是实用的,但因数很小。

编辑_添加: 实际上,与较低的范围相比,页面分段的 SoE 和 SoA 都不会继续在这些巨大范围内以线性时间运行,因为两者都遇到了“切换”和“剔除”操作由于不得不跳过而失去效率的问题每次跳转大量页面;这意味着两者都需要修改算法来处理这种“页面跳转”,如果这样做的方式有任何细微差别,那么对 SoA 的微小优势可能会被完全抹去。SoA 在“跳表”中的项比 SoE 多得多,大约是找到的素数数量与该范围的平方根之间的反比,这可能会在处理过程中为两者添加一个 O(log n) 项,但对于在“跳转表”中具有更多条目数的 SoA,会增加一个常数因子。这个额外的事实几乎完全抵消了 SoA 相对于 SoE 的任何优势,即使对于非常大的范围也是如此。此外,SoA 在实现三个单独的二次方程的条件以及“无素数平方”循环而不只是素数剔除循环的更多情况下需要更多单独循环的持续开销,因此永远不会有如此低的完全优化时,每次操作的平均计算时间为 SoE。 END_EDIT_ADD

EDIT_ADD2: 我认为,关于阿特金筛子的大部分混淆是由于问题中引用的维基百科文章中伪代码的缺陷,因此提出了一个更好的伪代码版本来解决至少与“x”和“y”变量的范围以及模12与模60混淆有关的一些缺陷如下:

limit ← 1000000000        // arbitrary search limit

// Initialize the sieve
for i in {7,11,13,17,19,23,29,31, 37,41,43,47,49,53,59,61,...}:
    is_prime(i) ← false

// Put in candidate primes:
// integers which have an odd number of
// representations by certain quadratic forms.
while n ≤ limit, n ← 4x²+y² where x ∈ {1,2,...} and y ∈ {1,3,...} // odd y's
    if n mod 60 ∈ {1,13,17,29,37,41,49,53}:
        is_prime(n) ← ¬is_prime(n)
while n ≤ limit, n ← 3x²+y² where x ∈ {1,3,...} and y ∈ {2,4,...} // only odd 
    if n mod 60 ∈ {7,19,31,43}:                            // x's and even y's
        is_prime(n) ← ¬is_prime(n)
while n ≤ limit, n ← 3x²-y² where x ∈ {2,3,...} and y ∈ {x-1,x-3,...,1} //all 
    if n mod 60 ∈ {11,23,47,59}:                   // even/odd odd/even combos
        is_prime(n) ← ¬is_prime(n)

// Eliminate composites by sieving, only for those occurrences on the wheel
for n² ≤ limit where n ∈ {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
        while c ≤ limit, c ← n² × k where
                      k ∈ {1,7,11,13,17,19,23,29, 31,37,41,43,47,49,53,59,...}:
            is_prime(c) ← false

output 2, 3, 5
for n ≤ limit when n ← 60 × k + x where
  k ∈ {0..} and
  x ∈ {7,11,13,17,19,23,29,31, 37,41,43,47,49,53,59,61}:
    if is_prime(n): output n

上面看起来很简单,是一个很好的算法,除了它仍然没有比使用相同 2;3;5 分解轮的基本 Eratosthenes 筛更快,因为它浪费了几乎一半的内部循环切换操作失败模测试。展示:

以下是在 15 个垂直“x”值(每个值)和 15 个水平“y”奇数值范围内的 60 模 60 的重复 4x²+y² 模式;两者都从一个开始。请注意,它们关于垂直轴对称,但对于完整的 30 个“x”值范围,只有 225 个中的 128 个或 450 个中的 256 个是有效的切换操作:

  0 13 29 53  0  0 53 49 53  0  0 53 29 13  0
 17  0 41  0 37 17  0  1  0 17 37  0 41  0 17
 37  0  1  0  0 37  0  0  0 37  0  0  1  0 37
  0 13 29 53  0  0 53 49 53  0  0 53 29 13  0
 41 49  0 29  1 41 29  0 29 41  1 29  0 49 41
  0  0 49 13  0  0 13  0 13  0  0 13 49  0  0
 17  0 41  0 37 17  0  1  0 17 37  0 41  0 17
 17  0 41  0 37 17  0  1  0 17 37  0 41  0 17
  0  0 49 13  0  0 13  0 13  0  0 13 49  0  0
 41 49  0 29  1 41 29  0 29 41  1 29  0 49 41
  0 13 29 53  0  0 53 49 53  0  0 53 29 13  0
 37  0  1  0  0 37  0  0  0 37  0  0  1  0 37
 17  0 41  0 37 17  0  1  0 17 37  0 41  0 17
  0 13 29 53  0  0 53 49 53  0  0 53 29 13  0
  1  0  0 49  0  1 49  0 49  1  0 49  0  0  1

以下是在 5 个垂直“x”奇数值和 15 个水平“y”偶数值范围内重复 3x²+y² 模式的限定模 60。请注意,它们关于水平轴是对称的,但对于完整的 30 个“x”值范围,只有 75 个中的 48 个或 450 个中的 144 个是有效的切换操作,因为 450 个无效操作中的所有 144 个具有偶数“x”和奇数'y' 已经被淘汰:

  7 19  0  7 43  0 19 19  0 43  7  0 19  7  0
 31 43  0 31  7  0 43 43  0  7 31  0 43 31  0
 19 31  0 19  0  0 31 31  0  0 19  0 31 19  0
 31 43  0 31  7  0 43 43  0  7 31  0 43 31  0
  7 19  0  7 43  0 19 19  0 43  7  0 19  7  0

以下是在 5 个垂直“x”奇数值和 15 个水平“y”偶数值范围内重复 3x²-y² 模式的限定模 60。请注意,它们关于水平轴对称,但对于完整的 30 个“x”值范围,只有 75 个中的 48 个或 450 个中的 224 个是有效的切换操作:

 59 47  0 59 23  0 47 47  0 23 59  0 47 59  0
 23 11  0 23 47  0 11 11  0 47 23  0 11 23  0
 11 59  0 11  0  0 59 59  0  0 11  0 59 11  0
 23 11  0 23 47  0 11 11  0 47 23  0 11 23  0
 59 47  0 59 23  0 47 47  0 23 59  0 47 59  0

以下是在 5 个垂直“x”偶数值和 15 个水平“y”奇数值范围内以 60 为模的重复 3x²-y² 模式。请注意,它们关于垂直轴对称,但对于完整的 30 个“x”值范围,只有 75 个中的 48 个或 450 个中的 224 个是有效的切换操作:

 11  0 47 23  0 11 23  0 23 11  0 23 47  0 11
 47  0 23 59  0 47 59  0 59 47  0 59 23  0 47
 47  0 23 59  0 47 59  0 59 47  0 59 23  0 47
 11  0 47 23  0 11 23  0 23 11  0 23 47  0 11
 59  0  0 11  0 59 11  0 11 59  0 11  0  0 59

通过查看上表可以确定,对于上面的伪代码,总共有30个x值(包括偶数和奇数)的重复范围,在总共1125个组合操作中只有688个有效操作;因此,由于非生产性的跳过操作是最内层循环的一部分,因此它浪费了几乎一半的处理来有条件地跳过值。有两种可能的方法可以低效地消除“命中”冗余,如下所示:

  1. Atkin 和 Bernstein 论文中概述的方法,它使用“x”和“y”的组合模中识别的模式来仅处理模“命中”,使用的事实是,一旦在模式上找到给定的模,然后在某个模数(等于轮位位置)处有一个无限的值序列,其中每个模式由一个易于计算的(可变)偏移量分隔,其形式为“当前位置加上 A 乘以 'x' 的当前值加上 B ”和“当前位置加上 C 乘以 'y' 加上 D 的当前值”,其中 A、B、C 和 D 是简单的常数(简单的意思是小容易操作)。Bernstein 在许多查找表 (LUT) 的额外帮助下使用了该方法。

  2. 创建整体模式状态查找表 (LUT) 的方法(一个用于上述四个表中的每一个,加上一个用于次要素数平方自由循环)由“x”的模当前值与“y”的模组合索引找到要跳过的 A、B、C 和 D 参数,不是跳到下一个模式,而是跳到水平序列上的下一个位置。这可能是更高性能的选项,因为它更容易使用展开循环代码的内联来减少每个操作的极端时钟周期,并且当页面分段时,由于每个循环的跳转平均较小,将为大范围生成更有效的代码. 这可能会使每个循环的时钟周期接近于高度优化的埃拉托色尼筛(如在素筛中),但由于必须计算可变步长而不是能够为 SoE 使用固定的素数偏移,因此不太可能达到那么低。

因此,在减少初筛运行时间方面存在三个管理目标,如下所示:

  1. 一个成功的筛子减少了操作的数量,即使是“命中优化”的 SoA 与高度分解的 SoE 相比,在大约几十亿的范围内也减少了大约 16.7%。

  2. 一个成功的筛子减少了每个操作的 CPU 时钟周期,与高度优化的 SoE(如 primesieve)相比,SoA 失败了,因为它的操作由于可变增量而更加复杂,同样可能减少 20% 到 60%。Atkin 和 Bernstein 的素数 (SoA) 大约需要 4.5 个时钟周期,而素数 (SoE) 每次操作大约需要 2.5 个时钟周期,每个操作的范围为 10 亿,但也许可以通过借用一些优化技术来加快 SoA初筛。

  3. 成功的筛子具有相当低的代码复杂度,因此可以使用诸如“桶筛子”和其他页面分段优化之类的技术更容易地将其扩展到更大的范围。为此,阿特金筛子惨遭失败,因为它在页面分割大数字范围时变得更加复杂。编写一个可以与伯恩斯坦的素数 (SoA) 竞争的 SoA 程序非常复杂,更不用说与素数相媲美了,而编写与素数接近相同性能的代码却很容易。

  4. 一个成功的筛子具有较低的经验复杂性,SoA 确实比 SoE 高出 (log (log n) 倍,其中 n 是要筛分的上限,但这个额外的小比率可能不足以弥补对于上述前两个综合损失率,这个因素随着范围的增加而变小 。END_EDIT_ADD2

所以回答“为什么要使用阿特金筛子”这个问题?“如果 SoE 是通过最大轮分解优化实现的,直到筛出的数字非常大(64 位数字范围和更高),那么根本没有理由使用它,然后 SoA 优势非常小并且可能无法实现一切都取决于相对优化中的微小调整。” .

作为对另一个类似的阿特金筛子问题的回答,我在 https://stackoverflow.com/a/20559687/549617上发布了一个按上述优化的 SoE 的页面分段 C# 版本。

于 2013-12-13T17:54:56.447 回答
16

在我对这个问题的原始回答中,为了帮助更好地理解阿特金筛算法,我扩展了维基百科阿特金筛 (SoA) 伪代码,以纠正给出的简化“模 12”算法中的一些缺陷,与完整的“模 60”算法。同样,原始的 Wikipedia 伪代码算法导致实现效率低下;具体来说,“x”和“y”变量的范围都超出了必要的范围,模 12 的简化会导致额外的 20% 冗余工作。但是,正如该答案中所讨论的那样,该伪代码仍然效率不高,因为模测试仍在最里面的二次解循环中,因此对于那些未通过模测试的解决方案,这些解循环中几乎有一半是无效的;

Atkin 和 Bernstein 在他们的“primegen”程序中使用的避免循环冗余的算法可以通过基于我之前回答中的“命中”表的以下观察结果更好地理解:

  1. “命中”表的对称性取决于“y”是奇数(垂直对称)还是偶数(水平对称)。

  2. 因此,如果我们希望模式的前进沿水平方向进行,如图所示,我们可以翻转“3x+”和“3x-;odd x -> even”的“x”和“y”循环的顺序y" 案例表示奇数的 'x' 案例。

  3. 现在很容易消除两个翻转的“3x”表的零情况,它们在外循环中的条件水平只有五次重复,因为“y mod 3 等于 0”的情况加上中间的三个额外情况“轴”列,其中“y mod 5 等于 0”。

  4. 对于最后一个“3x-;偶数 x -> 奇数 y”的情况,只需过滤那些“(y + 2) mod 3 等于 0”的情况加上额外的两个情况,就可以很容易地消除大部分零最后一行“(x mod 60)等于59”,其中“(y + 2)mod 5等于0”(对称垂直轴列的重复消除,始终消除)。尽管这些“y 测试”似乎必须发生在内部循环中,因为垂直对称轴的每侧只有五个“命中”,但这些可以由“x”循环内联编码。

  5. 最难消除的是模式更复杂的“4x”表:有五个不同的水平模式,每个都可以识别,因为“第一列的(x mod 60)”是不同的,其中一个领先上面第一个表中的零行模式的模数为 5,另一个模数为 25;现在可以通过将案例“内联”为每个模式的垂直对称轴两侧的偏移量来确定五种不同的模式。

  6. 对于第一个表,到下一行的间距正好是 4 × ((x + 1)² - x²) 或 8x + 4,对于后两个(交换的)新表,是 (y + 2)² - y² 或 4y + 4 , 最后一个表为 3 × ((x + 2)² - x²) 或 12x + 12。

  7. 对于所有(现在)垂直对称的表格,中心轴偏移量可以从第一列的偏移量确定为 (y + 16)² - y² 或 32y + 256(对于长行)和 3 × ((x + 4) ² - x²) 或 24x + 48 用于短行。

  8. 以类似的方式,垂直轴两侧的对称偏移可以很容易地从垂直轴的偏移中计算出来,例如通过对 (y + 2)² - y² 或 4 + 4x 与 (y - 2) ² - y² 或 4 - 4y 其中 y 是长排案例的中心轴偏移量,正负两个位置;短行“x”案例的情况类似,但只是将“3x”“x”值的整体规模乘以 3。

  9. 上述由水平变量条件确定的过滤器似乎打破了我们消除内部循环中的条件代码的目标,但是当确定模式不是内部循环而是开始一系列新的内部循环时,该规则并不适用。模式中每个有效元素的阵列中水平方向的模式。这是有效的,因为我们可以在数学上证明下一个模式中的每个相同等效位置(具有相同的模数)按如下方式分开:对于长行情况,水平模式前进 30,模式由 (x + 30)² - x² 或 60x + 900 所以 60 × (x + 15); 对于短排情况,水平图案前进 10,图案由 3 × ((x + 10)² - x²) 分隔,因此 60 * (x + 5)。由于两者的模 60 为零,这就是为什么在每个位置都有相同模数的重复模式的确切解释。因此,一旦我们找到每个左上角位置的偏移量,我们就可以确定非零模式位置的起始位置的模数和偏移量,并使用从该点到筛阵列限制的模式关系滚动一个内部循环.

  10. 所有上述列和行增量都可以在不使用模加法的“除法”操作的情况下完成,其中算法只需对模数/商进行“溢出检查和调整”操作即可跟踪商和 60 的模数对在每次操作之后配对,但是执行“检查和调整”的条件代码可能比编写的方式在计算上更昂贵,为此,大多数现代优化编译器将使用溢出截断特性生成乘法运算,以替换计算上昂贵的除法运算用于除以小整数,这通常比运算组合加上“检查和调整”方法所需的条件代码或使用直接除法运算所需的时间更少。

  11. 上述消除组对于候选素数的压缩位表示非常有效,因为 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,如下所示:

  1. SoA 代码要复杂得多,并且更难编写和维护,因此更难实现页面分段,这对于更大的范围非常重要,以便利用最大的 CPU 能力和内存访问速度。

  2. 由于稍微复杂的内部循环代码和更多数量的外部循环,SoA 可能会在 CPU 时钟周期上损失至少 20%。

  3. 在 1 到 40 亿的相当大的范围内,SoA 与高度分解的 SoE(比上述讨论的更高的分解)相比具有稍多的操作。

  4. 随着范围变大,SoA 通过具有 log (log n) 的改进因子在非常大的范围内获胜,但该因素仅允许 SoA 在非常高的数字范围(可能大约10 的 19 次方)如果有的话,那么只有在相对复杂性保持不变的情况下,这是不可能的。

所以我再说一遍:“当埃拉托色尼筛在几乎所有方面都更好,特别是不太复杂,同时具有较少的恒定计算因子开销时,为什么要使用阿特金筛呢?" 基于我在这里的结论,我认为对于任何合理的筛分范围,SoA 都是劣等筛分,并且可能永远不会编写阿特金筛分的页面分段版本,但已经用各种计算机语言为筛分埃拉托色尼. 虽然阿特金筛在智力和数学上很有趣,但它并不实用,而且阿特金和伯恩斯坦的比较在过度限制他们在比较中使用的埃拉托色尼筛的实现的轮因式分解方面存在缺陷,以显示优势大约 30%,这实际上应该是大约 60%,除了阿特金筛子的持续实施开销。

于 2014-03-04T02:22:04.590 回答