Pollard-Rho演算法 學習筆記
Miller Rabin 素性測試
假設我們要判定一個數 \(m\) 是否是素數。
Fermat素性測試
可以隨便找到一個數 \(a\)，然後判定 \(a^{m - 1} \bmod m\) 是否為 \(1\)。
如果我們找到 \(k\) 個 \(a\) 都滿足 \(a^{m - 1} \equiv 1\)，那麼我們就把 \(m\) 判定為素數。
這樣子的判定事實上是不對的，因為為有很多 (無限個) 合數滿足對於任意 \(a\)， \(a^{m - 1} \equiv 1 \pmod m\)
二次探測定理
對於任意奇素數 \(p\)， \(x^2 \equiv 1 \pmod p\) 的解只有 \(1\)
證明：\(x^2 - 1 \equiv 0 \pmod p\)，\((x + 1) (x - 1) \pmod p\)，因此要麼 \(p | x + 1\)，要麼 \(p | x - 1\)，因此 \(x\) 為 \(1\) 或 \(p - 1\)。
我們把二次探測定理和費馬小定理配合起來使用。
把 \(n - 1\) 分解成 \(t \times 2^k\)。我們判一下 \(k = a^t \pmod m\)。如果不是 \(1\) 且不是 \(m - 1\)，那麼繼續判斷。
如果 \(m\) 是素數，那麼平方之後為 \(1\) 且原先不是 \(1\) ，那麼一定是由 \(m - 1\)
一次成功判定的概率是 \(\frac{3}{4}\)。
（程式碼見下面 Pollard-Rho 演算法 的）
Pollard-Rho 演算法
生日悖論
隨機找 \(23\) 個人，他們中有人在同一天生日的概率是 \(50 \%\)
推廣：隨機寫在 \([1, n]\) 中的數，它們中有兩個數相同的期望次數是 \(\sqrt n\)
演算法
（假設我們要找 \(n\) 的因子）
我們生成一個序列 \(F\)，\(f_i = F(f_{i - 1})\)。\(F(x)\) 是一個函式，這裡 \(F(x) = ( x^2 + seed ) \bmod n\)
\(F(x)\) 是幾乎隨機的，\(F(x) \bmod q\) 和 \(F(x) \bmod n\) 也幾乎隨機。
假設 \(n\) 有一個因數 \(q\) ，那麼 \(F(x) \bmod q\) 必然會重複，因此有迴圈節，且迴圈節長度期望為 \(\sqrt q\)。
如果找到迴圈節，那麼用迴圈節中重複的兩個數 \(val1, val2\)，\(val1 \equiv val2 \pmod q\)，我們可以用 \(\gcd(n, |val1 - val2|)\) 來找到 \(q\)。
如何找迴圈節？
可以設兩個變數 \(a\) 和 \(b\) (初值相等)，每次讓 \(a = F(a), b = F(F(b))\)。最後他們會相遇，這樣子我們就找到了環。
找環的期望步數是 \(\sqrt q\) 的。\(q\) 可以取到 \(n\) 的最小因數，所以步數是期望 \(n^{\frac{1}{4}}\)。
但是我們找環的時候每次都要做一遍 \(\gcd\)，時間複雜度 \(\Theta(n^{\frac{1}{4}} \log n)\) 。
但是我們可以優化：分段處理。每 \(lim\) 個數處理一次，處理出 \(mul = \prod\limits_{i = 1}^{lim} p_i \bmod n\)，然後判一下 \(\gcd(mul, n)\) 是否為 \(1\)。
如果 \(\gcd\) 不為 \(1\) 就判一下找到環的位置。時間複雜度變成了 \(\Theta (n^{\frac{1}{4}})\) 。
程式碼：
#include<bits/stdc++.h>
#define L(i, j, k) for(int i = j, i##E = k; i <= i##E; i++)
#define R(i, j, k) for(int i = j, i##E = k; i >= i##E; i--)
#define ll long long
#define ull unsigned long long
#define db double
#define pii pair<int, int>
#define mkp make_pair
using namespace std;
mt19937_64 orz(233);
const int N = 1e5 + 7;
ll Sum(ll aa, ll bb, ll mod) {
aa += bb;
return aa >= mod ? aa - mod : aa;
}
ll Mul(ll a, ll b, ll mod) {
ll sav = (a * b - (ll)((long db) a * b / mod + 0.1) * mod) % mod;
return sav < 0 ? sav + mod : sav;
}
ll mypow(ll x, ll y, ll mod) {
ll res = 1;
for(; y; x = Mul(x, x, mod), y >>= 1) if(y & 1) res = Mul(res, x, mod);
return res;
}
const int testp[8] = {2, 3, 5, 7, 13, 19, 61, 233};
bool Miller_Rabin(ll x) {
if(x < 3) return x == 2;
if(x <= testp[7]) {
L(i, 2, sqrt(x)) if(x % i == 0) return 0;
return 1;
}
int k, j; ll a;
for(k = 0, a = x - 1; ! (a & 1); k ++) a >>= 1;
L(i, 0, 7) {
ll v = mypow(testp[i], a, x);
if(v == 1 || v == x - 1) continue;
for(j = 0; j < k; j ++) {
v = Mul(v, v, x);
if(v == x - 1) break;
}
if(j == k) return 0;
}
return 1;
}
int total;
ll f[N], seed;
ll F(ll x, ll mod) {
return Sum(Mul(x, x, mod), seed, mod);
}
ll Pollard_Rho(ll x) {
total = 0, seed = orz() % (x - 1) + 1;
ll val1 = rand() % (x - 1) + 1, val2 = val1, mul = 1, sav;
while(1) {
val1 = F(val1, x), val2 = F(F(val2, x), x);
f[++total] = abs(val1 - val2);
mul = Mul(mul, f[total], x);
if(total == 127) {
if(__gcd(mul, x) > 1) {
L(i, 1, total) if((sav = __gcd(f[i], x)) > 1) return sav;
return x;
}
total = 0;
}
}
}
int tot;
ll ans[N];
void edsolve(ll x) {
if(x < 2) return;
if(Miller_Rabin(x)) return ans[++tot] = x, void();
ll sav = Pollard_Rho(x);
while(sav == x) sav = Pollard_Rho(x);
edsolve(sav), edsolve(x / sav);
}
void solve(ll x) {
tot = 0;
L(i, 2, testp[7]) while(x % i == 0) ans[++tot] = i, x /= i;
edsolve(x), sort(ans + 1, ans + tot + 1);
}
ll n;
void Main() {
scanf("%lld", &n), tot = 0;
if(Miller_Rabin(n)) puts("Prime");
else solve(n), printf("%lld\n", ans[tot]);
}
int main() {
int T; scanf("%d", &T);
while(T--) Main();
return 0;
}