求十億內所有質數的和，怎麼做最快?

01-13

這個題目下的答案大致分為幾種：

Magic。例如 @陳碩， @渡子厄（半Magic，因為Wolfram Alpha並沒給出準確結果）。這個我就不評論了，因為沒說是什麼方法。
暴力1到10億，對每個數判斷是否素數（或暴力試除，或Miller Rabin）。方法太暴力，最快也不過 $Omega(nlog n)$ ，其中 $n$ 是上界10億。
篩出所有 $n$ 以內的素數然後加起來。有人用Eratosthenes篩（複雜度 $O(n log log n)$ ），有人用歐拉篩（最簡單的線性篩之一），也有人用Matlab等軟體。此方法也極其暴力。因為不超過n的素數有 $Theta(frac{n}{log n})$ 個，所以任何先找出所有素數再加起來的方法（即使使用亞線性的篩法，例如Atkin篩或Wheel篩）都不可能快於這個時間複雜度。
不知所云型。就不點名了。
理論上更優的演算法。目前只看到 @羅宸。他給出的代碼是來自於Project Euler第10題（Summation of primes ）的論壇中Lucy_Hedgehog給出的Python代碼的直接翻譯，注釋也都還在。接下來我就來介紹一下Lucy_Hedgehog給出的這個演算法。

==========正式回答分割線=============

首先來感受一下這個演算法的速度：

$ for n in 10000000 100000000 1000000000 10000000000 100000000000; do time python Main.py $n; done 3203324994356


real    0m0.049s

user    0m0.031s

sys     0m0.015s

279209790387276
real    0m0.117s

user    0m0.093s

sys     0m0.015s

24739512092254535
real    0m0.514s

user    0m0.468s

sys     0m0.046s

2220822432581729238
real    0m2.645s

user    0m2.625s

sys     0m0.031s

201467077743744681014

real 0m15.713s user 0m15.656s sys 0m0.031s

可以看到，即使是python這樣的解釋型動態語言，算出10億的結果也不過需要0.5s，算出1000億也只要15秒，是完虐以上各種暴力方法的。

整個演算法類似於Wikipedia中給出的求n以內素數個數的演算法（Prime-counting function），感興趣的也可以去看一下。

==========真正的！正式回答分割線 =============

定義 $S(v,p)$ 為 $2$ 到 $v$ 所有整數中，在普通篩法中外層循環篩完 $p$ 時仍然倖存的數的和。因此這些數要不本身是素數，要不其最小的素因子也大於 $p$ 。因此我們需要求的是 $S(n,lfloorsqrt n floor)$ ，其中 $n$ 是十億。

為了計算 $S(v,p)$ ，先考慮幾個特殊情況

$ple1$ 。此時所有數都還沒有被篩掉，所以 $S(v,p)=sum_{i=2}^{v}i=frac{(2+v)(v-1)}{2}$
$p$ 不是素數。因為篩法中 $p$ 早已被別的數篩掉，所以在這步什麼都不會做，所以此時 $S(v,p)=S(v,p-1)$ ；
$p$ 是素數，但是 $v<P>。因為每個合數都一定有一個不超過其平方根的素因子，如果篩到<img src=$ 時還沒篩掉一個數，那麼篩到 $p-1$ 時這個數也還在。所以此時也有 $S(v,p)=S(v,p-1)$ 。

現在考慮最後一種稍微麻煩些的情況： $p$ 是素數，且 $p^2le v$ 。

此時，我們要用素數 $p$ 去篩掉剩下的那些數中 $p$ 的倍數。注意到現在還剩下的合數都沒有小於 $p$ 的素因子。因此有：

$S(v,p)=S(v,p-1)-sum_{substack{2le k le v,\ pmbox{為}kmbox{的最小素因子}}}k$

後面那項中提取公共因子 $p$ ，有：

$S(v,p)=S(v,p-1)-p imessum_{substack{2le k le v,\ pmbox{為}kmbox{的最小素因子}}}frac{k}{p}$

