跳转至

篩法

素數篩法

引入

如果我們想要知道小於等於 \(n\) 有多少個素數呢？

一個自然的想法是對於小於等於 \(n\) 的每個數進行一次質數檢驗。這種暴力的做法顯然不能達到最優複雜度。

埃拉託斯特尼篩法

過程

考慮這樣一件事情：對於任意一個大於 \(1\) 的正整數 \(n\)，那麼它的 \(x\) 倍就是合數（\(x > 1\)）。利用這個結論，我們可以避免很多次不必要的檢測。

如果我們從小到大考慮每個數，然後同時把當前這個數的所有（比自己大的）倍數記為合數，那麼運行結束的時候沒有被標記的數就是素數了。

實現

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
vector<int> prime;
bool is_prime[N];

void Eratosthenes(int n) {
  is_prime[0] = is_prime[1] = false;
  for (int i = 2; i <= n; ++i) is_prime[i] = true;
  for (int i = 2; i <= n; ++i) {
    if (is_prime[i]) {
      prime.push_back(i);
      if ((long long)i * i > n) continue;
      for (int j = i * i; j <= n; j += i)
        // 因為從 2 到 i - 1 的倍數我們之前篩過了，這裏直接從 i
        // 的倍數開始，提高了運行速度
        is_prime[j] = false;  // 是 i 的倍數的均不是素數
    }
  }
}

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
prime = []
is_prime = [False] * N

def Eratosthenes(n):
    is_prime[0] = is_prime[1] = False
    for i in range(2, n + 1):
        is_prime[i] = True
    for i in range(2, n + 1):
        if is_prime[i]:
            prime.append(i)
            if i * i > n:
                continue
            for j in range(i * i, n + 1, i):
                is_prime[j] = False

以上為 Eratosthenes 篩法（埃拉託斯特尼篩法，簡稱埃氏篩法），時間複雜度是 \(O(n\log\log n)\)。

證明

現在我們就來看看推導過程：

如果每一次對數組的操作花費 1 個單位時間，則時間複雜度為：

\[ O\left(\sum_{k=1}^{\pi(n)}{\frac{n}{p_k}}\right)=O\left(n\sum_{k=1}^{\pi(n)}{\frac{1}{p_k}}\right) \]

其中 \(p_k\) 表示第 \(k\) 小的素數，\(\pi(n)\) 表示 \(\le n\) 的素數個數。\(\sum_{k=1}^{\pi(n)}\) 表示第一層 for 循環，其中累加上界 \(\pi(n)\) 為 if (prime[i]) 進入 true 分支的次數；\(\frac{n}{p_k}\) 表示第二層 for 循環的執行次數。

根據 Mertens 第二定理，存在常數 \(B_1\) 使得：

\[ \sum_{k=1}^{\pi(n)}{\frac{1}{p_k}}=\log\log n+B_1+O\left(\frac{1}{\log n}\right) \]

所以 Eratosthenes 篩法 的時間複雜度為 \(O(n\log\log n)\)。接下來我們證明 Mertens 第二定理的弱化版本 \(\sum_{k\le\pi(n)}1/p_k=O(\log\log n)\)：

根據 \(\pi(n)=\Theta(n/\log n)\)，可知第 \(n\) 個素數的大小為 \(\Theta(n\log n)\)。於是就有

\[ \begin{aligned} \sum_{k=1}^{\pi(n)}{\frac{1}{p_k}} &=O\left(\sum_{k=2}^{\pi(n)}{\frac{1}{k\log k}}\right) \\ &=O\left(\int_2^{\pi(n)}{\frac{\mathrm dx}{x\log x}}\right) \\ &=O(\log\log\pi(n))=O(\log\log n) \end{aligned} \]

當然，上面的做法效率仍然不夠高效，應用下面幾種方法可以稍微提高算法的執行效率。

篩至平方根

顯然，要找到直到 \(n\) 為止的所有素數，僅對不超過 \(\sqrt n\) 的素數進行篩選就足夠了。

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
vector<int> prime;
bool is_prime[N];

