title: 質數篩法
date: 2021-09-16
tags:

• 數論
• 學習筆記
• 演算法
  categories:
  學習筆記
  mathjax: true

簡介：質數篩

質數 。
一類很奇異的數字。
質數的定義是：
“在大於1的自然數中，除了1和它本身以外不再有其他因數的自然數。”

於是我們就可以想到一個簡單的演算法來求 \(n\) 是不是質數：

簡單的演算法

for(int i = 2; i <= sqrt(n); i++)
{
   if(n % i == 0)
   {
       cout << "いいえ" << endl;
       break;
   }
}
cout << "はい" << endl;
return 0;
 

這個演算法的時間複雜度是 \(O(\sqrt[]{n})\) 。

程式碼好寫歸好寫，但是時間複雜度真的是太高了，有時候並不能滿足我們的需要；更何況我們經常需要的操作是找出 \(n\)以內的所有質數，而不是簡單地判斷一個數是否是質數。
從\([1,n]\) 範圍內的所有整數裡挑質數其實很簡單。
很多人大概接觸過這樣一種方法：

• 首先列出所有數字並按照從大到小的順序排好：

• 然後拿出五顏六色的彩筆，開始對資料進行處理
• \(1\) 既不是質數也不是合數，提前把它剔除

• 然後看數列的第一個數： \(2\) 。拿筆圈出2，並且劃掉2的所有倍數（也就是所有剩下的偶數）。

• 數列的下一個數是 \(3\)
  。換個顏色的筆圈出3，並劃掉所有剩下的3的倍數。

• 越過被劃掉的4，下一個數是 \(5\) 。操作我就不用說了吧。

• 下一個是 \(7\) 。

• 越過8,9和10，我們來到了 \(11\) 。由於11大於 \(\sqrt[]{100}\) ，也就是10，我們就停止操作，換一個顏色的筆，圈出剩下的所有數。

• 擦掉所有劃去的數，我們得到了100以內的所有質數。

如果我們將這個操作用程式碼來實現，那就是：

埃氏篩

埃氏篩 ，或者叫做 埃拉託斯特尼篩法 ，利用的就是這種思想。
程式碼如下：

int n;
cin >> n;
bool visit[n];
memset(visit, false, n);
for(int i = 2; i <= n; i++)
{
   if(!visit[i])
   {
       if(n < i * i)break;//埃氏篩只需要剔除根號n以內的素數的倍數 
       for(int j = i * i; j <= n; j += i)visit[j] = true;//標記質數
   }
}
for(int i = 2; i <= n; i++)
{
   if(!visit[i])cout << i << endl;//輸出 
}
 

埃氏篩法的時間複雜度只有 \(O(n \log \log n)\) 。
時間比簡單方法的 \(O(\sqrt[]{n})\)，埃氏篩法明顯快了很多。

但是，埃氏篩法還可以進一步優化。
舉個栗子：
數字 \(223092870\) 可以被分解成 \(2 \times 3 \times 5 \times 7 \times 11 \times 13 \times 17 \times 19 \times 23\) ，在埃氏篩法中會被搜尋9次。
如果我們可以在這個數被搜尋到的第一次就標記好，並且使得這個數在之後的過程中不會被再次搜尋到，我們演算法的速度就會快很多。
那麼如何實現這一想法呢？
我們繼續拿 \(223092870\) 來舉例。如果我們在開始判斷之前就用 &1 的方法來去掉所有偶數呢？
顯然會快很多。但是，類似 \(4849845\) ( \(3 \times 5 \times 7 \times 11 \times 13 \times 17 \times 19\) )這樣的數仍然會被不可避免地搜尋到多次。
但是，如果我們保證 \(4849845\) 這個數只在被分解成 \(3 \times 1616615\) 的時候被搜尋到並標記為合數，就可以實現這個想法。

於是，我們就得到了一種演算法：

尤拉篩

尤拉篩 ，又叫 線性篩 ，在埃氏篩的基礎上進行了優化，使得每個數被且只被搜尋到一次。
尤拉篩利用合數的 最小質因子 進行篩選。
再舉一個例子：
數字 \(4849845\) 被搜尋到的路線是這樣的：
\(19\to 19\times 17=323\to 19\times 17\times 13=4199\to 19\times 17\times 13\times 11=46189\to 19\times 17\times 13\times 11\times 7=323323\)
\(\to 19\times 17\times 13\times 11\times 7\times 5=1616615\to 19\times 17\times 13\times 11\times 7\times 5\times 3=4849845\)

如果這還不夠直觀，我們可以打表看一下100以內的合數是怎麼被推出來的。

、

當 \(i>50\) 之後，因為 \(i\times 2\) 已經超出100這個範圍，無法推出任何更多的範圍內的合數。

