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

【模板】Pollard-Rho算法

题目链接：luogu P4718

题目大意

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

思路

Miller Rabin

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

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

费马小定理

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

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

二次探测定理

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

那转回到费马小定理：
(a^{p-1}equiv1(mod p))
(a^{{(frac{p-1}{2})}^2}equiv1(mod 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
");
			else printf("%lld
", maxn_pr); 
	}
	
	return 0;
}