计算特定数字以下质数和的几种算法

第一次关注到这个问题是在做project euler第10题的时候,原题目是要求两百万以内质数的和,知乎的题目把这个数字调到了10亿,事实证明这个规模调整是决定性的,很多在小规模可用的算法在10亿这个规模都不可用了。和其它欧拉工程的题目类似,这个题目存在一个很明显的暴力解法,但也存在一些效率更高的算法。暴力解法要不是通过对N以下的每个奇数做素性测试,要不是通过埃拉托斯特尼筛或者其它线性与亚线性筛得到N以下的所有质数然后相加,如菜鱼ftfish所言,这种暴力算法存在素数个数决定的时间复杂度下限,所以肯定还存在更优的算法。可以优化的根本原因在于,要计算N以下所有质数的和并不需要知道N以下的所有质数。在此基础上,我们可以使用各种技巧来提升算法的表现。

这个答案下目前最快的算法应该就是菜鱼ftfish所列的Lucy Hedgehog给出的算法,这个算法d在两个方面让人好奇,第一是效率极高,在时间和空间复杂度上的表现都极为优异,甚至对于python这种较慢的脚本语言计算十亿内的质数和都可以在一秒内出结果,而我自己用python写的的暴力算法甚至完全无法处理这个数据规模。第二是算法中采用了动态规划的思路,给出了一个求解质数和的递推式,让人非常好奇作者是怎么想到的,以及这个递推式背后有没有更为一般和深刻的原理。可惜的是这个代码可能由于过度优化导致可读性变得很差,原作者在论坛里也没有做详细的解释,因此在好奇心驱使下,我仔细做了一点研究,读了一些相关文献,对不同的算法做了尝试,最终实现了四个不同的算法版本。我想知道自己能不能写出一个更快的算法,因为我自己的主力语言也是python,和Hedgehog使用的语言相同,在同种语言下的比较应该是相对公平的。事实表明仅有一个算法在小规模上数据上的表现比Hedgehog的算法要更好,但在大规模数据上,Hedgehog的算法仍然是最稳健的。下面列的是我的四个算法和Hedgehog算法的对比:

img

图中横轴表示数据规模的对数,纵轴表示运行时间的对数。在我写的这四个算法中,表现最好的是hedgehog_recursive这个算法,使用带备忘录的自上而下动态规划方法实现了Hedgehog算法中的原理,奇怪的是这个算法虽然只是直接翻译了菜鱼ftfish的数学推导,却在小数据规模上表现到如此之好,基本上耗时都在Hedgehog算法的3%以内,这不点我也不是很理解。其次表现类似的是sum_primes_sieve和legendre这两个算法,前面这个算法使用了一个改进的埃拉托斯特尼筛,而后者则是对法国数学家勒让德提出的一个计算N以下素数个数的递推式的推广,使其可以计算N以下素数的和,我猜测这也是Hedgehog使用的递推式的灵感来源。表现再次的是meissel这个算法,它依据的是德国天文学家对勒让德计算N以下素数个数算法的改进,并将其推广到可以计算素数的和,理论上这个算法应比勒让德的算法更优,但实际算法表现并没有更好,可能是我的算法实现的原因。

下面我具体介绍一下以上四种算法的基本原理和代码实现,我相信在代码上还有很多优化的空间,大家如果有什么改进意见敬请提出来。

一、改进的埃拉托斯特尼筛

要求N以下的所有质数的和,一个显而易见的原理是用N以下所有自然数的和减去所有合数的和,自然数的和可以直接用求和公式计算,问题在于筛选出所有合数并求和,显然这里可以先用埃拉托斯特尼筛筛选出\(\sqrt{N}\)以内的所有质数,才依次筛选出这些质数在N以下的倍数并求和,这里的问题是有些倍数被重复计算了,可以用某些欧拉筛来避免重复筛选,或者也可以用python的集合来去重,我这里使用的是后者,因为经过尝试我发现用集合去重比用算法来避免重复筛选效率更高。

