今天這篇是演算法與資料結構專題的第23篇文章，我們繼續數論相關的演算法，來看看大名鼎鼎的埃式篩法。

我們都知道在數學領域，素數非常重要，有海量的公式和研究關於素數，比如那個非常著名至今沒有人解出來的哥德巴赫猜想。和數學領域一樣，素數在資訊領域也非常重要，有著大量的應用。舉個簡單的例子，很多安全加密演算法也是利用的質數。我們想要利用素數去進行各種計算之前，總是要先找到素數。所以這就有了一個最簡單也最不簡單的問題，我們怎麼樣來尋找素數呢？

判斷素數

尋找素數最樸素的方法當然是一個一個遍歷，我們依次遍歷每一個數，然後分別判斷是否是素數。所以問題的核心又回到了判斷素數上，那麼怎麼判斷一個數是不是素數呢？

素數的性質只有一個，就是隻有1和它本身這兩個因數，我們要判斷素數也只能利用這個性質。所以可以想到，假如我們要判斷n是否是素數，可以從2開始遍歷到n-1，如果這n-1個數都不能整除n，那麼說明n就是素數。這個我沒記錯在C語言的練習題當中出現過，總之非常簡單，可以說是最簡單的演算法了。

def is_prime(n):
    for i in range(2, n):
        if n % i == 0:
            return False
    return n != 1

顯然，這個演算法是可以優化的，比如當n是偶數的時候，我們根本不需要迴圈，除了2意外的偶數一定是合數。再比如，我們迴圈的上界其實也沒有必要到n-1，到

這個改進也很簡單，稍作改動即可：

def is_prime(n):
    if n % 2 == 0 and n != 2:
        return False
    for i in range(3, int(math.sqrt(n) + 1)):
        if n % i == 0:
            return False
    return n != 1

這樣我們把O(n)的演算法優化到了O( )也算是有了很大的改進了，但是還沒有結束，我們還可以繼續優化。數學上有一個定理，只有形如6n-1和6n+1的自然數可能是素數，這裡的n是大於等於1的整數。

這個定理乍一看好像很高階，但其實很簡單，因為所有自然數都可以寫成6n,6n+1,6n+2,6n+3,6n+4,6n+5這6種，其中6n,6n+2,6n+4是偶數，一定不是素數。6n+3可以寫成3(2n+1)，顯然也不是素數，所以只有可能6n+1和6n+5可能是素數。6n+5等價於6n-1，所以我們一般寫成6n-1和6n+1。

利用這個定理，我們的程式碼可以進一步優化：

def is_prime(n):
    if n % 6 not in (1, 5) and n not in (2, 3):
        return False
    for i in range(3, int(math.sqrt(n) + 1)):
        if n % i == 0:
            return False
    return n != 1

雖然這樣已經很快了，但仍然不是最優的，尤其是當我們需要尋找大量素數的時候，仍會消耗大量的時間。那麼有沒有什麼辦法可以批量查詢素數呢？

有，這個方法叫做埃拉託斯特尼演算法。這個名字念起來非常拗口，這是一個古希臘的名字。此人是個古希臘的大牛，是大名鼎鼎的阿基米德的好友。他雖然沒有阿基米德那麼出名，但是也非常非常厲害，在數學、天文學、地理學、文學、歷史學等多個領域都有建樹。並且還自創方法測量了地球直徑、地月距離、地日距離以及黃赤交角等諸多數值。要知道他生活的年代是兩千五百多年前，那時候中國還是春秋戰國時期，可以想見此人有多厲害。

埃式篩法

我們今天要介紹的埃拉託斯特尼演算法就是他發明的用來篩選素數的方法，為了方便我們一般簡稱為埃式篩法或者篩法。埃式篩法的思路非常簡單，就是用已經篩選出來的素數去過濾所有能夠被它整除的數。這些素數就像是篩子一樣去過濾自然數，最後被篩剩下的數自然就是不能被前面素數整除的數，根據素數的定義，這些剩下的數也是素數。

