[演算法] 線性篩

埃拉托斯特尼篩法

我們從大家國一就學過的埃拉托斯特尼篩法開始。
維護一個布林陣列not_prime，索引為$i$那欄表示數字$i$是否為質數，若為質數則該欄為0，若為合數則為1。

我之所以用0表示質數1表示合數是因為如此可省去初始化填入1的步驟。
而埃拉托斯特尼篩法的邏輯就是先預設所有數都是質數，再一步步把不是的劃掉

因為大家都學過這個篩法的步驟和原理(去問你國中數學老師)，以下直接上程式

std::vector<int> prime;
bool not_prime[N+1];

void Eratosthenes(int n) {
  not_prime[0] = not_prime[1] = true;
  for (int i=2; i<=n; i++) {
    if (!not_prime[i]) {
      prime.push_back(i); // 將質數加入 prime
      if (1LL*i*i > n) continue;
      for (int j=i*i; j<=n; j+=i) not_prime[j] = true;
      // 之所以從i*i開始是因為比這個數小的i的倍數i*k (k<i) 都已經被篩過了
    }
  }
}

可以用 Mertens 第二定理證明其時間複雜度為$O(n\log\log n)$

我不會證

線性篩(歐拉篩法)

因為我這裡卡了很久，所以希望能把它解釋清楚讓更多人能茅塞頓開。

線性篩是埃拉托斯特尼篩法的優化版本。
埃拉托斯特尼篩法中，一個合數會被多次篩到，以$60$為例，當迴圈跑到i=2時，會跑一次not_prime[60] = true;，i=3時也會再次執行not_prime[60] = true;，i=5時也會，但除了第一次，之後都只是花時間在做一件已經完成的事情(也就是重複的跑not_prime[60] = true;)。
線性篩的想法很簡單：

要是我們能讓每個合數都只被他最小的質因數篩掉就好了

如此，篩法的時間複雜度就會降到$O(n)$

原理解說

每個合數$K$都可以被分解成$p_K\cdot i_K$，其中$p_K$表示$K$的最小質因數。
由於$p_K$是$K$的最小質因數，為$K$最小的非1的因數，從而$i_K$一定大於$p_K$

我們以$i_K$對合數$K$進行分類，被分到同一類的合數，都能藉由同一個$i_K$，乘以一個「比$i_K$小的質數$p_K$」表示出來。
因此我們能用一個for迴圈，當$i=i_K$時，將所有可用$i_K$乘一個質數表示出來的合數$K=p\cdot i_K$找出來。
$p_K$比$i_K$小是這個演算法能成功跑下去的關鍵，能讓i隨著for迴圈慢慢變大的過程中，動態添加新質數，而不用擔心會漏篩。
什麼意思呢？我們手動試一下就知道了。

• 首先$i=2$，將2加入prime。此時prime = {2}。將所有可以被表示為$p\cdot 2$的合數標記出來，其中$p$為小於等於$i(=2)$的質數。因此只有4被標記為合數：not_prime[2*2] = 1。
• 接著$i=3$，將3加入prime，此時prime = {2,3}。將所有可以被表示為$p\cdot 3$的合數標記出來，其中$p\leq i(=3)$。因此6、9被標記為合數。
• 接著$i=4$，4已被標記為合數，因此不加入prime，此時prime = {2,3}。將所有可以被表示為$p\cdot 4$的合數標記出來，其中$p$為小於等於$i(=4)$的質數。因此8、12被標記為合數。(注意：此時我們先忽略後面要加入的break。實際線性篩中 i=4 只會對8標記，而不會對12標記)
• 接著$i=5$，5還未被標記為合數，故為質數，加入prime，此時prime = {2,3,5}。將所有可以被表示為$p\cdot 5$的合數標記出來，其中$p$為小於等於$i(=5)$的質數。因此10、15、25被標記為合數。

不斷重複，直到$i=N$。

這裡我們解釋兩個細節：

1. 迴圈跑到$i=i_K$時，如果$i_K$是合數，他必可拆成比$i_K$小的質數和另一數的乘積，從而在之前的步驟就被標記出來了，所以若$i_K$沒被標記為合數，他就一定是質數。稍微延伸一下可知此時所有小於$i_K$的質數此時也都被找出來了。
2. 由於$p_K$比$i_K$小，所以當$i$跑到$i_K$時，所有比$i_K$小的質數$p_K$都已經在prime裡了，所以不會漏篩合數。

從上我們知道了線性篩的步驟，也能理解線性篩是能正確運作的。現在剩最後一個問題：
如何確保每個數「只被最小質因數篩到」？
我們回到先前的例子，從$i=6$繼續跑下去：

• 接著$i=6$，6已被標記為合數，因此不加入prime，此時prime = {2,3,5}。將所有可以被表示為$p\cdot 6$的合數標記出來，其中$p$為小於等於$i(=6)$的質數。因此12, 18, 30被標記為合數。

我們發現12被重複標記了，往回看可以看到當$i=4$時，$12(=3\cdot 4)$被3這個質數篩掉了，但我們希望12是被他的最小質因數2篩掉啊！這代表我們在i=4時要做一些事，讓12這個數被留到i=6，才由2這個質數(因為$12=2\cdot 6$，2是12的最小質因數)篩掉。
我們可以注意到prime中的質數是由小到大排列的，所以我們可以在迴圈中多加一個判斷式：對於每個質數$p$，都判斷$p$是否整除$i$，這樣找出來的第一個質數便為$i$的最小質因數。一旦遇到了$i$的最小質因數，後面的數就不用再乘了，因為這些乘出來的數要留給他們的最小質因數篩。
看上去很抽象，我們一樣以小數字舉例：

