注:对知乎的公式编辑功能实在无力吐槽,用typora写的文章直接粘过来公式无法显示,只好又手工加上了全部公式,不过可能还是会有遗漏。大家可以点击这个链接 查看我的博客原文。以下是正文:
第一次关注到这个问题是在做project euler第10题的时候,原题目是要求两百万以内质数的和,知乎的题目把这个数字调到了10亿,事实证明这个规模调整是决定性的,很多在小规模可用的算法在10亿这个规模都不可用了。和其它欧拉工程的题目类似,这个题目存在一个很明显的暴力解法,但也存在一些效率更高的算法。暴力解法要不是通过对N以下的每个奇数做素性测试,要不是通过埃拉托斯特尼筛或者其它线性与亚线性筛得到N以下的所有质数然后相加,如菜鱼ftfish所言,这种暴力算法存在素数个数决定的时间复杂度下限,所以肯定还存在更优的算法。可以优化的根本原因在于,要计算N以下所有质数的和并不需要知道N以下的所有质数。在此基础上,我们可以使用各种技巧来提升算法的表现。
这个答案下目前最快的算法应该就是菜鱼ftfish所列的Lucy Hedgehog给出的算法,这个算法d在两个方面让人好奇,第一是效率极高,在时间和空间复杂度上的表现都极为优异,甚至对于python这种较慢的脚本语言计算十亿内的质数和都可以在一秒内出结果,而我自己用python写的的暴力算法甚至完全无法处理这个数据规模。第二是算法中采用了动态规划的思路,给出了一个求解质数和的递推式,让人非常好奇作者是怎么想到的,以及这个递推式背后有没有更为一般和深刻的原理。可惜的是这个代码可能由于过度优化导致可读性变得很差,原作者在论坛里也没有做详细的解释,因此在好奇心驱使下,我仔细做了一点研究,读了一些相关文献,对不同的算法做了尝试,最终实现了四个不同的算法版本。我想知道自己能不能写出一个更快的算法,因为我自己的主力语言也是python,和Hedgehog使用的语言相同,在同种语言下的比较应该是相对公平的。事实表明仅有一个算法在小规模上数据上的表现比Hedgehog的算法要更好,但在大规模数据上,Hedgehog的算法仍然是最稳健的。下面列的是我的四个算法和Hedgehog算法的对比:
图中横轴表示数据规模的对数,纵轴表示运行时间的对数。在我写的这四个算法中,表现最好的是hedgehog_recursive这个算法,使用带备忘录的自上而下动态规划方法实现了Hedgehog算法中的原理,奇怪的是这个算法虽然只是直接翻译了菜鱼ftfish的数学推导,却在小数据规模上表现到如此之好,基本上耗时都在Hedgehog算法的3%以内,这不点我也不是很理解。其次表现类似的是sum_primes_sieve和legendre这两个算法,前面这个算法使用了一个改进的埃拉托斯特尼筛,而后者则是对法国数学家勒让德提出的一个计算N以下素数个数的递推式的推广,使其可以计算N以下素数的和,我猜测这也是Hedgehog使用的递推式的灵感来源。表现再次的是meissel这个算法,它依据的是德国天文学家对勒让德计算N以下素数个数算法的改进,并将其推广到可以计算素数的和,理论上这个算法应比勒让德的算法更优,但实际算法表现并没有更好,可能是我的算法实现的原因。
下面我具体介绍一下以上四种算法的基本原理和代码实现,我相信在代码上还有很多优化的空间,大家如果有什么改进意见敬请提出来。
一、改进的埃拉托斯特尼筛
要求N以下的所有质数的和,一个显而易见的原理是用N以下所有自然数的和减去所有合数的和,自然数的和可以直接用求和公式计算,问题在于筛选出所有合数并求和,显然这里可以先用埃拉托斯特尼筛筛选出以内的所有质数,才依次筛选出这些质数在N以下的倍数并求和,这里的问题是有些倍数被重复计算了,可以用某些欧拉筛来避免重复筛选,或者也可以用python的集合来去重,我这里使用的是后者,因为经过尝试我发现用集合去重比用算法来避免重复筛选效率更高。
以上的算法明显还可以继续改进,首先想到所有除二以外的质数都为奇数,所以只需在奇数中筛选即可,再用所有奇数的和减去奇合数的和即可。进一步的,我们知道所有大于三的素数都可以表示成为
的形式,因此我们只需要列出所有
形式的数,用求和公式计算其总和,再筛选其中的合数并求和(同样用集合来去重),两者相减即为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的原始算法,这也是我感觉奇怪的地方。我猜测原因可能是在这里的动态规划中,有很多中间数据无需计算,自上而下的动态规划可以直接跳过这些数据,回到初始的边界条件,而自下而上的动态规划则必需一步步的计算,才能得到最终的计算结果。这个算法的最大问题在于处理大规模数据时递归深度过深的问题,根据这个算法,其递归树的深度约为