第一次關注到這個問題是在做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以下所有自然數的和減去所有合數的和，自然數的和可以直接用求和公式計算，問題在於篩選出所有合數並求和，顯然這里可以先用埃拉托斯特尼篩篩選出\(\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\)等等。則有：
其中\([\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\)的數中所有奇數的個數，則有：
這個公式同樣可以使用證明容斥原理時使用的數學歸納法加以證明，詳細的證明我就不寫了，只說明一下\(a=2\)的簡單情況：
有了這個遞推式和上面給出的邊界條件，我們就可以計算出\(\phi(x,a)\)。如果我們設\(a=\pi(\sqrt{x})\)，則\(\phi(x,a)\)實際上表示是\(\sqrt{x}\)到\(x\)之間素數的個數，則可以得到：
因而我們可據此算出\(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\)的數中所有奇數的和，則我們可以得到以下遞推式：
和上面類似，我們只說明一下\(a=2\)的簡單情況，定義\(\sigma(x)\)表示\([1,x]\)的自然數之和，如我們有\(\sigma(5)=1+2+3+4+5=15\)，則有：
如我們設\(\sum(x)\)表示小於等於\(x\)的所有素數的和，並設\(a=\pi(\sqrt{x})\)則根據和上面類似的原理，我們有：
據此我們可以求出小於等於\(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)\)表示\([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\)，據此我們展開上式有：
通過適當的選擇\(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})\)，則有：
經過不太復雜的推導，可以發現\(P_2(x,a)\)可通過下式計算：
其中\(a<\pi(\sqrt{x})=b\)，據此我們可以計算\(\pi(x)\)。使用和上面的類似的原理，並定義\(T_2(x,c)\)為區間\([1,x]\)中恰好有兩個素因子，且兩個素因子均大於\(p_c\)的數的和，則有：
其中\(c=\pi(x^{1/3})\)，且：
綜合上面兩個公式，我們也可以計算\(\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])。上面幾種算法的詳細的時間與空間復雜度分析可以參見這篇文章，這里列一個簡要的結果：
在此之后，Deleglise-Rivat以及Xavier Gourdon又對算法做出了改進。在github上作者Kim Walisch對以上算法做了實現，他的測試結果如下：
在2015年9月，他宣布利用自己的算法計算了\(\pi(10^{27})\)，計算花費了數年時間，最高峰時使用了235G的內存，之后他們又花了五個月時間進行重新計算驗證，確認的結果是：
這是目前[維基頁面](https://en.wikipedia.org/wiki/Prime-counting_function#Algorithms_for_evaluating_%CF%80(x)上列出的最大的\(\pi(x)\)值，至於是否有人算出了更大的\(\pi(x)\)值，我沒有查到確切的信息。
以上計算\(\pi(x)\)的算法中自Lehmer以后都變得非常復雜，至少我已經無法理解。同時由於使用的方法越來越復雜，其代碼實現也會變得更加麻煩。而且相關方法能否進行推廣用來計算N以下的質數和也未可知，這些問題感興趣的同學可以繼續探索。
最后，參考文獻[7]，設\(\pi(x)=n\)，則\(x\)以下質數的和的漸進估計是：
全文完。
六、參考文獻
- Crandall, R., & Pomerance, C. B. (2006). Prime numbers: a computational perspective (Vol. 182). Springer Science & Business Media.
- 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.
- Lagarias, J. C., Miller, V. S., & Odlyzko, A. M. (1985). Computing 𝜋 (𝑥): the Meissel-Lehmer method. Mathematics of Computation, 44(170), 537-560.
- Lagarias, J. C., & Odlyzko, A. M. (1987). Computing π (x): An analytic method. Journal of Algorithms, 8(2), 173-191.
- Lehmer, D. H. (1959). On the exact number of primes less than a given limit. Illinois Journal of Mathematics, 3(3), 381-388.
- Riesel, H. (2012). Prime numbers and computer methods for factorization (Vol. 126). Springer Science & Business Media.
- Sinha, N. K. (2010). On the asymptotic expansion of the sum of the first n primes. arXiv preprint arXiv:1011.1667.
- Xavier Gourdon, Computation of pi(x) : improvements to the Meissel, Lehmer, Lagarias, Miller, Odllyzko, Deléglise and Rivat method, February 15, 2001.