void Eratosthenes(int n) {
  is_prime[0] = is_prime[1] = false;
  for (int i = 2; i <= n; ++i) is_prime[i] = true;
  // i * i <= n 説明 i <= sqrt(n)
  for (int i = 2; i * i <= n; ++i) {
    if (is_prime[i])
      for (int j = i * i; j <= n; j += i) is_prime[j] = false;
  }
  for (int i = 2; i <= n; ++i)
    if (is_prime[i]) prime.push_back(i);
}

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
prime = []
is_prime = [False] * N

def Eratosthenes(n):
    is_prime[0] = is_prime[1] = False
    for i in range(2, n + 1):
        is_prime[i] = True
    # 讓 i 循環到 <= sqrt(n)
    for i in range(2, isqrt(n) + 1): # `isqrt` 是 Python 3.8 新增的函數
        if is_prime[i]:
            for j in range(i * i, n + 1, i):
                is_prime[j] = False
    for i in range(2, n + 1):
        if is_prime[i]:
            prime.append(i)

這種優化不會影響漸進時間複雜度，實際上重複以上證明，我們將得到 \(n \ln \ln \sqrt n + o(n)\)，根據對數的性質，它們的漸進相同，但操作次數會明顯減少。

只篩奇數

因為除 2 以外的偶數都是合數，所以我們可以直接跳過它們，只用關心奇數就好。

首先，這樣做能讓我們內存需求減半；其次，所需的操作大約也減半。

減少內存的佔用

我們注意到篩選時只需要 bool 類型的數組。bool 數組的一個元素一般佔用 \(1\) 字節（即 \(8\) 比特），但是存儲一個布爾值只需要 \(1\) 個比特就足夠了。

我們可以使用 位運算 的相關知識，將每個布爾值壓到一個比特位中，這樣我們僅需使用 \(n\) 比特（即 \(\dfrac n 8\) 字節）而非 \(n\) 字節，可以顯著減少內存佔用。

但是，這種稱為 位級壓縮 的方法會使這些位的操作複雜化。任何位上的讀寫操作都需要多次算術運算，最終會使算法變慢。因此，這種方法只有在 \(n\) 特別大，以至於我們不能再分配內存時才合理。在這種情況下，我們將犧牲效率，通過顯著降低算法速度以節省內存（減小到原來的 \(\dfrac n 8\)）。

值得一提的是，存在自動執行位級壓縮的數據結構，如 C++ 中的 vector<bool> 和 bitset<>（參見 bitset: 與埃氏篩結合）。

分塊篩選

由優化「篩至平方根」可知，不需要一直保留整個 is_prime[1...n] 數組。為了進行篩選，只保留到 \(\sqrt n\) 的素數就足夠了，即 prime[1...sqrt(n)]。並將整個範圍分成塊，每個塊分別進行篩選。這樣，我們就不必同時在內存中保留多個塊，而且 CPU 可以更好地處理緩存。

設 \(s\) 是一個常數，它決定了塊的大小，那麼我們就有了 \(\lceil {\frac n s} \rceil\) 個塊，而塊 \(k\)(\(k = 0 \dots \lfloor {\frac n s} \rfloor\)) 包含了區間 \([ks, ks + s - 1]\) 中的數字。我們可以依次處理塊，也就是説，對於每個塊 \(k\)，我們將遍歷所有質數（從 \(1\) 到 \(\sqrt n\)）並使用它們進行篩選。

值得注意的是，我們在處理第一個數字時需要稍微修改一下策略：首先，應保留 \([1, \sqrt n]\) 中的所有的質數；第二，數字 \(0\) 和 \(1\) 應該標記為非素數。在處理最後一個塊時，不應該忘記最後一個數字 \(n\) 並不一定位於塊的末尾。

以下實現使用塊篩選來計算小於等於 \(n\) 的質數數量。

實現 
 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
int count_primes(int n) {
  const int S = 10000;
  vector<int> primes;
  int nsqrt = sqrt(n);
  vector<char> is_prime(nsqrt + 1, true);
  for (int i = 2; i <= nsqrt; i++) {
    if (is_prime[i]) {
      primes.push_back(i);
      for (int j = i * i; j <= nsqrt; j += i) is_prime[j] = false;
    }
  }
  int result = 0;
  vector<char> block(S);
  for (int k = 0; k * S <= n; k++) {
    fill(block.begin(), block.end(), true);
    int start = k * S;
    for (int p : primes) {
      int start_idx = (start + p - 1) / p;
      int j = max(start_idx, p) * p - start;
      for (; j < S; j += p) block[j] = false;
    }
    if (k == 0) block[0] = block[1] = false;
    for (int i = 0; i < S && start + i <= n; i++) {
      if (block[i]) result++;
    }
  }
  return result;
}

