Atkin筛法 - 解释和Java示例

32

我在维基百科上读到了有关Atkin筛的内容,但目前维基百科的内容还比较有限。我在寻找一个高层次的Atkin筛解释,并附带一个Java示例。

谢谢。


这并不是你想要的,但至少是一个开端。 - keyser
可能是Sieve of Atkin explanation的重复问题。 - QuantumMechanic
@QuantumMechanic 好吧,他提到了维基百科,但没有Java示例。 - keyser
1个回答

134
您可能已经了解了这里介绍的一些关于质数、合数和筛法的基本概念,但这些内容对其他读者理解算法的性质会有所帮助。本回答中的一些内容危险地接近于属于数学版的StackOverflow,但我认为其中一些内容是必要的,可以建立算法的操作和实现之间的联系。
维基百科文章中的三个模运算和二次配对是从Atkin和Bernstein论文中得出的三个配对派生而来,该论文发表了这个筛法的核心定理(及其证明),并显示它们共同形成了一个质数筛法。任何一个单独使用都只能产生质数,但不能产生所有质数。需要同时使用这三个配对才能产生所有质数。
以下是一个实现维基百科算法的Java程序。我不保证实现效率,只是它是一个在Java中直接实现维基百科算法的工作程序。
// SieveOfAtkin.java

import java.util.Arrays;

public class SieveOfAtkin {
private static int limit = 1000;
private static boolean[] sieve = new boolean[limit + 1];
private static int limitSqrt = (int)Math.sqrt((double)limit);

public static void main(String[] args) {
    // there may be more efficient data structure
    // arrangements than this (there are!) but
    // this is the algorithm in Wikipedia
    // initialize results array
    Arrays.fill(sieve, false);
    // the sieve works only for integers > 3, so 
    // set these trivially to their proper values
    sieve[0] = false;
    sieve[1] = false;
    sieve[2] = true;
    sieve[3] = true;

    // loop through all possible integer values for x and y
    // up to the square root of the max prime for the sieve
    // we don't need any larger values for x or y since the
    // max value for x or y will be the square root of n
    // in the quadratics
    // the theorem showed that the quadratics will produce all
    // primes that also satisfy their wheel factorizations, so
    // we can produce the value of n from the quadratic first
    // and then filter n through the wheel quadratic 
    // there may be more efficient ways to do this, but this
    // is the design in the Wikipedia article
    // loop through all integers for x and y for calculating
    // the quadratics
    for (int x = 1; x <= limitSqrt; x++) {
        for (int y = 1; y <= limitSqrt; y++) {
            // first quadratic using m = 12 and r in R1 = {r : 1, 5}
            int n = (4 * x * x) + (y * y);
            if (n <= limit && (n % 12 == 1 || n % 12 == 5)) {
                sieve[n] = !sieve[n];
            }
            // second quadratic using m = 12 and r in R2 = {r : 7}
            n = (3 * x * x) + (y * y);
            if (n <= limit && (n % 12 == 7)) {
                sieve[n] = !sieve[n];
            }
            // third quadratic using m = 12 and r in R3 = {r : 11}
            n = (3 * x * x) - (y * y);
            if (x > y && n <= limit && (n % 12 == 11)) {
                sieve[n] = !sieve[n];
            } // end if
            // note that R1 union R2 union R3 is the set R
            // R = {r : 1, 5, 7, 11}
            // which is all values 0 < r < 12 where r is 
            // a relative prime of 12
            // Thus all primes become candidates
        } // end for
    } // end for
    // remove all perfect squares since the quadratic
    // wheel factorization filter removes only some of them
    for (int n = 5; n <= limitSqrt; n++) {
        if (sieve[n]) {
            int x = n * n;
            for (int i = x; i <= limit; i += x) {
                sieve[i] = false;
            } // end for
        } // end if
    } // end for
    // put the results to the System.out device
    // in 10x10 blocks
    for (int i = 0, j = 0; i <= limit; i++) {
        if (sieve[i]) {
            System.out.printf("%,8d", i);
            if (++j % 10 == 0) {
                System.out.println();
            } // end if
            if (j % 100 == 0) {
                System.out.println();
            } // end if
        } // end if
    } // end for
} // end main
} // end class SieveOfAtkin