因為 $p$ 整除 $k$ ，稍微變形一下，令 $t=frac{k}{p}$ ，有：

$S(v,p)=S(v,p-1)-p imessum_{substack{2le t le lfloorfrac{v}{p} floor,\ tmbox{的最小素因子}ge p}}t$

注意到一開始提到的 $S$ 的定義（「這些數要不本身是素數，要不其最小的素因子也大於(注意!) $p$ 」），此時 $p$ 後面這項可以用 $S$ 來表達：

$S(v,p)=S(v,p-1)-p imes(S(leftlfloorfrac{v}{p} ight floor,p-1)-{p-1mbox{以內的所有素數和}})$

再用 $S$ 替換素數和得到最終表達式：

$S(v,p)=S(v,p-1)-p imes(S(leftlfloorfrac{v}{p} ight floor,p-1)-S(p-1,p-1))$

我們最終的結果是 $S(n,lfloorsqrt n floor)$ 。計算時可以使用記憶化，也可以直接自底向上動態規劃。

至於演算法的複雜度就留作練習，是低於以上任何一種暴力方法的。

希望我講清楚了。代碼我就不放了，自己解決Project Euler 第10題之後可以去論壇里看Lucy_Hedgehog的代碼（在論壇第5頁）。也可以去看 @羅宸給出的代碼。

========== 更新：

@吳昌隆在評論中給出一個鏈接，其中包含一個理論上更優的演算法（Theoretical Computer Science Stack Exchange上Charles的回答）。我粗略讀了一下，認為理論上是可行的。因為本題中的答案只有 $O(frac{n^2}{log n})$ 級別，使用中國剩餘定理的話也只需要算到 $O(log n)$ 級別的素數，因此總複雜度為 $O(n^{1/2+epsilon}log n)$ 。

粗略讀了下其參考文獻（J. C. Lagarias and A. M. Odlyzko, Computing π(x): An analytic method, Journal of Algorithms8 (1987), pp. 173-191. ）後，我認為這演算法不可能像這個答案中的演算法一樣在十幾行代碼內實現，僅放出來供感興趣的朋友去看一下。

// 本代碼僅供娛樂 #include &


template&

struct ReverseBoolize {

    static const int i = 0;

};
template&<&>

struct ReverseBoolize&<0&> {

    static const int i = 1;

};
// N 是否整除 M

template&

struct Divisible {

    static const int i = ReverseBoolize&::i;

};
template&

struct CheckPrimeSum {

    static const int i = CheckPrimeSum&::i + Divisible&::i;

};
template&

struct CheckPrimeSum&<2, N&> {

    static const int i = Divisible&<2, N&>::i;

};
template&

struct IsPrime {

    static const int i = ReverseBoolize&::i&>::i;

};
template&<&>

struct IsPrime&<2&> {

    static const int i = 1;

};
template&

struct Select {

    static const int i = M;

};
template&

struct Select& {

    static const int i = 0;

};
template&

struct Sum {

    static const long i = Sum&::i + Select&::i&>::i;

};
template&<&>

struct Sum&<2&> {

    static const long i = 2;

};

int main(void) { std::cout &<&< Sum&<1000000000&>::i &<&< std::endl; return 0; }

sum(primes(10^9))

這就是我喜歡Matlab的原因。

Elapsed time is 18.441650 seconds.

題主也沒說是要程序運行的快還是要人寫得快對吧~

這道題最近討論過，存在比線性複雜度更優的演算法，js實現也能在2秒內跑完，占坑明天填...

----------------------------------------------------------------------------------------------------------------------------------------

先扔一個js版本（coffeescript寫的） Try Coffee （這個鏈接可能需要滑鼠在上面停幾秒鐘再點擊代碼才能載入出來），結果最後兩位不對，這是因為js用的double存整數，所以精度不夠，不過這沒關係，相同的演算法你用c/c++實現一下，用long long來算就好了。