分塊篩法的漸進時間複雜度與埃氏篩法是一樣的（除非塊非常小），但是所需的內存將縮小為 \(O(\sqrt{n} + S)\)，並且有更好的緩存結果。 另一方面，對於每一對塊和區間 \([1, \sqrt{n}]\) 中的素數都要進行除法，而對於較小的塊來説，這種情況要糟糕得多。 因此，在選擇常數 \(S\) 時要保持平衡。

塊大小 \(S\) 取 \(10^4\) 到 \(10^5\) 之間，可以獲得最佳的速度。

線性篩法

埃氏篩法仍有優化空間，它會將一個合數重複多次標記。有沒有什麼辦法省掉無意義的步驟呢？答案是肯定的。

如果能讓每個合數都只被標記一次，那麼時間複雜度就可以降到 \(O(n)\) 了。

實現 
 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
vector<int> pri;
bool not_prime[N];

void pre(int n) {
  for (int i = 2; i <= n; ++i) {
    if (!not_prime[i]) {
      pri.push_back(i);
    }
    for (int pri_j : pri) {
      if (i * pri_j > n) break;
      not_prime[i * pri_j] = true;
      if (i % pri_j == 0) {
        // i % pri_j == 0
        // 換言之，i 之前被 pri_j 篩過了
        // 由於 pri 裏面質數是從小到大的，所以 i 乘上其他的質數的結果一定會被
        // pri_j 的倍數篩掉，就不需要在這裏先篩一次，所以這裏直接 break
        // 掉就好了
        break;
      }
    }
  }
}

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
pri = []
not_prime = [False] * N

def pre(n):
    for i in range(2, n + 1):
        if not not_prime[i]:
            pri.append(i)
        for pri_j in pri:
            if i * pri_j > n:
                break
            not_prime[i * pri_j] = True
            if i % pri_j == 0:
                """
                i % pri_j == 0
                換言之，i 之前被 pri_j 篩過了
                由於 pri 裏面質數是從小到大的，所以 i 乘上其他的質數的結果一定會被
                pri_j 的倍數篩掉，就不需要在這裏先篩一次，所以這裏直接 break
                掉就好了
                """
                break

上面的這種 線性篩法 也稱為 Euler 篩法（歐拉篩法）。

Note

注意到篩法求素數的同時也得到了每個數的最小質因子。

篩法求歐拉函數

注意到在線性篩中，每一個合數都是被最小的質因子篩掉。比如設 \(p_1\) 是 \(n\) 的最小質因子，\(n' = \frac{n}{p_1}\)，那麼線性篩的過程中 \(n\) 通過 \(n' \times p_1\) 篩掉。

觀察線性篩的過程，我們還需要處理兩個部分，下面對 \(n' \bmod p_1\) 分情況討論。

如果 \(n' \bmod p_1 = 0\)，那麼 \(n'\) 包含了 \(n\) 的所有質因子。

\[ \begin{aligned} \varphi(n) & = n \times \prod_{i = 1}^s{\frac{p_i - 1}{p_i}} \\\\ & = p_1 \times n' \times \prod_{i = 1}^s{\frac{p_i - 1}{p_i}} \\\\ & = p_1 \times \varphi(n') \end{aligned} \]

那如果 \(n' \bmod p_1 \neq 0\) 呢，這時 \(n'\) 和 \(p_1\) 是互質的，根據歐拉函數性質，我們有：

\[ \begin{aligned} \varphi(n) & = \varphi(p_1) \times \varphi(n') \\\\ & = (p_1 - 1) \times \varphi(n') \end{aligned} \]

實現

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
vector<int> pri;
bool not_prime[N];
int phi[N];