附上程式碼：

long long n;
cin >> n;
int prime[n];
bool visit[n];
memset(prime, 0, n);
memset(visit, false, n);
prime[0] = 0;
for(long long i = 2; i <= n; i++)
{
   if(!visit[i])
   {
       prime[prime[0] + 1] = i;
       prime[0]++;//記錄質數 
   }
   for(long long j = 1; j <= prime[0] && i * prime[j] <= n; j++)
   {
       visit[i * prime[j]] = true;
       if(i % prime[j] == 0)
       {
           break;//判斷是否是最小質因子 
       }
   }
}
for(long long i = 1; i <= prime[0]; i++)
{
   if(n % prime[i] == 0)
   {
       visit[n] = true;
   }//最後一個一般篩不出來，特判一下 
}
for(long long i = 2; i <= n; i++)
{
   if(!visit[i])
   {
       cout << i << endl;//輸出 
   }
}

尤拉篩 的時間複雜度是 \(O(n)\) （所以這就是為什麼它也叫 線性篩 ），在資料量很大的時候，速度會比埃氏篩快很多。

如果需要更直觀一點的話，可以看這個動畫：

2倍速版：

當 \(n\) 很大的時候，找出\(n\)以內的所有質數已經不太現實了。但是，判斷 \(n\) 是否為質數仍然是一種 （看起來） 可行的操作。
這就用到了：

\(Miller-Rabin\) 定理

\(Miller-Rabin\) 定理涉及了兩個定理：

費馬小定理

費馬小定理是這樣的：

對於質數 \(p\) 和任意整數 \(a\) ，有 \(a^p\equiv a\pmod{p}\)。

這裡我們需要用的是它的逆定理，即

若滿足 \(a^p\equiv a\pmod{p}\) ， \(p\) 為質數。
此時我們兩邊同時約去一個 \(a\) ，得 \(a^{p-1} \equiv 1\pmod{p}\)。

很遺憾，費馬小定理的逆定理不是個真命題，但滿足條件的 \(p\) 還是有很大概率是質數的。
但是如果我們選取不同的 \(a\) 進行多次測試，那麼我們幾乎就可以斷定這個數是一個質數了。

二次探測定理

如果 \(p\) 是奇素數，則 \(x^2\equiv 1\pmod{p}\) 的解為 \(x\equiv 1\) 或 \(x\equiv p-1\pmod{p}\) 。
這是很容易證明的：

又∵ \(p\) 為奇素數，其因子只有 \(1\) 與 \(p\) 兩個。
∴只有兩解 \(x\equiv 1\) 與 \(x\equiv p-1 \pmod{p}\) 。

如果 \(a^{n−1}\equiv 1\pmod{n}\) 成立，\(Miller−Rabin\) 演算法不是立即找另一個 \(a\) 進行測試，而是看 \(n−1\) 是不是偶數。如果 \(n−1\) 是偶數，另 \(u=\dfrac{n-1}{2}\) ，並檢查是否滿足二次探測定理即 \(a^u\equiv 1\) 或 \(a^u\equiv n−1 \pmod{n}\) 。
將這兩條定理合起來，也就是最常見的 Miller−Rabin 測試 。
單次時間複雜度是 \(O(\log n)\) ，總複雜度就是 \(O(times\times \log n)\) 了。

參考程式碼：

#include<cstdio>
#define ll long long 
inline int read()
{
    char c = getchar(); int x = 0, f = 1;
    while(c < '0' || c > '9') { if(c == '-') f = -1; c = getchar(); }
    while(c >= '0' && c <= '9') x = x * 10 + c - '0', c = getchar();
    return x * f;
}
int N, M, Test[10] = { 2, 3, 5, 7, 11, 13, 17 };
int pow(int a, int p, int mod)
{
    int base = 1;
    for(; p; p >>= 1, a = (1ll * a * a) % mod)
        if(p & 1) base = (1ll * base * a) % mod;
    return base % mod;
}
bool Query(int P)
{
    if(P == 1) return 0;
    int t = P - 1, k = 0;
    while(!(t & 1)) k++, t >>= 1;
    for(int i = 0; i < 4; i++)
    {
        if(P == Test[i]) return 1;
        ll a = pow(Test[i], t, P), nxt = a;
        for(int j = 1; j <= k; j++)
        {
            nxt = (a * a) % P;
            if(nxt == 1 && a != 1 && a != P - 1) return 0;
            a = nxt;
        }
        if(a != 1) return 0;
    }
    return 1;
}
int main()
{
    N = read(); M = read();
    while(M--) puts(Query(read()) ? "Yes" : "No");
    return 0;
}

\(\blacktriangleright\)