----------------------------------------------------------------------------------------------------------------------------------------

發現菜魚ftfish 已經做了很詳細的講解了，那麼我就再補充一個更好理解一點的dp版本吧：Try Coffee （這個鏈接可能需要滑鼠在上面停幾秒鐘再點擊代碼才能載入出來），當然這個速度比較慢，僅供方便理解。

24739512092254535

這個問題問的如此直接粗暴，而碩爺的答案竟然更加直接粗暴。

我獃獃的停留在這個頁面大概 10 秒鐘，第一個想法竟然是按題主所言，動手寫個程序驗證下碩爺的答案是不是正確。

===========結果當然不言而喻啦（計時比較粗略哈）=============

說說思路：

1. 求和的前提是找出10億內的所有質數，那麼如何最快的找是關鍵。

2. 搜到尋找質數貌似最快的演算法：Sieve of Eratosthenes 。

3. 實現該演算法，封裝一下。

4. 把搜出來的全部質數都加起來。

代碼這樣：

#include &

#include & #include &


int main()

{

    primesieve::iterator pi;

    uint64_t sum = 0;

    uint64_t prime;
    clock_t begin = clock();
    while ((prime = pi.next_prime()) &< 1000000000ull)
        sum += prime;
    
    clock_t end = clock();
    double elapsed_secs = static_cast&(end - begin) / CLOCKS_PER_SEC;

std::cout &<&< "Sum of the primes below 10^9 = " &<&< sum &<&< " ("&<&< elapsed_secs &<&< " s)" &<&< std::endl; }

===========不要問我 primesieve.hpp 裡面的內容=============