以上的算法明显还可以继续改进,首先想到所有除二以外的质数都为奇数,所以只需在奇数中筛选即可,再用所有奇数的和减去奇合数的和即可。进一步的,我们知道所有大于三的素数都可以表示成为\(6k\pm1\)的形式,因此我们只需要列出所有\(6k\pm1\)形式的数,用求和公式计算其总和,再筛选其中的合数并求和(同样用集合来去重),两者相减即为N 以下所有质数的和。算法的代码实现比较简单,我这里就不多做解释了。

from sympy import primerange
from math import floor

def sum_primes_sieve(n=2e6):
    primes = list(primerange(2,n**0.5+1))
    j = floor(n/6)
    total_sum = 6*j*(j+1)
    res_set = set()
    for p in primes[2:]:
        k = p
        while p*k < n:
            res_set.add(p*k)
            k += 4 if k%6==1 else 2
    ans = total_sum - sum(res_set)
    return ans+5

二、Hedgehog算法的递归版本

菜鱼ftfish解释了Hedgehog算法的基本原理,其核心是答案中所列的递推式,使得我们可以用动态规划来解决质数和的问题。Hedgehog算法中显然使用的是自下而上的动态规划,我好奇用自上而下的动态规划会如何表现。因此,我使用递归函数直接翻译了菜鱼ftfish对这个算法的解释,并使用python中functools模块中的lru_cache装饰器,实现算法的记忆化,这里免去自己写备忘录代码的麻烦。这个算法是我所花时间最短,但却是表现最好的算法,甚至在小规模数据上要好于Hedgehog的原始算法,这也是我感觉奇怪的地方。我猜测原因可能是在这里的动态规划中,有很多中间数据无需计算,自上而下的动态规划可以直接跳过这些数据,回到初始的边界条件,而自下而上的动态规划则必需一步步的计算,才能得到最终的计算结果。这个算法的最大问题在于处理大规模数据时递归深度过深的问题,根据这个算法,其递归树的深度约为\(\sqrt{N}\),如果要计算十亿内质数的和,则递归深度要达到31622层,而在我的python版本允许的最大递归深度仅为3000层,虽然可以自己修改python允许的最大递归深度,但python仍然会报“超过最大递归深度”的错误,所以这个算法能够处理的最大问题规模大约就在两百万左右,再大就无法保证正确执行了。可能可以使用尾递归的方法来解决这个问题,但python对尾递归优化的支持并不好,我并没有尝试,如果有人尝试成功了,可以分享一下经验。

from functools import lru_cache
from sympy import isprime
from math import floor

@lru_cache(maxsize=32)
def s(v,p):
    if v == 1:
        return 0
    if v == 2:
        return 2
    if p == 1:
        return (2+v)*(v-1)/2
    if p**2<=v and isprime(p):
        return s(v,p-1)-p*(s(floor(v/p),p-1)-s(p-1,p-1))
    else:
        return s(v,p-1)

def hedgehog_recursive(n=2e6):
    p = int(n**0.5) + 1
    return s(n,p)

三、勒让德算法

计算N以下所有素数的和似乎在数论领域并不是一个重要的问题,我看了一些文献,只在部分文献里看过对这个质数和的渐进估计,但并没有看到给出确切的质数和的值的算法分析。但是计算N以下的所有素数的个数的问题则是数论中的热门话题了,相关文献连篇累牍,因为这个问题在数论领域有相当的重要性,甚至还有一个专门的函数\(\pi(x)\)表示小于等于\(x\)的所有素数的个数。高斯和勒让德通过经验统计的方式猜测\(\pi(x)\approx x/ln(x)\),这个猜想在1896年得到证明成为素数定理,这是解析数论领域的最重要的成就之一。黎曼猜想也是在改进对\(\pi(x)\)的估计中被提出来的,现在应该是数论领域最重要的未被证明的猜想。除开解析数论的进路以外,很多数学家也在不断改进计算\(\pi(x)\)的确切值的算法,相关的研究进展大家可以参见这个维基页面,更深入的研究可以参见我在文末列出的参考文献。

