【模板】Pollard-Rho算法
题目链接:luogu P4718
题目大意
要你求一个数的最大质因子。
思路
Miller Rabin
首先,这种题肯定要判断一个数是不是素数。
但普通的方法会爆炸。
于是我们要用 Miller Rabin。
这是一个 O ( l o g n ) O(logn) O(logn) 判断质数的东西,运用了费马小定律和二次探测定理。
费马小定理
a p − 1 ≡ 1 ( m o d p ) a^{p-1}\equiv1(\bmod\ p) ap−1≡1(mod p)
(条件是 p p p 是质数)
那我们可以根据这个来判断,随机几次 a a a,看是否成立,如果不成立就说明不行。
但也会有一些数,它是合数,但你就是找不到 a a a 使得式子不满足,这些数叫做 Carmichael 数,就比如 561 561 561。
这些数会使得判定出现误差,那我们考虑怎么搞。
二次探测定理
若 b 2 ≡ 1 ( m o d p ) b^2\equiv1(\bmod\ p) b2≡1(mod p), p p p 为质数,则 p p p 一定可以被 b + 1 b+1 b+1 和 b − 1 b-1 b−1 其中一个整除。
证明其实很简单:
b 2 ≡ 1 ( m o d p ) b^2\equiv1(\bmod\ p) b2≡1(mod p)
b 2 − 1 ≡ 0 ( m o d p ) b^2-1\equiv0(\bmod\ p) b2−1≡0(mod p)
( b + 1 ) ( b − 1 ) ≡ 0 ( m o d p ) (b+1)(b-1)\equiv0(\bmod\ p) (b+1)(b−1)≡0(mod p)
然后就很显然了, p p p 是质数,所以就一定有倍数关系了。
然后再想想就发现 b = 1 / p − 1 b=1/p-1 b=1/p−1。
然后如果 b = 1 b=1 b=1 那又可以继续二次探测!
那转回到费马小定理:
a p − 1 ≡ 1 ( m o d p ) a^{p-1}\equiv1(\bmod\ p) ap−1≡1(mod p)
a ( p − 1 2 ) 2 ≡ 1 ( m o d p ) a^{{(\frac{p-1}{2})}^2}\equiv1(\bmod\ p) a(2p−1)2≡1(mod p)(前提是 p − 1 p-1 p−1 是偶数)
然后就搞,然后就按这个来二次探测就可以了。
这么搞了之后,就很难会出错,大概 1 0 10 10^{10} 1010 才会有一个出错,所以可以忽略不计。
然后不难看到要用点快速幂快速乘。
Pollard-Rho 算法
但你只能判断质数还不行,你还要分解质因数啊。
这时候就要用 Pollard-Rho 算法了,可以快速的将一个你确实是合数的数找到它的其中一个因子。
GCD
我们先从 GCD 来看,因为是合数,所以必然会有小于 n \sqrt{n} n
的因子。那我们考虑随机一个数,看它与 n n n 的 gcd \gcd gcd 是否大于 1 1 1,得到的 gcd \gcd gcd 就是一个因子。那这样选到的可能就是 n \sqrt{n} n
。
但在 1 0 18 10^{18} 1018 面前,还是太弱小了,没有我们要的速度。
x^2+c
这是一个玄学的随机方法,可以使得你找到跟 n n n 有非 1 1 1 公约数的概率会增加。
这个方法是利用了生日悖论,它不符合生活经验,但是可以推导出来没问题的东西。
它的 x , c x,c x,c 一开始是随机的,然后让 y = x , x = x 2 + c y=x,x=x^2+c y=x,x=x2+c,然后每次拿 y − x y-x y−x 去跟 p p p 求 gcd \gcd gcd。
它能搞到的概率十分之高,只能用玄学来形容。
判环
它这种方法很优秀,但是你因为生成过程中会取模,搞了多次一定会有环,那你如果这样都没找到,那你一直找下去都不会有结果,还不如重新随机 x , c x,c x,c,重新开搞。
那你就要判环。
我们可以用倍增的思想, y y y 记住 x x x 位置, x x x 跑当前数量的一倍的次数,然后再记住,再跑。
因为有倍增,那 x , y x,y x,y 位置相同时, x x x 一定跑到了圈,而且不会多跑太多。
二次优化
然后你发现你每次判断都要求 gcd \gcd gcd,这会让复杂度由 O ( n 1 4 ) O(n^{\frac{1}{4}}) O(n41) 变成 O ( n 1 4 log n ) O(n^{\frac{1}{4}}\log n) O(n41logn)。
虽然 log n \log n logn 很小,但毕竟是 1 0 18 10^{18} 1018,会让整个变慢一些。
你发现取模的数不是质数。然后有这么一个质数:
如果模数和被模的数都有一个公约数,那取模结果也是这个公约数的倍数。
那我们可以把中途的 y − x y-x y−x 都乘起来,那只要其中有一个跟 n n n 有非 1 1 1 共约数,你乘起来的数跟 n n n 就有非 1 1 1 公约数。
那我们可以多算几次 y − x y-x y−x,再求一次 gcd \gcd gcd,经过计算和测试,一般连续算 127 127 127 次再 gcd \gcd gcd 一次是最优的,跑的很快。
但你会发现会有一点小问题。
我们可能跑了没有 127 127 127 次,就出现了环。然后就导致还没有求出答案,就退了出来。
或者由于算法过度优秀(bushi,乘积乘着乘着直接变 n n n 倍数,取模直接是 0 0 0,那没有必要干等着。
那处理第一个问题,我们可以用倍增来解决,在每个 i = j i=j i=j 要重新换位置的时候我们也做一次 gcd \gcd gcd,因为是倍增的时候才 gcd \gcd gcd 所以影响不大。
那第二个简单的也做一个判断就可以了。
然后就好啦。
代码
#include<cstdio>
#include<cstdlib>
#define ll long long
using namespace std;
int T;
ll n, maxn_pr;
ll Abs(ll x) {
return (x < 0) ? -x : x;
}
ll ksc(ll x, ll y, ll mo) {
x %= mo; y %= mo;
ll c = (long double)x * y / mo;
long double re = x * y - c * mo;
if (re < 0) re += mo;
else if (re >= mo) re -= mo;
return re;
}
ll ksm(ll x, ll y, ll mo) {
ll re = 1;
while (y) {
if (y & 1) re = ksc(re, x, mo);
x = ksc(x, x, mo);
y >>= 1;
}
return re;
}
ll gcd(ll x, ll y) {
ll tmp;
while (y) {
tmp = x; x = y; y = tmp % y;
}
return x;
}
bool mr(ll x, ll p) {
if (ksm(x, p - 1, p) != 1) return 0;
ll y = p - 1, z;
while (!(y & 1)) {
y >>= 1;
z = ksm(x, y, p);
if (z != 1 && z != p - 1) return 0;
if (z == p - 1) return 1;
}
return 1;
}
bool prime(ll now) {//Miller Rabin 判是否是质数
if (now < 2) return 0;
if (now == 2 || now == 3 || now == 5 || now == 7 || now == 43) return 1;
return mr(2, now) && mr(3, now) && mr(5, now) && mr(7, now) && mr(43, now);
}
ll rho(ll p) {
ll x, y, z, c, g;
int i, j;
while (1) {
x = y = rand() % p;
z = 1;
c = rand() % p;
i = 0; j = 1;
while (++i) {
x = (ksc(x, x, p) + c) % p;
z = ksc(z, Abs(y - x), p);
if (x == y || !z) break;
if (!(i % 127) || i == j) {//优化
g = gcd(z, p);
if (g > 1) return g;
if (i == j) {
y = x;
j <<= 1;
}
}
}
}
}
void work(ll now) {
if (now < 2) return ;
if (now <= maxn_pr) return ;//小小剪枝
if (prime(now)) {
maxn_pr = now;
return ;
}
ll l = rho(now);
while (now % l == 0) now /= l;//要一直除,把这个因子都除掉
work(l); work(now);//继续分解出来继续搞
}
int main() {
srand(19491001);
scanf("%d", &T);
while (T--) {
scanf("%lld", &n);
maxn_pr = 0;
work(n);
if (maxn_pr == n) printf("Prime\n");
else printf("%lld\n", maxn_pr);
}
return 0;
}