void pre(int n) {
  phi[1] = 1;
  for (int i = 2; i <= n; i++) {
    if (!not_prime[i]) {
      pri.push_back(i);
      phi[i] = i - 1;
    }
    for (int pri_j : pri) {
      if (i * pri_j > n) break;
      not_prime[i * pri_j] = true;
      if (i % pri_j == 0) {
        phi[i * pri_j] = phi[i] * pri_j;
        break;
      }
      phi[i * pri_j] = phi[i] * phi[pri_j];
    }
  }
}

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
pri = []
not_prime = [False] * N
phi = [0] * N

def pre(n):
    phi[1] = 1
    for i in range(2, n + 1):
        if not not_prime[i]:
            pri.append(i)
            phi[i] = i - 1
        for pri_j in pri:
            if i * pri_j > n:
                break
            not_prime[i * pri_j] = True
            if i % pri_j == 0:
                phi[i * pri_j] = phi[i] * pri_j
                break
            phi[i * pri_j] = phi[i] * phi[pri_j]

篩法求莫比烏斯函數

定義

根據莫比烏斯函數的定義，設 \(n\) 是一個合數，\(p_1\) 是 \(n\) 的最小質因子，\(n'=\frac{n}{p_1}\)，有：

\[ \mu(n)= \begin{cases} 0 & n' \bmod p_1 = 0\\\\ -\mu(n') & \text{otherwise} \end{cases} \]

若 \(n\) 是質數，有 \(\mu(n)=-1\)。

實現

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
vector<int> pri;
bool not_prime[N];
int mu[N];

void pre(int n) {
  mu[1] = 1;
  for (int i = 2; i <= n; ++i) {
    if (!not_prime[i]) {
      mu[i] = -1;
      pri.push_back(i);
    }
    for (int pri_j : pri) {
      if (i * pri_j > n) break;
      not_prime[i * pri_j] = true;
      if (i % pri_j == 0) {
        mu[i * pri_j] = 0;
        break;
      }
      mu[i * pri_j] = -mu[i];
    }
  }
}

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
pri = []
not_prime = [False] * N
mu = [0] * N

def pre(n):
    mu[1] = 1
    for i in range(2, n + 1):
        if not not_prime[i]:
            pri.append(i)
            mu[i] = -1
        for pri_j in pri:
            if i * pri_j > n:
                break
            not_prime[i * pri_j] = True
            if i % pri_j == 0:
                mu[i * pri_j] = 0
                break
            mu[i * pri_j] = -mu[i]

篩法求約數個數

用 \(d_i\) 表示 \(i\) 的約數個數，\(num_i\) 表示 \(i\) 的最小質因子出現次數。

約數個數定理

定理：若 \(n=\prod_{i=1}^m p_i^{c_i}\) 則 \(d_i=\prod_{i=1}^m (c_i+1)\)。

證明：我們知道 \(p_i^{c_i}\) 的約數有 \(p_i^0,p_i^1,\dots ,p_i^{c_i}\) 共 \(c_i+1\) 個，根據乘法原理，\(n\) 的約數個數就是 \(\prod_{i=1}^m (c_i+1)\)。

實現

因為 \(d_i\) 是積性函數，所以可以使用線性篩。

在這裏簡單介紹一下線性篩實現原理。

  1. 當 \(i\) 為質數時，\(\textit{num}_i \gets 1,\textit{d}_i \gets 2\)，同時設 \(q = \left\lfloor \dfrac {i}{p} \right\rfloor\)，其中 \(p\) 為 \(i\) 的最小質因子。
  2. 當 \(p\) 為 \(q\) 的質因子時，\(\textit{num}_i \gets \textit{num}_q + 1,\textit{d}_i \gets \dfrac{\textit{d}_q}{\textit{num}_i} \times (\textit{num}_i + 1)\)。
  3. 當 \(p,q\) 互質時，\(\textit{num}_i \gets 1,\textit{d}_i \gets \textit{d}_q \times (\textit{num}_i+1)\)。
 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
vector<int> pri;
bool not_prime[N];
int d[N], num[N];

void pre(int n) {
  d[1] = 1;
  for (int i = 2; i <= n; ++i) {
    if (!not_prime[i]) {
      pri.push_back(i);
      d[i] = 2;
      num[i] = 1;
    }
    for (int pri_j : pri) {
      if (i * pri_j > n) break;
      not_prime[i * pri_j] = true;
      if (i % pri_j == 0) {
        num[i * pri_j] = num[i] + 1;
        d[i * pri_j] = d[i] / num[i * pri_j] * (num[i * pri_j] + 1);
        break;
      }
      num[i * pri_j] = 1;
      d[i * pri_j] = d[i] * 2;
    }
  }
}

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
pri = []
not_prime = [False] * N
d = [0] * N
num = [0] * N