我们对计算N以下素数和算法来源于对勒让德对计算N以下素数个数的算法的推广。在勒让德以前,数学家们计算\(\pi(x)\)的方法就是筛选出\(x\)以下的所有素数然后数个数,勒让德首次指出,为了计算\(\pi(x)\),我们并不需要知道\(x\)以下的所有素数,只需要知道\(\pi(\sqrt{x})\)就可以了。它给出的算法基于容斥原理。我们设\(\phi(x,a)\)表示小于等于\(x\)的数中不能被\(p_1,p_2\cdots p_a\)整除的数的个数,其中\(p_1,p_2\cdots p_a\)表示前\(a\)个质数,如\(p_1=2,p_2=3,p_3=5\)等等。则有:

$$ \phi(x,a)=[x]-\sum_{i=1}^{a}[\frac{x}{p_i}]+\sum_{i<j}^a[\frac{x}{p_ip_j}]-\sum_{i<j<k}^a[\frac{x}{p_ip_jp_k}]+\cdots $$

其中\([\cdots]\)表示下取整。这个公式的意思是为了计算小于等于\(x\)不能被\(p_1,p_2\cdots p_a\)整除的数的个数,我们从\(x\)中减去可以被\(p_1,p_2\cdots p_a\)整除的数的个数,但是这样同时被两个质数整除的数就重复减去了,所以我们需要把它们的个数加回来,但是这样又会导致对可以同时被三个质数整除重复加入了,所以我们需要减去这样的数的个数,之后依次类推,显然这只是容斥原理的一个简单应用。如果真的要用这个公式来计算\(\phi(x,a)\)仍然显得比较复杂,通过仔细分析\(\phi(x,a)\)的算法,我们可以发现一个递推式。我们定义\(\phi(x,0)=0\),而\(\phi(x,1)\)显然表示小于等于\(x\)的数中所有奇数的个数,则有:

$$ \phi(x,a)=\phi(x,a-1)-\phi(x/p_a,a-1) $$

这个公式同样可以使用证明容斥原理时使用的数学归纳法加以证明,详细的证明我就不写了,只说明一下\(a=2\)的简单情况:

$$ \begin{aligned} \phi(x,2)&=[x]-([\frac{x}{p_1}]+[\frac{x}{p_2}])+[\frac{x}{p_1p_2}]\\ \phi(x,1)&=[x]-[\frac{x}{p_1}]\\ \phi(x/p_2,1)&=[\frac{x}{p_2}]-[\frac{x}{p_{1}p_2}]\\ \Rightarrow \phi(x,a)-\phi(x/p_2,1)&=[x]-([\frac{x}{p_1}]+[\frac{x}{p_2}])+[\frac{x}{p_1p_2}]=\phi(x,2) \end{aligned} $$

有了这个递推式和上面给出的边界条件,我们就可以计算出\(\phi(x,a)\)。如果我们设\(a=\pi(\sqrt{x})\),则\(\phi(x,a)\)实际上表示是\(\sqrt{x}\)\(x\)之间素数的个数,则可以得到:

$$ \pi(x)=\phi(x,a)+\pi(\sqrt{x})-1 $$

因而我们可据此算出\(x\)以下的素数个数。可以看出,勒让德的算法将\(x\)以下所有质数分成了两部分,然后分别计算它们的个数,加起来即为\(x\)以下所有质数的个数。我们计算N以下质数和的算法也是基于同样的原理,我们设\(S(x,a)\)表示小于等于\(x\)的数中不能被\(p_1,p_2\cdots p_a\)整除的数的和,其中\(p_1,p_2\cdots p_a\)表示前\(a\)个质数,则\(S(x,1)\)显然表示小于等于\(x\)的数中所有奇数的和,则我们可以得到以下递推式:

$$ S(x,a)=S(x,a-1)-p_a\times S(\frac{x}{p_a},a-1) $$

和上面类似,我们只说明一下\(a=2\)的简单情况,定义\(\sigma(x)\)表示\([1,x]\)的自然数之和,如我们有\(\sigma(5)=1+2+3+4+5=15\),则有:

$$ \begin{aligned} S(x,2)&=\sigma(x)-(p_1\sigma([x/p_1])+p_2\sigma([x/p_2]))+p_1p_2\sigma(x/p_1p_2)\\ S(x,1)&=\sigma(x)-p_1\sigma([x/p_1])\\ S(x/p_2,1)&=\sigma([x/p_2])-p_1\sigma([x/p_1p_2])\\ \Rightarrow p_2S(x/p2,1)&=p_2\sigma([x/p_2])-p_1p_2\sigma([x/p_1p_2])\\ \Rightarrow S(x,2)&=s(x,1)-p_2S(x/p_2,1) \end{aligned} $$

如我们设\(\sum(x)\)表示小于等于\(x\)的所有素数的和,并设\(a=\pi(\sqrt{x})\)则根据和上面类似的原理,我们有:

$$ \sum(x)=S(x,a)+\sum(\sqrt{x})-1 $$

据此我们可以求出小于等于\(x\)的所有质数的和。可以看到这里的递推式和Hedgehog算法的递推式非常相似,区别在于这里的递推式中\(p\)表示素数,则不需要像Hedgehog算法那样需要判断是否为素数。更重要的区别在于如果使用递归实现Hedgehog算法,则其递归深度为\(\sqrt{x}\),而在这里的算法中,递归深度为\(\pi(\sqrt{x})\),前者明显大于后者,因而这里的算法可以更快的达到边界条件,并且因为素数越大则分布密度越低,则两者在大规模数据中差距会更加明显。如当\(x=10000\)时,\(\sqrt{x}=100\),而\(\pi(\sqrt{x})=\pi(100)=25\),前者的递归深度是后者的四倍。

