Pollard-Rho 是一种期望
O
(
n
1
4
)
O(n^{\frac{1}{4}})
O(n41) 时间分解质因数的算法。
Pollard-Rho 用到了 Miller-Rabin,Miller-Rabin 是一种 O ( log n ) O(\log n) O(logn) 时间判素数的方法。
首先来讲 Miller-Rabin,它用到了一种叫二次探测的东西和费马小定理。
Miller-Rabin
二次探测
有这样一个方程 x 2 ≡ 1 ( m o d p ) x^2 \equiv 1 \,(\mathrm{mod}\,p) x2≡1(modp), x x x 的根分别为 1 1 1 和 p − 1 p-1 p−1。
证明:
原方程即为 ( x − 1 ) ( x + 1 ) ≡ 0 ( m o d p ) (x-1)(x+1) \equiv 0\,(\mathrm{mod}\,p) (x−1)(x+1)≡0(modp)。
于是 x = 1 x=1 x=1 或 x = − 1 x=-1 x=−1。
也就是 x = 1 x=1 x=1 或 x = p − 1 x=p-1 x=p−1。
我们可以利用这个性质来判断 x x x 是否为质数。
费马小定理
x p − 1 ≡ 1 ( m o d p ) x^{p-1}\equiv 1\,(\mathrm{mod}\,p) xp−1≡1(modp)
这个不用多说。
快速乘
因为要判断素数的数和要分解质因数的数都很大,乘一次就超出了 long long 范围,所以要用快速乘。
模板:
ll qpow(ll a, ll b, ll p) {
return (a * b - (ll)((long double)a / p * b) * p + p) % p;
}
原理是余数等于被除数减商乘除数。
Miller-Rabin 算法
我们选取了一个质数 p p p,首先对 x x x 判断费马小定理是否成立,不成立就一定不是素数,返回 false。
如果成立,且 2 ∣ ( p − 1 ) 2\mid(p-1) 2∣(p−1)。
那么有 ( x p − 1 2 ) 2 ≡ 1 ( m o d p ) (x^{\frac{p-1}{2}})^2\equiv 1\,(\mathrm{mod}\,p) (x2p−1)2≡1(modp)。
这个式子长得和二次探测很像。
若 x p − 1 2 m o d p = 1 x^{\frac{p-1}{2}}\,\mathrm{mod}\, p=1 x2p−1modp=1,那么就可以继续判断,因为这样式子右边还是 1 1 1,还可以用二次探测。
若 x p − 1 2 m o d p = p − 1 x^{\frac{p-1}{2}}\,\mathrm{mod}\, p=p-1 x2p−1modp=p−1,那么无法判断 x x x 是不是质数,返回 true。
其他情况, x x x 一定不是质数,返回 false。
当 2 ∤ ( p − 1 ) 2\nmid(p-1) 2∤(p−1) 时,停止循环,返回 true。
虽然 Miller-Rabin 不一定正确,但可以多选取几个质数,这样正确率就很高了。
Miller-Rabin 代码如下:
ll qmul(ll a, ll b, ll p) { return (a * b - (ll)((long double)a / p * b) * p + p) % p; }
ll qpow(ll x, ll k, ll p) {
x %= p;
ll res = 1;
while(k) {
if(k & 1) res = qmul(res, x, p);
x = qmul(x, x, p), k >>= 1;
}
return res;
}
const int pcnt = 12;
const ll p0[pcnt + 1] = { 0, 2, 3, 5, 7, 11, 13, 17, 19, 61, 2333, 4567, 24251 };
bool miller_rabin(ll x, ll p) {
if(qpow(x, p - 1, p) != 1) return false;
ll k = p - 1;
while(!(k & 1)) {
k >>= 1;
ll t = qpow(x, k, p);
if(t != 1 && t != p - 1) return false;
if(t == p - 1) return true;
}
return true;
}
bool test(ll p) {
if(p == 1) return false;
lp(i, 1, pcnt) if(p == p0[i]) return true;
lp(i, 1, pcnt) if(!miller_rabin(p0[i], p)) return false;
return true;
}
Pollard-Rho
生日悖论
有 n n n 个人,一年有 d d d 天,存在两个人生日相同的概率是多少。
不难看出答案是 p ( n ) = 1 − d ! ( d − n ) ! n d p(n)=1-\frac{d!}{(d-n)!n^d} p(n)=1−(d−n)!ndd!,若 d d d 为 365,那么 n n n 为 23 23 23 时 p ( n ) p(n) p(n) 就超过了 1 2 \frac{1}{2} 21。
这启示我们如果选出若干个数,将它们两两作差,得到要分解的数 n n n 的因数的概率很大。考虑 n = p 2 n=p^2 n=p2 时(最坏情况)。
设选出 i i i 个数后使其存在两个数在模 p p p 意义下相等的期望步数为 E ( i ) E(i) E(i)。
那么有
E
(
p
)
=
1
E
(
i
)
=
p
−
i
p
(
E
(
i
+
1
)
+
1
)
+
i
p
=
p
−
i
p
E
(
i
+
1
)
+
1
\begin{aligned} E(p) & = 1 \\E(i) & = \frac{p-i}{p}(E(i+1)+1)+\frac{i}{p} \\ &=\frac{p-i}{p}E(i+1)+1 \end{aligned}
E(p)E(i)=1=pp−i(E(i+1)+1)+pi=pp−iE(i+1)+1
E ( 0 ) E(0) E(0) 大概是 p = n 1 4 \sqrt p=n^{\frac{1}{4}} p=n41,怎么得出来的呢?懒得推了,写个程序跑一遍就知道了。
但是这样做依旧不行,因为如果两两之间判断的话复杂度就又回到了 O ( n 1 2 ) O(n^{\frac{1}{2}}) O(n21),不是想要的结果。
伪随机数序列
Pollard-Rho 使用伪随机数生成序列(因为它有个特殊性质),假设现在生成了随机数 x x x,设 f ( x ) = ( x 2 + C ) m o d n f(x)=(x^2+C)\,\mathrm{mod}\,n f(x)=(x2+C)modn,那么 x x x, f ( x ) f(x) f(x), f ( f ( x ) ) f(f(x)) f(f(x))… 构成了伪随机数序列。
它会形成一个形似 状的序列,即从某一点之后形成了循环。
它还有一个很好的性质。
若 x ∣ ∣ a − b ∣ x\mid |a-b| x∣∣a−b∣,那么 f ( a ) − f ( b ) = a 2 − b 2 = ( a + b ) ( a − b ) f(a)-f(b)=a^2-b^2=(a+b)(a-b) f(a)−f(b)=a2−b2=(a+b)(a−b)。
因此 x ∣ ∣ f ( a ) − f ( b ) ∣ x\mid |f(a)-f(b)| x∣∣f(a)−f(b)∣。
这个性质接下来会用到。
Pollard-Rho 算法
终于来了!
基于 Floyd 判环的 Pollard-Rho 算法
Pollard-Rho 通过不断将 x x x 分解为两个因数的乘积,再将这两个因数继续分解,直到分解为质数。
每次执行 Pollard-Rho 时,先用 Miller-Rabin 判断要分解的数是否为质数,若是则说明找到了一个质因子。
Pollard-Rho 使用伪随机数序列,初始先随机出两个数 t 1 t_1 t1 和 t 2 t_2 t2,令 t 1 ← f ( t 2 ) t_1 \leftarrow f(t_2) t1←f(t2),即先让 t 2 t_2 t2 走一步。接下来每次循环 t 1 t_1 t1 走一步, t 2 t_2 t2 走两步。这样当 t 1 = t 2 t_1=t_2 t1=t2 时说明已经走完了一圈。
每次循环,判断是否 g c d ( ∣ t 1 − t 2 ∣ , x ) ≥ 2 gcd(|t_1-t_2|,x)\ge2 gcd(∣t1−t2∣,x)≥2,若是的话,说明找到了两个 x x x 的因子,递归下去继续分解这两个因子。
由于之前说的性质,假设 t 1 t_1 t1 和 t 2 t_2 t2 之间距离为 d d d,判断了 g c d ( ∣ t 1 − t 2 ∣ , x ) gcd(|t_1-t_2|,x) gcd(∣t1−t2∣,x),就相当于判断了序列上其他距离为 d d d 的两个点,而 t 1 t_1 t1 和 t 2 t_2 t2 间的距离每次增加 1 1 1,因此相当于判断了序列上任意两个数。
另外若在一个环上没能分解 x x x,就再随机一次继续找下去。
这样做期望时间复杂度为 O ( n 1 4 log n ) O(n^{\frac{1}{4}}\log n) O(n41logn)。
void pollard_rho(ll x) {
if(test(x)) {
// 此时 x 为其中一个质因子,在这里做你想做的事
return;
}
while(true) {
ll t1 = rd() % (x - 1) + 1, t2 = t1, b = rd() % (x - 1) + 1;
t2 = (qmul(t2, t2, x) + t2 + b) % x;
while(t1 != t2) {
ll t = gcd(abs(t1 - t2), x);
if(t != 1) {
pollard_rho(t), pollard_rho(x / t);
return;
}
t1 = (qmul(t1, t1, x) + t1 + b) % x;
t2 = (qmul(t2, t2, x) + t2 + b) % x;
t2 = (qmul(t2, t2, x) + t2 + b) % x;
}
}
}
倍增优化
我们可以将 ∣ t 1 − t 2 ∣ |t_1-t_2| ∣t1−t2∣ 的结果乘起来,每隔 1,2,4,8…… 次循环求一次 gcd,因为如果 gcd ( t , x ) ≥ 2 \gcd(t,x)\ge2 gcd(t,x)≥2,那么有 gcd ( k t m o d x ) = gcd ( ( k m o d x ) ( t m o d x ) , x ) ≥ 2 \gcd(kt\,\mathrm{mod}\,x)=\gcd((k\,\mathrm{mod}\,x)(t\,\mathrm{mod}\,x),x)\ge2 gcd(ktmodx)=gcd((kmodx)(tmodx),x)≥2。
通常我们可以设一个最大距离,倍增超过这个距离后就一直保持这个距离。
特别的,当累乘的结果要变成 0 0 0 时(因为是在模 x x x 意义下的),就退出累乘,直接计算公因数。
期望时间复杂度 O ( n 1 4 ) O(n^{\frac{1}{4}}) O(n41) 。
void pollard_rho(ll x) {
if(test(x)) return;
while(true) {
ll t1 = rd() % (x - 1) + 1, b = rd() % (x - 1) + 1, t2 = F(t1, b, x);
int lim = 1;
while(t1 != t2) {
ll cnt = 1;
lp(i, 1, lim) {
ll tmp = qmul(cnt, abs(t1 - t2), x);
if(!tmp) break;
cnt = tmp;
t1 = F(t1, b, x), t2 = F(F(t2, b, x), b, x);
}
ll t = gcd(cnt, x);
if(t != 1) return pollard_rho(t), pollard_rho(x / t), void();
lim = min(lim << 1, 128);
}
}
}
模板代码:
ll qmul(ll a, ll b, ll p) { return (a * b - (ll)((long double)a / p * b) * p + p) % p; }
ll qpow(ll x, ll k, ll p) {
x %= p;
ll res = 1;
while(k) {
if(k & 1) res = qmul(res, x, p);
x = qmul(x, x, p), k >>= 1;
}
return res;
}
const int pcnt = 12;
const ll p0[pcnt + 1] = { 0, 2, 3, 5, 7, 11, 13, 17, 19, 61, 2333, 4567, 24251 };
bool miller_rabin(ll x, ll p) {
if(qpow(x, p - 1, p) != 1) return false;
ll k = p - 1;
while(!(k & 1)) {
k >>= 1;
ll t = qpow(x, k, p);
if(t != 1 && t != p - 1) return false;
if(t == p - 1) return true;
}
return true;
}
bool test(ll p) {
if(p == 1) return false;
lp(i, 1, pcnt) if(p == p0[i]) return true;
lp(i, 1, pcnt) if(!miller_rabin(p0[i], p)) return false;
return true;
}
ll gcd(ll a, ll b) { return a ? gcd(b % a, a) : b; }
mt19937 rd(time(NULL));
ll F(ll x, ll b, ll p) { return (qmul(x, x, p) + x + b) % p; }
void pollard_rho(ll x) {
if(test(x)) return;
while(true) {
ll t1 = rd() % (x - 1) + 1, b = rd() % (x - 1) + 1, t2 = F(t1, b, x);
int lim = 1;
while(t1 != t2) {
ll cnt = 1;
lp(i, 1, lim) {
ll tmp = qmul(cnt, abs(t1 - t2), x);
if(!tmp) break;
cnt = tmp;
t1 = F(t1, b, x), t2 = F(F(t2, b, x), b, x);
}
ll t = gcd(cnt, x);
if(t != 1) return pollard_rho(t), pollard_rho(x / t), void();
lim = min(lim << 1, 128);
}
}
}
模板题:洛谷 [P4718]