舉個例子，比如我們要篩選出100以內的所有素數，我們知道2是最小的素數，我們先用2可以篩掉所有的偶數。然後往後遍歷到3，3是被2篩剩下的第一個數，也是素數，我們再用3去篩除所有能被3整除的數。篩完之後我們繼續往後遍歷，第一個遇到的數是7，所以7也是素數，我們再重複以上的過程，直到遍歷結束為止。結束的時候，我們就獲得了100以內的所有素數。

如果還不太明白，可以看下下面這張動圖，非常清楚地還原了這整個過程。

這個思想非常簡單，理解了之後寫出程式碼來真的很容易：

def eratosthenes(n):
    primes = []
    is_prime = [True] * (n + 1)
    for i in range(2, n+1):
        if is_prime[i]:
            primes.append(i)
            # 用當前素數i去篩掉所有能被它整除的數
            for j in range(i * 2, n+1, i):
                is_prime[j] = False
    return primes

我們執行一次程式碼看看：

和我們的預期一樣，獲得了小於100的所有素數。我們來分析一下篩法的複雜度，從程式碼當中我們可以看到，我們一共有了兩層迴圈，最外面一層迴圈固定是遍歷n次。而裡面的這一層迴圈遍歷的次數一直在變化，並且它的運算次數和素數的大小相關，看起來似乎不太方便計算。實際上是可以的，根據素數分佈定理以及一系列複雜的運算（相信我，你們不會感興趣的），我們是可以得出篩法的複雜度是 。

極致優化

篩法的複雜度已經非常近似 了，因為即使在n很大的時候，經過兩次ln的計算，也非常近似常數了，實際上在絕大多數使用場景當中，上面的演算法已經足夠應用了。

但是仍然有大牛不知滿足，繼續對演算法做出了優化，將其優化到了 的複雜度。雖然從效率上來看並沒有數量級的提升，但是應用到的思想非常巧妙，值得我們學習。在我們理解這個優化之前，先來看看之前的篩法還有什麼可以優化的地方。比較明顯地可以看出來，對於一個合數而言，它可能會被多個素數篩去。比如38，它有2和19這兩個素因數，那麼它就會被置為兩次False，這就帶來了額外的開銷，如果對於每一個合數我們只更新一次，那麼是不是就能優化到 了呢？

怎麼樣保證每個合數只被更新一次呢？這裡要用到一個定理，就是每個合數分解質因數只有的結果是唯一的。既然是唯一的，那麼一定可以找到最小的質因數，如果我們能夠保證一個合數只會被它最小的質因數更新為False，那麼整個優化就完成了。

那我們具體怎麼做呢？其實也不難，我們假設整數n的最小質因數是m，那麼我們用小於m的素數i乘上n可以得到一個合數。我們將這個合數消除，對於這個合數而言，i一定是它最小的質因數。因為它等於i * n，n最小的質因數是m，i 又小於m，所以i是它最小的質因數，我們用這樣的方法來生成消除的合數，這樣來保證每個合數只會被它最小的質因數消除。

根據這一點，我們可以寫出新的程式碼：

def ertosthenes(n):
    primes = []
    is_prime = [True] * (n+1)
    for i in range(2, n+1):
        if is_prime[i]:
            primes.append(i)
        for j, p in enumerate(primes):
            # 防止越界
            if p > n // i:
                break
            # 過濾
   is_prime[i * p] = False
            # 當i % p等於0的時候說明p就是i最小的質因數
            if i % p == 0:
                break
                
    return primes

總結

到這裡，我們關於埃式篩法的介紹就告一段落了。埃式篩法的優化版本相對來說要難以記憶一些，如果記不住的話，可以就只使用優化之前的版本，兩者的效率相差並不大，完全在可以接受的範圍之內。

篩法看著程式碼非常簡單，但是非常重要，有了它，我們就可以在短時間內獲得大量的素數，快速地獲得一個素數表。有了素數表之後，很多問題就簡單許多了，比如因數分解的問題，比如資訊加密的問題等等。我每次回顧篩法演算法的時候都會忍不住感慨，這個兩千多年前被髮明出來的演算法至今看來非但不過時，仍然還是那麼巧妙。希望大家都能懷著崇敬的心情，理解演算法當中的精髓。

如果喜歡本文，可以的話，請點個關注，給我一點鼓勵，也方便獲取更多文章。

本文使用 mdnice 排版