我有一份Atkin(与Bernstein合著)原始论文的副本,其中描述了筛法构建的定理。该论文可在此处获取:http://www.ams.org/mcom/2004-73-246/S0025-5718-03-01501-1/S0025-5718-03-01501-1.pdf。对于非数学家而言,这是一篇密集的阅读材料,具有美国数学学会论文的简洁性。
以下是从描述和Atkin与Bernstein的论文中推导算法的更深入解释。
Atkin和Bernstein(理所当然地)假设其读者已充分理解素数筛法、模运算和使用模运算进行轮子分解。不幸的是,Wikipedia文章的描述和算法也假设了类似的事情,尽管程度稍微轻微些。Atkin和Bernstein并没有声称他们的三个轮子分解和不可约二次式是唯一可用的,并且只是简单地提供了其他可用对的例子,而不进一步解释如何使用它们。因此,Atkin和Bernstein给出定理和证明的三个对是基于他们工作的算法中使用的三个对。Atkin和Bernstein还没有声称他们的三个对是最优的,它们显然是方便的。
在此类讨论中非常有用的特殊数学符号在这里不可用。为了回答问题,我将使用以下表示:

{一些列举的集合或定义其中一个的属性}

表示一个集合;

Nat0

表示包括零在内的自然数集合,即Nat0 = {0, 1, 2, ...};

Nat

表示不包括零的自然数集合,即Nat = {1, 2, 3, ...},以及以下构造用于定义集合和其中一个元素的符号:

{集合中元素的符号:通过该符号定义集合的条件}

#

表示集合基数,即集合中元素的数量;

^

表示乘方,即x的平方写作x^2。
轮子分解中使用的模运算出现在两种等效形式中:

