【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;
}

__EOF__

本文作者：あおいSakura
本文链接：https://www.cnblogs.com/Sakura-TJH/p/luogu_P4718.html
关于博主：评论和私信会在第一时间回复。或者直接私信我。
版权声明：本博客所有文章除特别声明外，均采用 BY-NC-SA 许可协议。转载请注明出处！
声援博主：如果您觉得文章对您有帮助，可以点击文章右下角【推荐】一下。您的鼓励是博主的最大动力！