我是站在巨人的肩膀上~ (羞走

用java實現的，不到1s

/* * sum: 24739512092254535 time:904 ms */ public static long primeSum(final int range) { int i, k; HashMap& S = new HashMap&(); int r = (int) Math.sqrt(range); int p = range / r; int[] V = new int[r + p - 1]; k = r + 1; for (i = 1; i &< k; i++) { V[i - 1] = range / i; } int count = 1; for (i = r + p - 2; i &>= r; i--) { V[i] = count++; } // ArrayUtils.printArray(V); for (i = 0; i &< V.length; i++) { S.put(V[i], ((long) V[i] * (V[i] + 1) / 2 - 1)); } // System.out.println(S); Long sp, p2; for (p = 2; p &< r + 1; p++) { if (S.get(p) &> S.get(p - 1)) { sp = S.get(p - 1); p2 = (long) (p * p); for (i = 0; i &< V.length; i++) { if (V[i] &< p2) { break; } S.put(V[i], S.get(V[i]) - p * (S.get(V[i] / p) - sp)); } } } return S.get(range); }

隨手寫的一個問題看著還挺火的，鄙視一下菜魚，逼本群換了個簡單的題目

import java.util.BitSet;


public class Main {
    public static void main(String[] args) {

        runTask("compute", () -&> System.out.println(compute(10_0000_0000)));

    }

private static long compute(int n) { if (n &< 2) { return 0; } BitSet bitSet = new BitSet(n + 1); final int sqrt = (int) Math.sqrt(n); long result = 2; int index = 3; while (index &<= n) { if (!bitSet.get(index)) { result += index; int add = index * 2; if (index &<= sqrt) { for (int temp = index * index; temp &<= n; temp += add) { bitSet.set(temp); } } } index += 2; } return resu< } public static void runTask(String taskName, Runnable runnable) { long start = System.currentTimeMillis(); runnable.run(); long end = System.currentTimeMillis(); System.out.printf("task: %s, cost %d ms. ", taskName, end - start); } }

輸出： 24739512092254535 task: compute, cost 6709 ms.

def P10(n): r = int(n**0.5) # assert r*r &<= n and (r+1)**2 &> n V = [n//i for i in range(1,r+1)] # print V V += list(range(V[-1]-1,0,-1)) # print V S = {i:i*(i+1)//2-1 for i in V} # print S for p in range(2,r+1): if S[p] &> S[p-1]: # p is prime sp = S[p-1] # sum of primes smaller than p p2 = p*p for v in V: if v &< p2: break S[v] -= p*(S[v//p] - sp) return S[n] N = 1000000000 print P10(N)

24739512092254535

[Finished in 0.8s]

題主不是想要代碼么？那我給題主發一份吧：

Public Function CalcSumOfPrime(Optional ByVal MaxTo As Long = 1000000000, Optional ByVal BufferSize As Long = -1) As String Dim Div1 As Long Dim Start As Double, Finish As Double Dim SumH As Long, SumL As Long


    Div1 = Sqr(MaxTo) + 1
    If BufferSize = -1 Then

        BufferSize = 1

        While BufferSize &< Div1
            BufferSize = BufferSize * 2
        Wend
    End If

    BufferSize = BufferSize * 8
    
    Start = Timer
    
    Dim I As Long, J As Long, K As Long, L As Long
    Dim Div2 As Long
    Dim DivArr() As Long
    Dim Count As Long
    Count = 1
    
    For I = 3 To Div1
        Div2 = Fix(Sqr(I)) + 1
        For J = 2 To Div2
            If I Mod J = 0 Then
                Exit For
            End If
        Next J
        If J &> Div2 Then

            Count = Count + 1

            ReDim Preserve DivArr(1 To Count) As Long

            DivArr(Count) = I

            SumL = SumL + I

        End If

    Next I

    DivArr(1) = 2

    SumL = SumL + 2

    Dim Map() As Boolean

    ReDim Map(0 To BufferSize - 1) As Boolean

    Dim LoopCount As Long

    For I = Div1 + 1 To MaxTo Step BufferSize

        If I + BufferSize &> MaxTo Then

            LoopCount = MaxTo - I

        Else

            LoopCount = BufferSize

        End If

        For J = 0 To BufferSize - 1

            Map(J) = False

        Next J

        For K = 1 To Count

            If I Mod DivArr(K) = 0 Then

                J = 0

            Else

                J = DivArr(K) - (I Mod DivArr(K))

            End If

            For J = J To LoopCount - 1 Step DivArr(K)

                Map(J) = True

            Next J

        Next K
        For J = 0 To LoopCount - 1

            If Map(J) = False Then

                SumL = SumL + I + J

                If SumL &> 1000000000 Then

                    SumH = SumH + 1

                    SumL = SumL - 1000000000

                End If

            End If

        Next J

    Next I

Finish = Timer MsgBox "Cost time = " (Finish - Start) " Sum = " SumH Format(SumL, "000000000") CalcSumOfPrime = SumH Format(SumL, "000000000") End Function

結果：

本機環境：i7 4770K，沒超頻。

解釋一下：

1、看不懂是什麼嗎？這是VB6的代碼，沒有環境也可以用VBA來跑（效率低，在我本機上是40-50s），比如Excel、Word什麼的都可以，OutLook也可以；

2、效率不算最快，但從內存使用率來說不算高，如果題主你肯把它改造成C語言版本的，應該能快一倍多；

3、思路嘛，就是篩法，分段篩而已，分多少段其實很有學問的；

4、我最擔心的就是我寫的這個東西，知乎上能看懂的好像不太多，VB6太小眾了。

題主你還不快過來感謝我。

上面已經有同學說了 Sieve of Eratosthenes ，於是我去看了下還有個叫做 Miller-Rabin primarily test 的演算法，試之，得結果

Sum is 24739512092254535, used PT394.850955S

T_T，為什麼他們的這麼快？

https://gist.github.com/zonyitoo/50736ea1709533844852

放代碼炸一下 Rust 粉

。。。。0.5s 你是用天河二號跑的吧！！

你說你是不是用的天河二號

是不是！

我應該是最次的吧。。。。

========= 經過@張永亮的提示修改後，能再節約20S=========

#include& #include& using namespace std; const unsigned long long int maxn = 1000000000; bool a[maxn]; int main(){ a[2]=true; for( int i = 3 ; i &< maxn ; i++ ) if(i%2!=0) a[i]=true; else a[i]=false; for( int i = 3 ; i * i &< maxn ; i++ ){ if( a[i] ) for( int j = i*i ; j &< maxn ; j+=i*2 ) a[j] = false; } unsigned long long int sum=0; for( int i = 2 ; i &<= maxn ; i++ ){ if( a[i] ){ sum+=i; } } printf( "%llu " , sum ); }

優化了一下少了25s 。。。。

i3 跑了44S i5跑了31S 你們可以試試。。。

據說這是一個python群的入群問題 @劉奕聰

用bash:

#!/bin/bash

print 24739512092254535

exit 0

逗逼們，這是python群的入群驗證問題，都被你們回答了，題目看來要改成100億了

線性篩法：

一邊生成質數表，一邊刪掉每個數的質數倍。複雜度顧名思義，O(N)

#include & #include & using namespace std;


const long long N = 1000000000;
bool sieve[N];
long long linear_sieve() {

    vector& prime;

for( int i = 2; i &< N; i++ ) { if( !sieve[i] ) prime.push_back( i ); for( int j = 0; i * prime[j] &< N; j++ ) { sieve[i * prime[j]] = true; if( i % prime[j] == 0 ) break; } } long long sum = 0; for( int i = 0; i &< ( int )prime.size(); ++i ) { sum += prime[i]; } return sum; } int main( void ) { printf( "%lld", linear_sieve() ); return 0; }

電腦CPU很水...i3-3217U

#define MAXVALUE 1000000000ull

int main() { bool * data = new bool[MAXVALUE]; for (int i = 0; i &< MAXVALUE; i++) { data[i] = true; } long long sum = 0; clock_t begin = clock(); int i, j; for (int i = 2; i &< MAXVALUE; i++) { if (data[i]) { sum += i; for (int j = i + i; j &< MAXVALUE; j += i) { data[j] = false; } } } clock_t end = clock(); double elapsed_secs = static_cast&(end - begin) / CLOCKS_PER_SEC; cout &<&< sum &<&< endl; cout &<&< elapsed_secs &<&< endl; delete[] data; system("pause"); }

答案和大神們的一樣。。問題是。。為什麼我這麼慢T T。。。

大圖預告。

唉各位大神請無視我。我賣個萌就去睡。

這個問題如果要寫程序的話，自然是會慢半拍啦。不如我們直接問wolfram alpha吧！

Dear Wolfram, could you please tell me what the sum of all primes less than 10^9 is?

matlab一句話確實很棒。但是開matlab也要30秒啊。

其餘的非腳本語言還要編譯。。。

當然啦，這裡有個問題：wolfram alpha需要付費才能用更多的機時來計算並查看詳細結果。

記得前幾年wolfram alpha剛推出的時候，每天可以詳細查看三個問題的答案。再多就要付費了。那時候剛上大學，年輕，看得我一愣一愣的。

感覺就好像在山洞裡撿到了一盞神燈，按回車後蹦出來一個燈神（指不定還是充氣的那種），說：

唉最近我的畫風是不是變得好贊？這時候應該推送一下微信號之類的。掃一掃關注以下微信號獲取Wolfram Alpha免費試用：

Wolfram|Alpha: Computational Knowledge Engine

這個答案肯定不是最快的，不過最近剛學haskell，剛好碰到了一樣的問題。

這裡可以感受一下lazy evaluation的魅力， input是所有的positive integers，但是一樣可以出結果

sieve :: [Integer] -&> [Integer] sieve [] = [] sieve (x:xs) = x : sieve (filter (y -&> y `mod` x /= 0) xs)

let primes = sieve [2..] sum (takeWhile (&<1000000000) primes)