n = (k * m) + r,其中k属于Nat0,r属于R = {r: r属于Nat0且r

n mod m = r,其中r属于R = {r: r属于Nat0且r 以下是论文中定理给出的定义以及一些关于模算术形式的注释:

1. 当且仅当有奇数个(x, y)对满足二次方程n = (4 * x^2) + (y^2),其中x和y为大于等于1的整数,n mod 4 = 1且n没有完全平方因子时,#{(x, y) : n = (4 * x^2) + (y^2), n in {n : (Nat0 * 4) + 1}, where x and y >= 1 and n has no perfect square factor}为奇数,则n始终为质数。请注意,这里的形式n mod m = r,其中r在R中,m = 4且R = {r : 1}。
2. 当且仅当有奇数个(x, y)对满足二次方程n = (3 * x^2) + (y^2),其中x和y为大于等于1的整数,n mod 6 = 1且n没有完全平方因子时,#{(x, y) : n = (3 * x^2) + (y^2), n in {n : (Nat0 * 6) + 1}, where x and y >= 1 and n has no perfect square factor}为奇数,则n始终为质数。请注意,这里的形式n mod m = r,其中r在集合R中,m = 6且R = {r : 1}。
3. 当且仅当有奇数个(x, y)对满足二次方程n = (3 * x^2) - (y^2),其中x,y为整数,x > y >= 1,n mod 12 = 11且n没有完全平方因子时,#{(x, y) : n = (3 * x^2) - (y^2), {n : (Nat0 * 12) + 11}, x > y >= 1 and n has no perfect square factor}为奇数,则n始终为质数。请注意,这里的形式n mod m = r,其中r在集合R中,m = 12且R = {r : 11}。
论文和维基百科文章假设读者熟悉轮筛法的一个属性,即可以使用模算术来选择只有某些质因数的整数。在形式上,
n mod m = r,其中r在R中,R = {r : Nat0, r < m},
如果只选择与m互质的R元素,则满足表达式的所有整数n要么是质数,要么与m互质。
相对质数指它们没有大于1的公共整数因子。相对质数的例子有:2和3互质,4和9互质,9和14互质。理解这一点对于理解模算术中使用的剩余(余数)以及它们在各种版本的解释和算法中的等效性是至关重要的。
现在将解释定理、算法和解释之间的关系。
对于第一个二次方程,n =(4 * x ^ 2)+(y ^ 2):
该定理采用以下形式:
n =(k * 4)+ r,其中r属于R1 = {r:1},k属于Nat0
这与写作如下相同:
n mod 4 = r,其中r属于R1 = {r:1}
请注意,它定义了n必须是从1开始的每个奇数自然数,即{1、5、9、13、...}。
对于算法,可以选择各种m值,并且通过正确设置R集合,可以保留定理所示的属性。本文和维基百科文章的作者假设读者已经知道所有这些,并会立即认识到它。对于论文和维基百科文章中使用的其他m值,当等效项为以下内容时:
n mod 12 = r,其中r属于R1a = {r:1、5、9}
n mod 60 = r,其中r属于R1b = {r:1、5、9、13、17、21、25、29、33、37、41、45、49、53、57}
R1a和R1b集合中的某些元素可能因为后面要解释的两个原因而被删除,但定理仍然适用。
对于第二个二次方程,n =(3 * x ^ 2)+(y ^ 2):
该定理采用以下形式:
n =(k * 6)+ r,其中r属于R2 = {r:1},k属于Nat0
同样,这等于
n mod 6 = r,其中r属于R2 = {r:1}
请注意,在自然数中,这是从1开始的每第三个奇数,即{1、7、13、19、...}
文章中的等效项如下:

当n除以12的余数为r,其中r属于R2a = {r : 1, 7}时:

n mod 60 = r,其中r属于R2b = {r : 1, 7, 13, 19, 25, 31, 37, 43, 49, 55}。

同样,集合R2a和R2b中的值可以因为后面将要解释的两个原因而被删除,但定理依旧适用。

对于第三个二次方程,(3 * x^2) - (y^2):

该定理使用以下形式:

n = k * 12 + r,其中k为自然数0,r属于R3a = {r : 11}。

同样,这等同于:

n mod 12 = r,其中r属于R3a = {r : 11}。

请注意,这是从11开始的自然数奇数序列中每六个数字中的一个,即{11、23、35、47、...}。

在论文和文章中等效的n为:

n mod 60 = r,其中r属于R3b = {r : 11, 23, 35, 47, 59}。

同样,集合R3b中的值可以因为后面将要解释的原因而被删除,但定理依旧适用。

在各种算法和解释中,对于m = 12和m = 60的值,可以删除集合R的元素,而不影响定理或算法的有效性。一些值在集合R中被抛弃的原因有两个。

第一个原因是,如果集合R中的任何r值与其配对的m不互质,则只会包括具有一个或多个m的素因子的复合整数n,其中没有一个是素数。这种模运算的特点是,轮筛法可用于从进一步检测是否是素数的大量非素数中过滤出来,通常需要更复杂且效率较低的测试。在此筛法中,更复杂的测试是特定不可约二次方程的解数是否为奇数。这意味着我们可以立即丢弃该算法中集合R中与使用该集合R的m值不互质的所有值。

第二个原因是,在文献中,轮子分解会创建重叠的整数集,包括重叠素数。虽然它们很方便,而且在定理中重叠并不重要,但在算法中,如果能够轻松避免重叠,则是浪费的。在这种情况下,可以轻松地避免。另外,如果来自轮子分解的整数集重叠,则一个二次方程中的奇数解的数量加上另一个二次方程中的奇数解的数量将变成累计偶数解(奇数加奇数始终是偶数)。在许多实现中,包括维基百科实现,这将确定质数不是质数,因为类似维基百科这样的实现从所有二次方程的累积解来确定质性,并且不将每个二次方程的解隔离。在这些情况下,从轮子分解中获取的整数必须是互斥的子集,这是至关重要的。
实现不希望在不必要的情况下在多个二次方程中测试相同的数字,也没有必要。对于在三个二次方程中使用的一组R中的r值,假设使用相同的m,需要仅出现在其中一个中。如果它出现在多个位置,则n的相同值将显示多次,并使用多个二次方程进行测试。对于正在使用的相同m值,确保R的相同元素不会出现在多个二次方程的R中将消除重叠。在维基百科文章的情况下,轮子分解所防止的重叠防止了累积二次方程解,在单个二次方程中为奇数,但在两个二次方程中积累到偶数。
不同的算法可能会在计算二次方程之前避免重叠。在对二次方程和轮子分解进行经验测试时,其中m = 12的轮子分解比解决二次方程的n值要少得多。使用m = 60的轮子分解将显着增加差异。如果针对特定n值的二次方程解算法非常高效,则仅使用来自轮子分解的值来测试二次方程将带来显着改进。
以下是删除不相对质的元素后的轮子分解。对于第一个二次方程:
n mod 12 = r,其中r在R1a = {1:1,5}(9具有与12公共的因子3)
n mod 60 = r,其中r在R1b = {r:1,13,17,29,37,41,49,53}(5、25和45具有与60公共的因子5;9、21、33、45和57与60具有因子3),这是Atkin和Bernstein论文中算法的形式。
对于第二个二次方程:
第一类二次剩余方程:
当m=12时,R1a={1, 5} (R中元素与12没有公共因子)
当m=60时,R1b={1, 13, 17, 29, 37, 41, 49, 53} (25和55的公共因子为5)
第二类二次剩余方程:
当m=12时,R2a={7} (R中元素与12没有公共因子,1已在R1a中出现)
当m=60时,R2b={7, 19, 31, 43} (1、13、37和49已在R1b中出现,它们被删除了)
第三类二次剩余方程:
当m=12时,R3a={11} (R中元素与12没有公共因子)
当m=60时,R3b={11, 23, 47, 59} (35的公共因子为5)
一个仍未解决的问题是为什么m的值范围为4、6、12和60。这与要排除多少复合数(即非质数)参与更复杂的使用二次方程进行素数测试以及使用轮筛法的复杂度有关。
所使用的m值可以确定哪些复合数可以立即被排除而不会排除质数。如果m = 4且R1 = {r:1},如第一个二次方程定理中所述,所有具有2个质因数的数字都会被排除,因为1与所有数字互质,而4具有2个质因数。需要注意的是,由于3不在此集合R中,使用m = 4和集合R1的轮筛法也将排除大量的质数,可能是其中一半。
如果m = 6且R2 = {r:1},如第二个二次方程定理中所述,所有具有2或3个质因数的复合数都会被排除,因为1与所有数字互质,而6具有2和3个质因数。同样,使用m = 6和集合R2,不包含5,也将排除大量的质数,可能是其中一半。
如果m = 12且R3 = {r:11},如第三个二次方程定理中所述,所有具有2或3个质因数的复合数都会被排除,因为11与12互质且12具有2和3个质因数。同样,使用m = 12和集合R3,不包含1、5或7,也将排除大量的质数,可能超过其中一半。
Atkin和Bernstein在其论文中非正式地展示了一个事实,即虽然定理中的轮筛法分别从它们各自的二次方程中排除质数,但是这三个轮筛法共同允许所有质数,并且在第一和第二个二次方程的轮筛法中,存在显着的重叠。尽管他们的算法没有在m = 60时消除重叠,但维基百科文章中在算法中设置m = 12,在文章解释中设置m = 60。
对于Atkin和Bernstein在他们的定理中使用的二次型,如果削弱与其相关的轮因式分解,则会使它们不符合定理操作。然而,我们可以通过某些方式加强它们,仅删除更多组合数,但仍保留完全相同的质数。对于其中m = 4(4 = 2 * 2)的形式,每个偶数都会被过滤。对于其中m = 12(12 = 2 * 2 * 3)的形式,每个具有2或3为质因数的整数都会被过滤。对于其中m = 60(60 = 2 * 2 * 3 * 5)的形式,每个具有2、3或5为质因数的整数都会被过滤。我们可以使用m = 6的过滤器实现与m = 12相同的效果,使用m = 30实现与m = 60相同的效果,但必须小心创建的内容是否与定理中使用的内容等效。
以下是关于轮因式分解的一些有用统计信息。自然数中50%的整数是偶数,除了2之外都不是质数。自然数中33%的整数具有3作为质因数且不是质数。自然数中20%的整数具有5作为质因数且不是质数。总体而言,自然数中67%的整数具有2或3的质因数且不是质数。总体而言,约75%的自然数具有2、3或5的质因数且不是质数。消除Nat0中1/2、2/3或3/4的复合数字以进行更复杂的素数测试的简单方法非常诱人,并且使用轮因式分解作为素数筛选器的初步过滤器的动机。这也是使用具有伴随集R的m值的动机,该集合可以过滤掉代表大量组合数的所有复合数。
作为绝对最小值,我们希望删除具有2(即所有偶数)的质因子的复合数,然后在最后添加2。我们至少希望删除具有2或3的质因子的复合数。最好,我们希望删除具有2、3或5个质因子的复合数。超出此范围,统计数据显示收益递减。 m = 4与R1实现了最低限度要求。 m = 12与R1a、R2a和R3a实现了最少所需的要求。 m = 60与R1b、R2b和R3b实现了非常理想的结果。
在使用m和集合R的值时还有一些需要考虑的事情。请注意,对于第一个二次型,这两种形式不等价:
n mod 12 = r,其中r在R1a = {r:1,5}中

n mod 6 = r,其中r在R = {r:1,5}中
因为当m=6时的形式与n mod 4 = r(其中r属于R1={r:1})不等价。注意到,n mod 6 = r(其中r属于R={r:1})等价于n mod 12 = r(其中r属于R={r:1,7}),且当第二个平方剩余被使用时,可以使用m=6的形式。实际上,在定理中使用了第二个平方剩余的形式。若我们使用该形式而非已考虑的例子,则在第一个平方剩余的m=12时,可将元素1从集合R中移除以消除重叠。
在调整算法时,必须谨慎以保持定理所要求的条件。在选择m的值和R集合时,必须确保形式是等价的,并且不会向平方剩余引入任何由定理中轮筛法未生成的新值。
对于实现,根据数据结构、算术、处理器能力(特别是乘除法)、可用处理器缓存、内存和数据结构大小等复杂性和效率进行选择。在需要检查R集合中的剩余数的值和数量之间存在权衡。一些原因可能导致解释中出现m=60和算法中出现m=12。Atkin和Bernstein在其论文中使用m=60的形式,这明显是他们使用算法的原因。
在Atkin和Bernstein的论文中,他们还提供了用格子找到特定n值的平方剩余解的算法。这些额外的算法使Atkin和Bernstein能够创建同时对平方剩余和轮筛法进行筛选的筛法算法。在维基百科文章的算法中,没有考虑用轮筛法得来的平方剩余的格子导出解的任何算法。在维基百科文章中,使用了一种穷举x、y值技术与平方剩余,并在平方剩余返回n值后应用轮筛法。同样,这是实现者必须决定的效率问题。

1
尽管它遵循了旧的WP文章(已经改进)中的SoA伪代码,但你从(旧的)WP文章中提取的Java代码并不是真正的SoA,因为你将x和y运行到平方根,而它们需要在每个三次方程中的不同点上提前终止,以便具有线性O(n)渐近计算复杂度和范围。我也在我的另一个SoA问题的答案中讨论了这个算法。 - GordonBGood
2
在你批评别人和他们的回答之前,请先阅读问题,然后再看答案。我没有声称我的Java实现是Bernstein论文中给出的SoA算法的实现。我已经从他们的论文中实现了格点算法,并将其放在一个子版本库中。它很复杂。这里的问题是关于WP算法的,就像在问题被提出时那样存在。根据Stack Overflow的惯例,这里的答案直接回答了这里的问题。 - Jim
此外,在告诉别人他们对答案做出了虚假的声明之前,请确实阅读答案。我没有声称这个特定的实现是Bernstein论文中呈现的算法,甚至也不高效。相反,我的回答第3段非常明确地说明了这一点。如果你的评论目的是为了膨胀自己的自尊心,那么从语气上看,我相信你可能已经实现了这一点。 - Jim
最后,关于算法效率(例如减少切换),它在计算机中可能不如一个更简单的算法那样高效,后者具有增加的切换。伯恩斯坦论文中晶格算法的计算成本可能超过了明显更高的数学效率。简而言之,在真实的编程语言和真实的机器上,一种带有一些冗余操作的简单算法可能比一种没有冗余的复杂算法更有效。我当时比较了伯恩斯坦晶格解法和WP解法,但我不记得结果了。 - Jim
2
@Jim,没有必要对一些建设性的批评如此防御。 - wvdz
显示剩余6条评论

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