• i=4時，首先乘以p=2得到8，篩掉，沒問題，接著判斷發現$2\mid4$，所以直接break;(即不用考慮後面的$3\cdot 4=12$，因為$i(=4)$的最小質因數已經為$2$，所以$p\cdot i$的最小質因數不會比$2$還大，故不該被比2大的質數(如這裡的3)給篩掉)
• i=6時，首先乘以p=2得到12，篩掉，沒問題，接著判斷發現$2\mid6$，所以直接break;(後面的$3\cdot 6$跟$5\cdot 6$的最小質因數都2，因此不該在這時由3、5篩掉)

Code

#include<bits/stdc++.h>

std::vector<int> prime;
bool not_prime[N+1];

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

線性篩求歐拉函數

在線性篩的過程中結合動態規劃，我們能達到意想不到的效果。首先是「歐拉函數」的計算。
歐拉函數$\phi(n)$表示小於等於$n$且與$n$互質的正整數個數。比如$\phi(10)=4$
因為1~10中，1,3,7,9共四個數與10互質。

歐拉函數的性質

1. 若$p$為質數，$\phi(p)=p-1$
2. 若$p$為質數，$\phi(p^k)=p^k-p^{k-1}=p^k\left(1-\frac{1}{p}\right)$
3. 歐拉函數為積性函數：若$(a,b)=1$，則$\phi(ab)=\phi(a)\phi(b)$
4. 歐拉函數一般式：$\phi(n)=n\prod\limits_{p\in S}\left({1-\frac{1}{p}}\right)$，$S$為$n$的所有質因數的集合

原理分析

我們的目標是結合動態規劃(用前一個算後面的)，並運用到線性篩能找出最小質因數的特性。
設$p_n$為$n$的最小質因數，令$n^\prime=\frac{n}{p_n}$。
分為以下兩種情形：
情形一：$p_n \mid n^\prime$
此時$n$的所有質因數皆亦為$n^\prime$的，因此

$$\phi(n)=n\prod\left({1-\frac{1}{p_i}}\right)=p_nn^\prime\prod\left({1-\frac{1}{p_i}}\right)=p_n\phi(n^\prime)$$

情形二：$p_n \nmid n^\prime$
此時$n^\prime$和$n$互質，因此

$$\phi(n)=\phi(p_n)\phi(n^\prime)=(p_n-1)\phi(n^\prime)$$

我們知道線性篩能夠將每個數都標記一次(看他是質數 / 合數)，因此我們多加一個步驟：在標記同時記錄那個數的歐拉函數值是多少。
和之前相同，我們用i乘以質數表中的數字p來生成出所有的合數。
最後一個問題是，我們要確定在計算$\phi(n)$時，$\phi(n^\prime)$已經算好了，解釋也很簡單
由於$n=p_nn^\prime$，$p_n$為$n$的最小質因數，所以「計算$\phi(n)$」發生在迴圈中i跑到$n^\prime$時。
我們再來看$\phi(n^\prime)$是在何時被算出來的。
若$n^\prime$為合數，則$n^\prime$的最小質因數至少為2，故$\phi(n^\prime)$必定在迴圈執行到$i=\lceil\frac{n^\prime}{2}\rceil\leq n^\prime$之前就算出來了。
若$n^\prime$為質數，則我們只要讓i跑到$n^\prime$時，先計算$\phi(n^\prime)=n^\prime-1$，再計算$\phi(n)$，就能保證$\phi(n^\prime)$在$\phi(n)$之前算出來。

程式實作

#include<bits/stdc++.h>

std::vector<int> prime;
bool not_prime[N+1];
int phi[N+1];

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

線性篩求正因數個數

說明

用類似的分析方法，這次我們維護三個陣列：not_prime、num、d，其中num紀錄num[k]表示$k$的最小質因數在$k$的標準分解式的次數，d[k]表示$k$的正因數個數。

1. 對於質數$p$，$num(p) = 1, d(p) = 2$
2. 對於合數$n=p_nn^\prime$，其中$p_n$為$n$的最小質因數
   1. 若$p_n \mid n^\prime$，則$num(n)=num(n^\prime)+1$、$d(n)=d(n^\prime)\frac{num(n)+1}{num(n^\prime)+1}$
   2. 若$p_n \nmid n^\prime$，則$num(n)=1$、$d(n)=d(n^\prime)(num(n)+1)=2d(n^\prime)$

程式

#include<bits/stdc++.h>

std::vector<int> prime;
bool not_prime[N+1];
int num[N+1];
int d[N+1];

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

線性篩求積性函數

從上面兩個例子，我們其實可以將其推廣到更多的積性函數$f(x)$(如正因數和、莫比烏斯函數等)，分析步驟如下：

1. 考慮怎麼算$f(p)$，其中$p$為質數
2. 考慮將合數拆解成$n=p_nn^\prime$，其中$p_n$為$n$的最小質因數，怎麼在已知$f(n^\prime)$下計算出$f(n)$
   • 分別考慮$p_n \mid n^\prime$和$p_n \nmid n^\prime$的情形

接著照樣造程式即可！

參考資料

• 筛法 - OI Wiki
• Linear Sieve - NTHU CPP
• 線性篩法 (不專業講解) - terryobeyes的創作 - 巴哈姆特