CHANGE LOG
- 2022.7.13：重構文章，更新 PR 模板代碼。
- 2023.1.23：對文章進行修補。
1. Miller-Rabin
Miller-Rabin 素性測試是一種具有隨機性的素數判定方法。它有一定概率將合數判定為素數，但不會將素數判定為合數。
素數判定的基本思路為根據所有質數但很少合數具有的性質，檢查被判定的數是否具有這些性質。若不具有，則該數是合數，否則該數大概率是質數。
1.1 費馬素性檢驗
當 \(p\) 是素數時，對於任意 \(a \perp p\) 均有 \(a ^ {p - 1}\equiv 1 \pmod p\)。相反，當 \(a ^ {p - 1} \equiv 1\pmod p\) 時，是否有 \(p\) 是素數？
可惜命題並不成立。有極小概率使得 \(a \perp p\)，\(p\) 是合數且 \(a ^ {p - 1} \equiv 1\pmod p\)，如當 \(a = 2\)，\(p = 341\) 時，\(2 ^ {340} \equiv 1 \pmod {341}\)，稱 \(341\) 是以 \(2\) 為底的偽素數。\(341\) 是最小的偽素數基數。
若 \(a ^ {p - 1}\not \equiv 1\pmod p\)，則 \(p\) 必然不是質數。多選幾個與 \(p\) 互質的數檢查，可以排除大部分合數。這被稱為費馬素性檢驗，它具有隨機性。
當面對形如 \(561\) 的卡邁克爾數時，費馬素性檢驗就相當劣了。卡邁克爾數 \(p\) 滿足所有與 \(p\) 互質的數的 \(p - 1\) 次方模 \(p\) 均為 \(1\)。
1.2 二次探測定理
根據二次剩余部分的知識，當 \(p\) 為奇素數時，\(x ^ 2\equiv 1\pmod p\) 有且僅有解 \(\pm 1\)。因此，若存在 \(a \neq \pm 1\) 滿足 \(a ^ 2 \equiv 1\pmod p\)，則 \(p\) 必然不是質數。
這被稱為 二次探測定理。
1.3 算法介紹
結合費馬素性檢驗與二次探測定理。
根據二次探測定理，當 \(a ^ {p - 1} \equiv 1\pmod p\) 時，若 \(p - 1\) 是 \(2\) 的倍數，則 \(a ^ {\frac {p - 1} 2}\) 必須是 \(\pm 1\)。若 \(\dfrac {p - 1} 2\) 仍是 \(2\) 的倍數且 \(a ^ {\frac {p - 1} 2}\equiv 1\pmod p\)，則 \(a ^ {\frac {p - 1} 4}\) 必須是 \(\pm 1\)，以此類推。
一般地，若 \(p - 1\) 是 \(2 ^ r\) 的倍數且 \(a ^ {\frac {p - 1}{2 ^ {r - 1}}} \equiv 1\pmod p\)，則 \(a ^ {\frac {p - 1}{2 ^ r}}\) 等於 \(\pm 1\)，否則 \(p\) 不是素數。例如 \(a = 2\) 且 \(p = 341\)。\(2 ^ {340} \equiv 1 \pmod {341}\)，\(2 ^ {170} \equiv 1\pmod {341}\)，但 \(2 ^ {85} \equiv 32\pmod {341}\)。這說明 \(341\) 不是質數。
這樣檢測的准確率很高。隨機選擇 \(k\) 個底數，Miller-Rabin 算法的正確率（不會將合數誤判成素數的概率）大於 \(1 - 4 ^ {-k}\)。算法時間復雜度為 \(\mathcal{O}(k \log ^ 2 n)\)。
Miller-Rabin 的效率與選取底數個數有關，我們希望減少底數並保證一定正確性。以下是常用底數，來自 wangrx 的博客。
- 對 \(2 ^ {32}\) 以內的數判素，使用 \(2, 7, 61\) 三個底數。
- 對 \(2 ^ {64}\) 以內的數判素，使用 \(2, 325, 9375, 28178, 450775, 9780504, 1795265022\) 七個底數。
- 使用前 \(12\) 個素數作為底數可對 \(318665857834031151167460(3\times 10 ^ {23} \approx 2 ^ {78})\) 以內的數判素。詳見 A014233 - OEIS。
- 固定底數時需特判底數。若以 \(2, 7, 61\) 作為底數，則當 \(n = 2, 7\) 或 \(61\) 時直接通過檢驗，因為 \(p\) 無法通過以 \(p\) 為底的費馬素性檢驗。
1.4 復雜度優化
Miller-Rabin 的復雜度和正確性足夠優秀，但注意到整個過程中我們多次使用快速冪計算 \(a\) 的冪，且指數每次除以 \(2\)，很浪費。
考慮將整個過程反過來，即預先處理好 \(p - 1 = r \cdot 2 ^ d\)，計算 \(a ^ r\) 並執行 \(d\) 次平方操作，即可得到每個 \(a ^ {r \cdot 2 ^ i}\)。時間復雜度優化為 \(\mathcal{O}(k\log n)\)。
進一步地，先判掉 \(a ^ r\equiv 1\pmod p\)，此時 \(p\) 通過檢驗。否則若 \(p\) 是質數，說明在 \(i\) 從 \(d\) 減小到 \(0\) 的過程中，\(a ^ {r \cdot 2 ^ i}\bmod p\) 存在從 \(1\) 變為 \(-1\) 的過程，因此只需判斷是否存在 \(0\leq i < d\) 使得 \(a ^ {r \cdot 2 ^ i} \equiv -1 \pmod p\)。容易證明這是 \(p\) 通過本輪素性檢驗的充要條件。
代碼見 2.3 小節 P4718。
2. Pollard-Rho
分解質因數一般使用的試除法時間復雜度為 \(\mathcal{O}(\sqrt n)\)，因為必須枚舉到 \(\sqrt n\) 才能確定 \(n\) 為質數。若預先篩出 \(\sqrt n\) 以內的質數，則復雜度除以 \(\log n\)。
當 \(n\) 較小但數據組數 \(T\) 較多時，可預先使用線性篩 \(\mathcal{O}(n)\) 篩出值域內所有數的最小質因子，單次分解質因數 \(\mathcal{O}(\log n)\)。
Pollard-Rho 為時間復雜度又開了一次平方，它可以在期望 \(\sqrt[4]{n}\log n\) 的時間復雜度內求出 \(n\) 的一個非平凡因子，因此使用 Pollard-Rho 分解質因數的時間復雜度為 \(\sqrt [4]{n}\log ^ 2 n\)。
2.1 生日悖論
從 \(1\sim n\) 的正整數中 \(k\) 次等概率隨機選擇一個數，則所有數互不相同的概率為
公式含義為從 \(n\) 個數中有序選擇 \(k\) 個互不相同的數的方案數除以總方案數。
- 直觀認知：當 \(k\) 增大時，\(P\) 衰減很快，因為每次 \(\dfrac {n - k + 1} n\) 以相乘的方式作用在 \(P\) 上。可以理解為指數衰減，但比指數衰減慢。如下圖。
手玩函數圖像后，我們發現使得 \(P = \dfrac 1 2\) 的 \(k\) 似乎是 \(\sqrt n\) 級別的。
根據 \(1 + x \leq e ^ x\)，可知 \(P \leq e ^ {\sum_{i = 1} ^ {k} \frac {n - i + 1} n - 1} = e ^ {-\frac{k(k - 1)}{2n}}\)。令 \(e ^ {-\frac{k(k - 1)}{2n}} = \dfrac 1 2\)，解得 \(k\) 在 \(\sqrt {n\ln 4}\) 附近。
因此，不嚴謹地，在 \(n\) 個數中等概率隨機選擇，使得選出的數中存在兩個相同的數的期望次數為 \(\mathcal{O}(\sqrt n)\)。
2.2 算法介紹
首先，我們必須有能力快速判斷待分解的數是否為質數：Miller-Rabin。
Pollard-Rho 算法的精髓在於構造偽隨機函數 \(f(x) = x ^ 2 + c\)。
因 \(f\) 僅含一個變量，故對於相同的 \(x\)，\(f(x)\) 的返回值相同。在 \(f\) 進入循環前，它不斷迭代得到的數可視為隨機。其隨機性尚未被證明。
根據生日悖論，模 \(n\) 意義下，若給定 \(c\) 和初始值 \(x_0\)，則 \(f\) 在不斷迭代的過程中期望迭代 \(\mathcal{O}(\sqrt {n})\) 次進入長度為 \(\mathcal{O}(\sqrt {n})\) 的循環。一般以 \(0\) 作為初始值，則迭代過程形如 \(x_1 = f(x_0) = c\)，\(x_2 = f(x_1) = c ^ 2 + c\)，\(\cdots\)，\(x_i = f ^ i(x_0)\)。
因為整條路徑類似希臘字母 \(\rho\)，算法得名 Pollard-Rho，如下圖。
若 \(n\) 為合數，則其最小非平凡因子 \(m\) 不超過 \(\sqrt n\)。因此模 \(m\) 意義下期望迭代 \(\mathcal{O}(\sqrt [4]{n})\) 次進入長度為 \(\mathcal{O}(\sqrt [4]{n})\) 的循環。
同時，若模 \(m\) 一旦進入循環即存在 \(x_i \equiv x_j\pmod m\)，即可通過計算 \(|x_i - x_j|\) 和 \(n\) 的 \(\gcd\) 求出 \(m\) 或其倍數，前提為 \(x_i \not\equiv x_j\pmod n\)，否則 \(\gcd\) 結果為 \(n\)，平凡。
記 \(x_i\) 在模 \(n\) 意義下形成的路徑為 \(\rho_n\)，在模 \(m\) 意義下形成的路徑為 \(\rho_m\)。問題為求出一組 \(i, j\) 使得其處於 \(\rho_m\) 同一點，但處於 \(\rho_n\) 不同點。這說明 \(i\) 需要足夠大使得跳出 \(\rho_m\) 的尾巴，且 \(j - i\) 恰好是 \(\rho_m\) 循環節的倍數。據分析，最小的 \(i\) 和 \(j\) 的級別均為 \(\mathcal{O}(\sqrt [4]{n})\)。
注意，整個過程中我們 不知道 \(m\)，但分析可證求 \(m\) 的期望復雜度為 \(\mathcal{O}(\sqrt [4]{n}\log n)\)。求得的非平凡因子 \(m\) 不一定最小，因此還需繼續分解。
考慮從 \(i = 1\) 開始計算 \(d = \gcd(|x_{2i} - x_i|, n)\) 直到該值不等於 \(1\)，此時有兩種情況：
- \(d < n\)，說明 \(d\) 即 \(n\) 的非平凡因子，直接返回。
- \(d = n\)，說明進入 \(\rho_n\) 的循環節，大概率因為 \(\rho_n\) 的環長等於 \(\rho_m\) 的環長，也可能因為 \(\rho_m\) 的尾巴過長，使得第一次枚舉到 \(\rho_m\) 環長的倍數使得 \(i\) 跳出 \(\rho_m\) 的尾巴時就枚舉到了 \(\rho_n\) 的環長。無論如何，應結束本次失敗的 Pollard-Rho，調整參數 \(c\) 重新分解。
上述算法稱為基於 Floyd 判環的 Pollard-Rho 算法，期望時間復雜度 \(\mathcal{O}(\sqrt [4]{n} \log n)\)。
- 注意：筆者實現 Floyd 判環后，發現若令 \(d = \gcd(|x_{2i} - x_i|, n)\) 則 \(n = 4\) 無論 \(c\) 取何值均無法分解。需要特判 \(n = 4\) 或從 \(i = 0\) 開始計算 \(d = \gcd(|x_{2i + 1} - x_i|, n)\)。
- 優化：\(\gcd\) 次數過多降低效率。朴素 Floyd 判環法無法通過模板題。考慮設置 樣本累計上限 \(T\)，將 \(T\) 組 \(\gcd\) 測試打包，將常數減小 \(T\) 倍。\(T = 10\) 在模板題數據下表現優秀。其正確性基於 \(\gcd(a, n) \mid \gcd(ab\bmod n, n)\)。
- 優化：二分未知上界的數的最好方法是倍增。\(n\) 較大時 \(T = 10\) 太小了。考慮每檢查一次就將 \(T\) 乘以 \(2\)。同時為防止 \(T\) 過大無法及時檢查，令 \(T\) 對 \(C\) 取較小值。一般取 \(C = 127\)。
如果讀者擔心 \(n\) 較小時算法的正確性，可預處理答案。讀者需要認識到 PR 本身是隨機性算法的事實，沒有絕對正確的寫法，只有效率和正確性相對優秀的寫法。對於廣為流傳的 PR 寫法，前人已經驗證其在一定范圍內的正確性，可放心大膽使用。
總之，\(d\) 的計算方式多種多樣，但樣本累計優化不可少。檢查太耗時則打包，這樣的思想不僅可用於優化 Pollard-Rho，也可以應用在其它領域。
2.3 例題
P4718【模板】Pollard-Rho 算法
#include <bits/stdc++.h>
using namespace std;
using ll = long long;
mt19937 rnd(chrono::steady_clock::now().time_since_epoch().count());
ll rd(ll l, ll r) {return rnd() % (r - l + 1) + l;}
ll ksm(ll a, ll b, ll p) {
ll s = 1;
while(b) {
if(b & 1) s = (__int128) s * a % p;
a = (__int128) a * a % p, b >>= 1;
}
return s;
}
bool Miller(ll n) {
if(n < 3 || n % 2 == 0) return n == 2;
ll r = n - 1, d = 0;
while(r & 1 ^ 1) r >>= 1, d++;
for(int _ = 0; _ < 10; _++) {
ll a = rd(2, n - 1), v = ksm(a, r, n);
if(v == 1) continue;
for(int i = 0; i <= d; i++) {
if(i == d) return 0;
if(v == n - 1) break;
v = (__int128) v * v % n;
}
}
return 1;
}
ll Pollard(ll n) {
ll c = rd(1, n - 1), s = c, t = 0;
auto f = [&](ll x) {return ((__int128) x * x + c) % n;};
ll acc = 0, prod = 1, d, limit = 1;
while(s != t) {
prod = (__int128) prod * abs(s - t) % n;
if(++acc == limit) {
if((d = __gcd(prod, n)) > 1) return d;
acc = 0, limit = min(127ll, limit << 1);
}
s = f(f(s)), t = f(t);
}
if((d = __gcd(prod, n)) > 1) return d;
return n;
}
ll mxp(ll n) {
if(Miller(n)) return n;
if(n == 1) return 1;
ll d = Pollard(n);
while(d == n) d = Pollard(n);
while(n % d == 0) n /= d;
return max(mxp(d), mxp(n));
}
int main() {
ios::sync_with_stdio(0);
ll T, n;
cin >> T;
while(T--) {
cin >> n;
if(Miller(n)) cout << "Prime\n";
else cout << mxp(n) << "\n";
}
return 0;
}
參考資料
第一章：
第二章：