Atkin筛法的分段算法,可行吗?

6

我知道埃拉托色尼筛法可以被实现为无上限的连续寻找质数(分段筛法)。

我的问题是,阿特金/伯恩斯坦筛法能否以同样的方式实现?

相关问题:C#:如何使阿特金筛法增量

然而相关问题只有一个回答,它说“所有筛法都不可能”,这显然是不正确的。


2
我已经研究了Atkin/Bernstein筛法多年,但从未想到如何将其分段 - 我的意思是,从某个任意大的数字开始,并进行一些较小的预计算。如果有人能做到这一点,我会很感兴趣。 - librik
2个回答

4

Atkin/Bernstein在他们的原始论文的第5节中提供了一个分段版本。据推测,Bernstein的primegen程序使用了这种方法。


谢谢,我感到很惭愧,因为我没有从头到尾地阅读原始论文。 - K.Steff
2
阿特金和伯恩斯坦在他们的论文中没有提到的是,在使用大范围所需的页面分割时,维护Atkin筛法的效率极其困难:1)素数平方自由步骤很快变得比一页要大得多,需要一些复杂性(并增加了低效率)来处理这些问题。2)随着范围的增加,每个用于二次解的“x”和“y”变量的新页面起点的计算变得越来越耗时,因此O(n)相对性能从未发生。 - GordonBGood

2
事实上,可以实现一个无界的Atkin筛法(SoA),而不必像我在F#中所做的那样使用分段。请注意,这是一个纯函数版本,甚至不使用数组来组合二次方程和无平方数过滤器的解决方案,因此比更加命令式的方法要慢得多。
Berstein使用查找表进行最佳32位范围的优化会使代码变得极为复杂,不适合在此展示,但很容易调整我的F#代码,使序列从设定的下限开始,并仅在一定范围内使用,以实现分段版本,并/或将相同的技术应用于使用数组的更加命令式的方法。
请注意,即使是Berstein的SoA实现,也并不比Sieve of Eratosthenes的所有可能的优化快,如Kim Walisch's primesieve所述,但只比他的实现中所选范围内的Sieve of Eratosthenes的等效优化版本快。
“EDIT_ADD: 对于那些不想深入研究Bernstein的伪代码和C代码的人,我在这个答案中添加了一种使用SoA的伪代码方法,以在从LOW到HIGH的范围内使用它,其中从LOW到HIGH + 1的增量可能被限制为偶数模60,以便使用模(和潜在的位打包只有2,3,5车轮上的条目)优化。”
这基于可能的实现,使用SoA二次方程组(4*x^2 + y^2)、(3*x^2 + y^2)和(3*x^2 -y^2),将它们表示为一系列数字,每个序列的x值固定在1到SQRT((HIGH-1)/4)、SQRT((HIGH-1)/3)之间的值,分别对 2*x^2 + 2*x - HIGH - 1 = 0 求解 x = (SQRT(1 + 2 * (HIGH + 1)) - 1) / 2。在我的F#代码中,根据顶部链接来表达这些序列。优化这些序列时,当仅筛选奇合数时,“4x”序列中的y值只需要是奇数,“3x”序列则需要在x为偶数和偶数为x时只使用奇数y值。进一步的优化通过观察上述序列的模式重复超小范围的x和y值来减少二次方程的解数(=序列中的元素数),Berstein代码中使用了这种方法,但在我的F#代码中尚未实现。

我也没有包括可以应用于素数“无平方因子筛法”的众所周知的优化,例如使用轮式分解和计算起始段地址的计算,这些在我的分段筛法实现中使用。

因此,为了计算“4x”、“3x+”和“3x-”(或者像我在F#代码中做的那样将“3x+”和“3x-”组合起来)的序列起始段地址,并根据上面的内容计算每个x的范围,伪代码如下:

  1. 计算范围LOW - FIRST_ELEMENT,其中FIRST_ELEMENT是每个给定的x值或y = x-1的最低适用的y值。

  2. 对于计算在此范围内有多少元素的任务,这归结为一个问题:有多少个(y1)^2 + (y2)^2 + (y3)^2...,其中每个y数字相差两个,以产生所需的偶数或奇数'y'。通常在平方序列分析中,我们观察到平方之间的差异具有恒定的增量,例如delta(9-1)为8,delta(25-9)为16,增加了8,delta(49-25)为24,进一步增加了8,等等。因此,对于n个元素,对于此示例,最后一个增量为8 * n。使用这个表达式来表示元素的序列,我们得到它是一个(或者选择作为第一个元素的任何一个)加上八倍于类似于(1 + 2 + 3 + ...+ n)的东西的序列。现在标准的线性序列约简适用于这个和为(n + 1) * n / 2或n^2/2 + n/2。我们可以通过解二次方程n^2/2 + n/2 - range = 0或n = (SQRT(8 * range + 1) - 1) / 2来解决在范围内有多少个n元素。

  3. 现在,如果FIRST_ELEMENT + 4 * (n + 1) * n不等于LOW作为起始地址,则将n加1,并使用FIRST_ELEMENT + 4 * (n + 2) * (n + 1)作为起始地址。如果使用进一步的优化将轮因子分解消除应用于序列模式,则可以使用查找表数组来查找满足条件的最接近的n值。

  4. 可以直接计算起始元素的模12或60,也可以通过基于模重复性质的查找表来产生。

  5. 然后,每个序列都用于切换到HIGH限制的组合状态。如果添加了额外的逻辑以在每个序列中仅跳转到适用的元素之间的值,则无需再使用模条件。

  6. 上述操作针对每个“4x”序列,接着是“3x+”和“3x-”序列(或者将“3x+”和“3x-”组合成一个“3x”序列),直到先前计算的x限制或循环测试的限制。

这就是它:如果使用适当的筛选范围分段方法,最好使用与CPU缓存大小相关的固定大小以实现最佳内存访问效率,则可以像Bernstein一样对SoA进行分段的方法,但表达方式要简单一些,因为它提到了但没有结合模数运算和位打包。

网页内容由stack overflow 提供, 点击上面的
可以查看英文原文,
原文链接