def pre(n):
    d[1] = 1
    for i in range(2, n + 1):
        if not not_prime[i]:
            pri.append(i)
            d[i] = 2
            num[i] = 1
        for pri_j in pri:
            if i * pri_j > n:
                break
            not_prime[i * pri_j] = True
            if i % pri_j == 0:
                num[i * pri_j] = num[i] + 1
                d[i * pri_j] = d[i] // num[i * pri_j] * (num[i * pri_j] + 1)
                break
            num[i * pri_j] = 1
            d[i * pri_j] = d[i] * 2

篩法求約數和

\(f_i\) 表示 \(i\) 的約數和，\(g_i\) 表示 \(i\) 的最小質因子的 \(p^0+p^1+p^2+\dots p^k\).

實現

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
vector<int> pri;
bool not_prime[N];
int g[N], f[N];

void pre(int n) {
  g[1] = f[1] = 1;
  for (int i = 2; i <= n; ++i) {
    if (!not_prime[i]) {
      pri.push_back(i);
      g[i] = i + 1;
      f[i] = i + 1;
    }
    for (int pri_j : pri) {
      if (i * pri_j > n) break;
      not_prime[i * pri_j] = true;
      if (i % pri_j == 0) {
        g[i * pri_j] = g[i] * pri_j + 1;
        f[i * pri_j] = f[i] / g[i] * g[i * pri_j];
        break;
      }
      f[i * pri_j] = f[i] * f[pri_j];
      g[i * pri_j] = 1 + pri_j;
    }
  }
}

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
pri = []
not_prime = [False] * N
f = [0] * N
g = [0] * N

def pre(n):
    g[1] = f[1] = 1
    for i in range(2, n + 1):
        if not not_prime[i]:
            pri.append(i)
            g[i] = i + 1
            f[i] = i + 1
        for pri_j in pri:
            if i * pri_j > n:
                break
            not_prime[i * pri_j] = True
            if i % pri_j == 0:
                g[i * pri_j] = g[i] * pri_j + 1
                f[i * pri_j] = f[i] // g[i] * g[i * pri_j]
                break
            f[i * pri_j] = f[i] * f[pri_j]
            g[i * pri_j] = 1 + pri_j

一般的積性函數

假如一個 積性函數 \(f\) 滿足：對於任意質數 \(p\) 和正整數 \(k\)，可以在 \(O(1)\) 時間內計算 \(f(p^k)\)，那麼可以在 \(O(n)\) 時間內篩出 \(f(1),f(2),\dots,f(n)\) 的值。

設合數 \(n\) 的質因子分解是 \(\prod_{i=1}^k p_i^{\alpha_i}\)，其中 \(p_1<p_2<\dots<p_k\) 為質數，我們在線性篩中記錄 \(g_n=p_1^{\alpha_1}\)，假如 \(n\) 被 \(x\cdot p\) 篩掉（\(p\) 是質數），那麼 \(g\) 滿足如下遞推式：

\[ g_n= \begin{cases} g_x\cdot p & x\bmod p=0\\\\ p & \text{otherwise} \end{cases} \]

假如 \(n=g_n\)，説明 \(n\) 就是某個質數的次冪，可以 \(O(1)\) 計算 \(f(n)\)；否則，\(f(n)=f(\frac{n}{g_n})\cdot f(g_n)\)。

本節部分內容譯自博文 Решето Эратосфена 與其英文翻譯版 Sieve of Eratosthenes。其中俄文版版權協議為 Public Domain + Leave a Link；英文版版權協議為 CC-BY-SA 4.0。

本页面最近更新：，更新历史
发现错误？想一起完善？ 在 GitHub 上编辑此页！
本页面贡献者：inkydragon, TravorLZH, YOYO-UIAT, wood3, shuzhouliu, Mr-Python-in-China, HeRaNO
本页面的全部内容在 CC BY-SA 4.0 和 SATA 协议之条款下提供，附加条款亦可能应用