勒让德算法的缺陷和第二个算法类似,都会在大规模数据上因为迭代深度的问题而无法计算,虽然勒让德算法已经大大减少了递归函数的递归深度,但减少的仍然不够。如要计算十亿以内的素数和,则勒让德算法的递归深度为\(\pi(\sqrt{10^9}=\pi(31623)=3401\),仍然超过python允许的最大递归深度。因此我们需要进一步缩减递归深度,meissel算法就是一个有趣的深度。

勒让德算法的代码实现如下:

from sympy import primerange,isprime
from functools import lru_cache

def legendre(n=2e6):
    primes = list(primerange(1,int(n**0.5)+1))

    @lru_cache(maxsize=8192)
    def sigma(x,a):
        if a == 1:
            return (x//2)**2 if x%2==0 else (x//2+1)**2
        elif x <= primes[a-1]:
            return 1
        else:
            return sigma(x,a-1) - primes[a-1]*sigma(x//primes[a-1],a-1)

    a = len(primes)
    res = sigma(n,a) + sum(primes) - 1
    return res

四、meissel算法

19世纪晚期,德国天文学家E. Meissel以上提到的勒让德算法进行了改进,进一步提升了计算\(\pi(x)\)的效率,他使用自己改进的算法计算了\(\pi(10^9)\),虽然比正确值小了56,不过考虑到他完全依靠手工计算,这个准确度已经非常惊人了。Meissel对勒让德算法的主要改进是加入了一个新项\(P_2(x,a)\),从而使得算法的时间复杂度从勒让德算法的\(O(x)\)改进到\(O(x/log^{3}x)\),空间复杂度从\(O(\sqrt{x})\)改进至\(O(\sqrt{x}/log\ x)\)。这里我们对Meissel的算法做了推广,使其可以计算N以下的质数和。首先我们对Meissel的算法做一个简单介绍:

定义\(P_k(x,a)\)表示\([1,x]\)中恰好拥有\(k\)个素因子,且这\(k\)个素因子\(p_i,p_j\cdots p_k\)均大于\(p_a\)的数的个数,则我们有:

$$ \phi(x,a)=\sum_{k=0}^{+\infty}P_k(x,a) $$

这个公式我们可以这么理解:\(\phi(x,a)\)表示\([1,x]\)中其最小素因子都大于\(p_a\)的数的个数,这些数里面包括只有一个大于\(p_a\)的素因子的数,实际上也就是大于\(p_a\)的素数;也包括恰好有两个素素因子且这两个素因子都大于\(p_a\),也包括恰好有三个素素因子且这三个素因子都大于\(p_a\),依次类推以至无穷就可以得到所有其最小素因子都大于\(p_a\)的数的个数,也就是\(\phi(x,a)\)。我们定义\(P_0(x,a)=1\),而\(P_1(x,a)\)表示\([1,x]\)中大于\(p_a\)的素数的个数,则\(P_1(x,a)=\pi(x)-a\),据此我们展开上式有:

$$ \phi(x,a)=1+\pi(x)-a+P_2(x,a)+P_3(x,a)+\cdots $$

通过适当的选择\(a\)的数值,我们可以让\(P_k(x,a)\)及以后的项变为零。如我们选择\(a=\pi(\sqrt{x})\),则有\(p_{a+1}>\sqrt{x}\),此时\(P_2(x,a)\)及以后的变为零,上式变成之前提到的勒让德的公式,证明勒让德公式只是这个公式的一个特例。如我们选择\(a=\pi(x^{1/3})\),则有\(x^{1/3}<p_{a+1}<x^{1/2}\),此时\(P_3(x,a)\)及以后的项变为零。一般地,选择\(a=\pi(x^{1/r})\),则有\(P_r(x,a)=0\)\(P_{r-1}(x,a)\ne0\)。假设我们选择选择\(a=\pi(x^{1/3})\),则有:

$$ \phi(x,a)=1+\pi(x)-a+P_2(x,a)\Rightarrow \pi(x)=\phi(x,a)+a-1-P_2(x,a) $$

经过不太复杂的推导,可以发现\(P_2(x,a)\)可通过下式计算:

$$ P_2(x,a)=\sum_{i=a+1}^{b}\bigg\{\pi\bigg(\frac{x}{p_i}\bigg)-(i-1)\bigg \}=-\frac{(b-a)(b+a-1)}{2}+\sum_{i=a+1}^{b}\pi\bigg(\frac{x}{p_i}\bigg) $$

其中\(a<\pi(\sqrt{x})=b\),据此我们可以计算\(\pi(x)\)。使用和上面的类似的原理,并定义\(T_2(x,c)\)为区间\([1,x]\)中恰好有两个素因子,且两个素因子均大于\(p_c\)的数的和,则有:

$$ \sum(x)=S(x,c)+\sum(x^{1/3})-T_2(x,c)-1 $$

其中\(c=\pi(x^{1/3})\),且:

$$ T_2(x,c)=\sum_{i=c+1}^{b}\bigg\{p_{i}\times\bigg[\sum(x/p_i)-\sum_{j=1}^{j=i-1}p_j\bigg]\bigg\} $$

综合上面两个公式,我们也可以计算\(\sum(x)\)。这个算法相对于勒让得算法的最显著的优势是需要递归的次数更少,如当\(x=10^9\)时,\(c=\pi((10^9)^{1/3})=\pi(10^3)=168\),因此最大递归深度只有168层。但在我自己实现这个算法时,它的表现并没有比勒让得算法更优,我猜测原因是计算\(T_2(x,c)\)时消耗了过多资源,不过也有可能是我自己算法实现的原因。如果大家有提升这个算法的方法,还请指出。这个算法的代码实现如下:

def meissel(n=2e6):
    primes = list(primerange(1,int(n**(1/2))+1))
    cr_primes = [x for x in primes if x<n**(1/3)]

    def prime_sum(n):
        ps = list(primerange(1,n+1))
        return sum(ps)

    @lru_cache(maxsize=32)
    def sigma(x,a):
        if a == 1:
            return (x//2)**2 if x%2==0 else (x//2+1)**2
        else:
            return sigma(x,a-1) - primes[a-1]*sigma(x//primes[a-1],a-1) 

    def total(x,a):
        res,index = 0,a
        for p in primes[a:]:
            res += p * (prime_sum(x//p) - sum(primes[:index]))
            index += 1
        return res

    a = len(cr_primes)
    ans = sigma(n,a) + sum(cr_primes) - total(n,a) - 1
    return ans

经过尝试,至少到我写这篇文章为止,我自己并没有能够实现一个在效率表现上和处理大规模问题上比Hedgehog算法更优的算法。但这种尝试仍是有意义的,至少可以让我更加清楚的理解Hedgehog算法的原理,并且对质数计数的各种算法和推广加深了了解。我相信这些算法仍有优化的空间,如果找到的话我再来做补充。

五、延伸阅读

前面已经提到,质数计数是数论中一个非常重要的问题,既有使用解析数论等方法对质数整体分布规律的研究,也有各种不断改进的计算\(\pi(x)\)的确切值的算法。在以上我提到的Meissel算法之后约半个世纪,Lehmer[5]对这个算法进行了进一步改进,他进一步计算了\(P_3(x,a)\),并使用IBM 701计算机计算了\(\pi(10^{10})\),他的计算结果只比正确结果大一,以上从勒让德到Lehmer的算法实质都是对勒让德最初提出算法的某种改进,统称为计算\(\pi(x)\)的组合方法(Combinatorial method)。之后在1987年,Lagarias and Odlyzko[4]提出了一种从从解析数论角度计算\(\pi(x)\)的方法,算法中需要用到黎曼Zeta函数的某些性质,并采用数值积分的方法,这种方法被称为计算\(\pi(x)\)的分析方法(见参考文献[1]),这个算法虽然具备更好的渐近性,它在性能上比不上作者们在1985年提出的对Meissel-Lehmer算法的改进(见[3])。上面几种算法的详细的时间与空间复杂度分析可以参见这篇文章,这里列一个简要的结果:

img

在此之后,Deleglise-Rivat以及Xavier Gourdon又对算法做出了改进。在github上作者Kim Walisch对以上算法做了实现,他的测试结果如下:

img

在2015年9月,他宣布利用自己的算法计算了\(\pi(10^{27})\),计算花费了数年时间,最高峰时使用了235G的内存,之后他们又花了五个月时间进行重新计算验证,确认的结果是:

$$ \pi(10^{27})=16,352,460,426,841,680,446,427,399 $$

这是目前维基页面上列出的最大的\(\pi(x)\)值,至于是否有人算出了更大的\(\pi(x)\)值,我没有查到确切的信息。

以上计算\(\pi(x)\)的算法中自Lehmer以后都变得非常复杂,至少我已经无法理解。同时由于使用的方法越来越复杂,其代码实现也会变得更加麻烦。而且相关方法能否进行推广用来计算N以下的质数和也未可知,这些问题感兴趣的同学可以继续探索。

最后,参考文献[7],设\(\pi(x)=n\),则\(x\)以下质数的和的渐进估计是:

$$ \sum_{p<x}p_i=\frac{n^2}{2}(ln\ n+ln\ ln\ n-3/2+o(1)) $$

全文完。

六、参考文献

  1. Crandall, R., & Pomerance, C. B.(2006). Prime numbers: a computational perspective (Vol. 182). Springer Science & Business Media.
  2. Deléglise, M., & Rivat, J. (1996). Computing 𝜋 (𝑥): the Meissel, Lehmer, Lagarias, Miller, Odlyzko method. Mathematics of Computation of the American Mathematical Society, 65(213), 235-245.
  3. Lagarias, J. C., Miller, V. S., & Odlyzko, A. M.(1985). Computing 𝜋 (𝑥): the Meissel-Lehmer method. Mathematics of Computation, 44(170), 537-560.
  4. Lagarias, J. C., & Odlyzko, A. M.(1987). Computing π (x): An analytic method. Journal of Algorithms, 8(2), 173-191.
  5. Lehmer, D. H.(1959). On the exact number of primes less than a given limit. Illinois Journal of Mathematics, 3(3), 381-388.
  6. Riesel, H. (2012). Prime numbers and computer methods for factorization (Vol. 126). Springer Science & Business Media.
  7. Sinha, N. K.(2010). On the asymptotic expansion of the sum of the first n primes. arXiv preprint arXiv:1011.1667.
  8. Xavier Gourdon, Computation of pi(x) : improvements to the Meissel, Lehmer, Lagarias, Miller, Odllyzko, Deléglise and Rivat method, February 15, 2001.