【luogu P4718】【模板】Pollard-Rho算法（数学）（也含 Miller Rabin）

【模板】Pollard-Rho算法

题目链接：luogu P4718

题目大意

要你求一个数的最大质因子。

思路

Miller Rabin

首先，这种题肯定要判断一个数是不是素数。
但普通的方法会爆炸。
于是我们要用 Miller Rabin。

这是一个 \(O(logn)\) 判断质数的东西，运用了费马小定律和二次探测定理。

费马小定理

\(a^{p-1}\equiv1(\bmod\ p)\)
（条件是 \(p\) 是质数）

那我们可以根据这个来判断，随机几次 \(a\)，看是否成立，如果不成立就说明不行。
但也会有一些数，它是合数，但你就是找不到 \(a\) 使得式子不满足，这些数叫做 Carmichael 数，就比如 \(561\)。
这些数会使得判定出现误差，那我们考虑怎么搞。

二次探测定理

若 \(b^2\equiv1(\bmod\ p)\)，\(p\) 为质数，则 \(p\) 一定可以被 \(b+1\) 和 \(b-1\) 其中一个整除。
证明其实很简单：
\(b^2\equiv1(\bmod\ p)\)
\(b^2-1\equiv0(\bmod\ p)\)
\((b+1)(b-1)\equiv0(\bmod\ p)\)
然后就很显然了，\(p\) 是质数，所以就一定有倍数关系了。
然后再想想就发现 \(b=1/p-1\)。
然后如果 \(b=1\) 那又可以继续二次探测！

那转回到费马小定理：
\(a^{p-1}\equiv1(\bmod\ p)\)
\(a^{{(\frac{p-1}{2})}^2}\equiv1(\bmod\ p)\)（前提是 \(p-1\) 是偶数）
然后就搞，然后就按这个来二次探测就可以了。

这么搞了之后，就很难会出错，大概 \(10^{10}\) 才会有一个出错，所以可以忽略不计。
然后不难看到要用点快速幂快速乘。

Pollard-Rho 算法

但你只能判断质数还不行，你还要分解质因数啊。
这时候就要用 Pollard-Rho 算法了，可以快速的将一个你确实是合数的数找到它的其中一个因子。

GCD

我们先从 GCD 来看，因为是合数，所以必然会有小于 \(\sqrt{n}\) 的因子。那我们考虑随机一个数，看它与 \(n\) 的 \(\gcd\) 是否大于 \(1\)，得到的 \(\gcd\) 就是一个因子。那这样选到的可能就是 \(\sqrt{n}\)。
但在 \(10^{18}\) 面前，还是太弱小了，没有我们要的速度。

x^2+c

这是一个玄学的随机方法，可以使得你找到跟 \(n\) 有非 \(1\) 公约数的概率会增加。
这个方法是利用了生日悖论，它不符合生活经验，但是可以推导出来没问题的东西。

它的 \(x,c\) 一开始是随机的，然后让 \(y=x,x=x^2+c\)，然后每次拿 \(y-x\) 去跟 \(p\) 求 \(\gcd\)。
它能搞到的概率十分之高，只能用玄学来形容。

判环

它这种方法很优秀，但是你因为生成过程中会取模，搞了多次一定会有环，那你如果这样都没找到，那你一直找下去都不会有结果，还不如重新随机 \(x,c\)，重新开搞。

那你就要判环。
我们可以用倍增的思想，\(y\) 记住 \(x\) 位置，\(x\) 跑当前数量的一倍的次数，然后再记住，再跑。
因为有倍增，那 \(x,y\) 位置相同时，\(x\) 一定跑到了圈，而且不会多跑太多。

二次优化

然后你发现你每次判断都要求 \(\gcd\)，这会让复杂度由 \(O(n^{\frac{1}{4}})\) 变成 \(O(n^{\frac{1}{4}}\log n)\)。
虽然 \(\log n\) 很小，但毕竟是 \(10^{18}\)，会让整个变慢一些。

你发现取模的数不是质数。然后有这么一个质数：
如果模数和被模的数都有一个公约数，那取模结果也是这个公约数的倍数。
那我们可以把中途的 \(y-x\) 都乘起来，那只要其中有一个跟 \(n\) 有非 \(1\) 共约数，你乘起来的数跟 \(n\) 就有非 \(1\) 公约数。
那我们可以多算几次 \(y-x\)，再求一次 \(\gcd\)，经过计算和测试，一般连续算 \(127\) 次再 \(\gcd\) 一次是最优的，跑的很快。

但你会发现会有一点小问题。
我们可能跑了没有 \(127\) 次，就出现了环。然后就导致还没有求出答案，就退了出来。
或者由于算法过度优秀（bushi，乘积乘着乘着直接变 \(n\) 倍数，取模直接是 \(0\)，那没有必要干等着。

那处理第一个问题，我们可以用倍增来解决，在每个 \(i=j\) 要重新换位置的时候我们也做一次 \(\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;
}