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\) 和 \(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\) 平方而来。所以我们只要平方，然后如果没有遇到 \(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\)，\(seed\) 是随机选择的一个数。
\(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;
}
浙公网安备 33010602011771号