java循环1000000000_求十亿内所有质数的和，怎么做最快?
注：對知乎的公式編輯功能實在無力吐槽，用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的原始算法，這也是我感覺奇怪的地方。我猜測原因可能是在這里的動態規劃中，有很多中間數據無需計算，自上而下的動態規劃可以直接跳過這些數據，回到初始的邊界條件，而自下而上的動態規劃則必需一步步的計算，才能得到最終的計算結果。這個算法的最大問題在于處理大規模數據時遞歸深度過深的問題，根據這個算法，其遞歸樹的深度約為
?，如果要計算十億內質數的和，則遞歸深度要達到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以下的所有素數的個數的問題則是數論中的熱門話題了，相關文獻連篇累牘，因為這個問題在數論領域有相當的重要性，甚至還有一個專門的函數
?表示小于等于
?的所有素數的個數。高斯和勒讓德通過經驗統計的方式猜測
?，這個猜想在1896年得到證明成為素數定理，這是解析數論領域的最重要的成就之一。黎曼猜想也是在改進對
?的估計中被提出來的，現在應該是數論領域最重要的未被證明的猜想。除開解析數論的進路以外，很多數學家也在不斷改進計算?
的確切值的算法，相關的研究進展大家可以參見這個維基頁面，更深入的研究可以參見我在文末列出的參考文獻。
我們對計算N以下素數和算法來源于對勒讓德對計算N以下素數個數的算法的推廣。在勒讓德以前，數學家們計算?
的方法就是篩選出
?以下的所有素數然后數個數，勒讓德首次指出，為了計算?
，我們并不需要知道
?以下的所有素數，只需要知道?
就可以了。他給出的算法基于容斥原理。我們設
?表示小于等于?
的數中不能被?
整除的數的個數，其中
?表示前?
個質數，如
?等等。則有：
其中
?表示下取整。這個公式的意思是為了計算小于等于
?不能被?
整除的數的個數，我們從
?中減去可以
被?整除的數的個數，但是這樣同時被兩個質數整除的數就重復減去了，所以我們需要把它們的個數加回來，但是這樣又會導致對可以同時被三個質數整除重復加入了，所以我們需要減去這樣的數的個數，之后依次類推，顯然這只是容斥原理的一個簡單應用。如果真的要用這個公式來計算
?仍然顯得比較復雜，通過仔細分析
?的算法，我們可以發現一個遞推式。我們定義
?，而?
顯然表示小于等于?
的數中所有奇數的個數，則有：
這個公式同樣可以使用證明容斥原理時使用的數學歸納法加以證明，詳細的證明我就不寫了，只說明一下
?的簡單情況：
有了這個遞推式和上面給出的邊界條件，我們就可以計算出
?。如果我們設
?，則
?實際上表示是
?到?
之間素數的個數，則可以得到：
因而我們可據此算出
?以下的素數個數。可以看出，勒讓德的算法將?
以下所有質數分成了兩部分，然后分別計算它們的個數，加起來即為
?以下所有質數的個數。我們計算N以下質數和的算法也是基于同樣的原理，我們設?
表示小于等于
?的數中不能被
?整除的數的和，其中
?表示前
?個質數，則?
顯然表示小于等于?的數中所有奇數的和，則我們可以得到以下遞推式：
和上面類似，我們只說明一下
?的簡單情況，定義
?表示?
的自然數之和，如我們有
?，則有：
如我們設
?表示小于等于
?的所有素數的和，并設
?則根據和上面類似的原理，我們有：
據此我們可以求出小于等于
?的所有質數的和。可以看到這里的遞推式和Hedgehog算法的遞推式非常相似，區別在于這里的遞推式中
?表示素數，則不需要像Hedgehog算法那樣需要判斷是否為素數。更重要的區別在于如果使用遞歸實現Hedgehog算法，則其遞歸深度為?
，而在這里的算法中，遞歸深度為
?，前者明顯大于后者，因而這里的算法可以更快的達到邊界條件，并且因為素數越大則分布密度越低，則兩者在大規模數據中差距會更加明顯。如當
?時，
?，而
?，前者的遞歸深度是后者的四倍。
勒讓德算法的缺陷和第二個算法類似，都會在大規模數據上因為迭代深度的問題而無法計算，雖然勒讓德算法已經大大減少了遞歸函數的遞歸深度，但減少的仍然不夠。如要計算十億以內的素數和，則勒讓德算法的遞歸深度為
?，仍然超過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以上提到的勒讓德算法進行了改進，進一步提升了計算
?的效率，他使用自己改進的算法計算了
?，雖然比正確值小了56，不過考慮到他完全依靠手工計算，這個準確度已經非常驚人了。Meissel對勒讓德算法的主要改進是加入了一個新項
?，從而使得算法的時間復雜度從勒讓德算法的
?改進到
?，空間復雜度從
?改進至
?。這里我們對Meissel的算法做了推廣，使其可以計算N以下的質數和。首先我們對Meissel的算法做一個簡單介紹：
定義
?表示?
中恰好擁有?
個素因子，且這
?個素因子
?均大于
?的數的個數，則我們有：
這個公式我們可以這么理解：?
表示?
中其最小素因子都大于
?的數的個數，這些數里面包括只有一個大于?
的素因子的數，實際上也就是大于
?的素數；也包括恰好有兩個素素因子且這兩個素因子都大于
?，也包括恰好有三個素素因子且這三個素因子都大于
?，依次類推以至無窮就可以得到所有其最小素因子都大于
?的數的個數，也就是
?。我們定義?
，而
?表示?
中大于
?的素數的個數，則
?，據此我們展開上式有：
通過適當的選擇
?的數值，我們可以讓
?及以后的項變為零。如我們選擇
?，則有
?，此時
?及以后的變為零，上式變成之前提到的勒讓德的公式，證明勒讓德公式只是這個公式的一個特例。如我們選擇
?，則有?
，此時
?及以后的項變為零。一般地，選擇
?，則有?
而
?。假設我們選擇
?，則有：
經過不太復雜的推導，可以發現
?可通過下式計算：
其中
?，據此我們可以計算?
。使用和上面的類似的原理，并定義
?為區間?
中恰好有兩個素因子，且兩個素因子均大于?
的數的和，則有：
其中
?，且：
綜合上面兩個公式，我們也可以計算
?。這個算法相對于勒讓得算法的最顯著的優勢是需要遞歸的次數更少，如當
?時，
?，因此最大遞歸深度只有168層。但在我自己實現這個算法時，它的表現并沒有比勒讓得算法更優，我猜測原因是計算?
時消耗了過多資源，不過也有可能是我自己算法實現的原因。如果大家有提升這個算法的方法，還請指出。這個算法的代碼實現如下：
def meissel(n=2e6):
primes = list(primerange(1,int(n**(1/2))+1))
cr_primes = [x for x in primes if x
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算法的原理，并且對質數計數的各種算法和推廣加深了了解。我相信這些算法仍有優化的空間，如果找到的話我再來做補充。
五、延伸閱讀
前面已經提到，質數計數是數論中一個非常重要的問題，既有使用解析數論等方法對質數整體分布規律的研究，也有各種不斷改進的計算
?的確切值的算法。在以上我提到的Meissel算法之后約半個世紀，Lehmer[5]對這個算法進行了進一步改進，他進一步計算了
?，并使用IBM 701計算機計算了
?，他的計算結果只比正確結果大一，以上從勒讓德到Lehmer的算法實質都是對勒讓德最初提出算法的某種改進，統稱為計算
?的組合方法(Combinatorial method)。之后在1987年，Lagarias and Odlyzko[4]提出了一種從從解析數論角度計算
?的方法，算法中需要用到黎曼Zeta函數的某些性質，并采用數值積分的方法，這種方法被稱為計算?的分析方法(見參考文獻[1])，這個算法雖然具備更好的漸近性，但在性能上比不上作者們在1985年提出的對Meissel-Lehmer算法的改進(見[3])。上面幾種算法的詳細的時間與空間復雜度分析可以參見這篇文章，這里列一個簡要的結果：
在此之后，Deleglise-Rivat以及Xavier Gourdon又對算法做出了改進。在github上作者Kim Walisch對以上算法做了實現，他的測試結果如下：
在2015年9月，他宣布利用自己的算法計算了
?，計算花費了數年時間，最高峰時使用了235G的內存，之后他們又花了五個月時間進行重新計算驗證，確認的結果是：
這是目前維基頁面上列出的最大的
?值，至于是否有人算出了更大的?
值，我沒有查到確切的信息。
以上計算?的算法中自Lehmer以后都變得非常復雜，至少我已經無法理解。同時由于使用的方法越來越復雜，其代碼實現也會變得更加麻煩。而且相關方法能否進行推廣用來計算N以下的質數和也未可知，這些問題感興趣的同學可以繼續探索。
最后，參考文獻[7]，設
?，則
?以下質數的和的漸進估計是：
全文完。
六、參考文獻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.
總結
以上是生活随笔為你收集整理的java循环1000000000_求十亿内所有质数的和，怎么做最快?的全部內容，希望文章能夠幫你解決所遇到的問題。
- 上一篇： PyQt特殊对话框介绍
- 下一篇： c语言sizeof测量字符组长度,C语言