【学习笔记】Pollard-Rho算法

前言

第一次打这个题还是 $2019/9$ 的事，当时一直没过。现在才算搞懂一些了。

日子过的真快啊。正贴合了《红楼梦》里的话：真把光阴虚度，岁月空添！

学习资料

这篇题解绝对是最好的一份！看了他的文章，我的文章也就是只剩下一些细节上的补充了。

其实，抛下 $Ak$ 的心思，学习本身是充满乐趣的。题解时常会有真知灼见，有真正的通晓者，有那些见多识广的前辈，看他们的文章就让人恍然大悟。有的文章就比较水……

算法解析

算法的目的是，快速找到 $n$ 的一个因数。——比试除法 $\mathcal{O}(\sqrt{n})$ 快得多！

碰运气

随机一个数，看看它是不是 $n$ 的因数。这样的正确性极低。

考虑利用 $\gcd$，随机到 $n$ 的任一因数的倍数即可。概率确实提升了，但是也没提升多少……

随机序列

考虑一个迭代生成的序列 $x_{i+1}=(x_i^2+1)\mathrm{m}\mathrm{o}\mathrm{d}n$，由于取值有限，必然会出现重复。而一旦出现 $x_i=x_j(i<j)$，它们 之后的结果是完全相同的，也就是出现了一个 循环节。

假设 这个 $x_i$ 是等概率随机的，根据我们熟知的 生日悖论，循环节与循环节之前的部分，都应该是期望长度 $\mathcal{O}(\sqrt{n})$ 的。这个在 $\text{OI Wiki}$ 上有讲解：用 $1+x⩽e^x$ 进行放缩（因为 $1+x$ 是 $e^x$ 的一条切线）。第 $k$ 次不重复的概率是 $1-\frac{k-1}{n}⩽\exp \left[\frac{-(k-1)}{n}\right]$，那么能够坚持 $k$ 步的概率就不超过 $\exp \left[\frac{-k(k-1)}{2n}\right]$ 。对它积分就可以得到期望。$\mathrm{w}\mathrm{o}\mathrm{l}\mathrm{f}\mathrm{r}\mathrm{a}\mathrm{m}\mathrm{a}\mathrm{l}\mathrm{p}\mathrm{h}\mathrm{a}$ 告诉我，它是 $\sqrt{\frac{n\pi }{2}}\sqrt[8n]{e}\mathrm{erf}(\frac{2k-1}{2\sqrt{2n}})$，所以就是 $\sqrt{n}$ 了。当然，打表也发现期望是 $\frac{5}{4}\sqrt{n}$ 左右，与结论相符。

然后，考虑 $n$ 的一个质因数 $p$ 。如果把递推公式的模数换成 $p$，那么得到的序列 $x_i^{\prime }=x_i\mathrm{m}\mathrm{o}\mathrm{d}p$，就是在原有的 $x_i$ 基础上再对 $p$ 取模。这是显而易见的。所以期望仍然是 $\mathcal{O}(\sqrt{p})$ 步走进一个循环节。模 $p$ 意义下的环，即 $x_a\equiv x_b(\mathrm{m}\mathrm{o}\mathrm{d}p)$，只要找到它们，就会发现 $(x_a-x_b)$ 是 $p$ 的倍数，然后就会找到这个质因数了！

显然 $n$ 的最小质因数 $p⩽\sqrt{n}$，然后环长是 期望 $\mathcal{O}(\sqrt{p})$ 的，所以复杂度是 期望 $\mathcal{O}(n^{1/4})$ 的。然而这只是期望啊！

注意：不能使用梅森旋转算法生成序列，因为它不是多项式型变换，两个模 $p$ 意义下相同的值 $x_1,x_2$ 的下一个值是可能不同的！所以此时数值会重复，但是并不会构成严格意义上的环，完整的 “环” 长度期望仍然是 $\mathcal{O}(\sqrt{n})$ 而非 $\mathcal{O}(\sqrt{p})$ 的！

$\text{Floyd}$ 判环

怎么找环？只需要让其中一个以 $1$ 单位速度前进，另一个以 $2$ 单位速度前进，二者相遇就是环了。巧妙的是，它同样可以（尽力）保证 $x_a\ne x_b$ 。如果 $x_a=x_b$，那就说明 $\mathrm{m}\mathrm{o}\mathrm{d}$ $n$ 和 $\mathrm{m}\mathrm{o}\mathrm{d}$ $p$ 意义下的环是等长的！可是那样就永远找不到了，怎么办？

解决方案是，再做一次，调整参数。当然，这种情况是比较少的，毕竟 $\sqrt{n}$ 和 $\sqrt{p}$ 差挺多的……

积累求 $\gcd$

求 $\gcd$ 是很慢的操作，可以将很多个乘在一起，然后再求一次 $\gcd$ 。

如果累积的太多，就不能立刻得到结果。比如你积累$\sout{10^6}$个样本只检测一次。

样本数量为 $100\sim 128$ 个时，效率比较高。其实这个 $128$ 就约略是 $2\log n$，也就是 $\gcd$ 的复杂度；这里相当于一个小的分块，做了一个复杂度的平衡。

代码实现

奇妙情况

有一种很奇妙的情况：积累的样本的乘积是 $n$ 的倍数。不妨设它是 $a\times b$ 导致的。由于 $b\in [1,n)$，并且这是第一次出现（即 $1⩽a<n$ 成立），只有可能是 $a,b$ 都与 $n$ 不互质。

那么直接 $\mathtt{b}\mathtt{r}\mathtt{e}\mathtt{a}\mathtt{k}$ 让 $a$ 去得到答案就行了。

示例代码

#include <cstdio> // XJX yyds!!!
#include <iostream>
#include <cstring>
#include <algorithm>
#include <cctype>
#include <vector>
#include <random>
using namespace std;
# define rep(i,a,b) for(int i=(a); i<=(b); ++i)
# define drep(i,a,b) for(int i=(a); i>=(b); --i)
typedef long long llong;
inline int readint(){
	int a = 0, c = getchar(), f = 1;
	for(; !isdigit(c); c=getchar())
		if(c == '-') f = -f;
	for(; isdigit(c); c=getchar())
		a = (a<<3)+(a<<1)+(c^48);
	return a*f;
}
void writeint(int x){
	if(x > 9) writeint(x/10);
	putchar(char((x%10)^48));
}

inline llong modMul(const llong &a,const llong &b,const llong &MOD){
	llong c = a*b-llong((long double)a/MOD*b+1e-8)*MOD;
	return ((c %= MOD) < 0) ? (c+MOD) : c;
}
inline llong qkpow(llong b,llong q,llong MOD){
	llong a = 1; b %= MOD;
	for(; q; q>>=1,b=modMul(b,b,MOD))
		if(q&1) a = modMul(a,b,MOD);
	return a;
}
bool is_prime(const llong &p){
	const static int BASE[] = {2,325,9375,28178,450775,9780504,1795265022};
	llong phi = p-1; int8_t x = 0;
	for(; !(phi&1); phi>>=1,++x);
	for(int i=0; i!=7; ++i){
		if(BASE[i]%p == 0) continue;
		llong v = qkpow(BASE[i],phi,p);
		if(v == 1) continue; // just fine
		for(int8_t j=0; j!=x; ++j){
			llong nxt = modMul(v,v,p);
			if(nxt == 1){
				if(v != p-1) return false;
				v = 1; break; // become 1
			}
			v = nxt; // just square
		}
		if(v != 1) return false;
	}
	return true;
}
inline llong getGcd(llong a,llong b){
	while(b && (a%b)) a %= b, b %= a;
	return b ? b : a;
}

std::mt19937_64 rnd;
llong ans;
void pollard(llong n){
	if(n <= ans) return ; // cut branch
	if((n&(n-1)) == 0) return ;
	if(is_prime(n) == true)
		return void(ans = n);
	__bad_luck:
	llong c = rnd()%(n-3)+3, d = 1;
	llong x = rnd()&INT_MAX, y = (x*x+c)%n;
	for(const int w=64; d==1&&x!=y; ){
		llong now = 1, tmp;
		for(int i=0; i!=w&&x!=y; ++i){
			tmp = modMul(now,y-x,n);
			if(tmp != 0) now = tmp;
			else break; // we found it!
			x = (modMul(x,x,n)+c)%n;
			y = (modMul(y,y,n)+c)%n;
			y = (modMul(y,y,n)+c)%n;
		}
		d = getGcd(n,now);
	}
	if(d == 1) goto __bad_luck;
	while(n%d == 0) n /= d;
	return pollard(n), pollard(d);
}

int main(){
	llong n;
	for(int T=readint(); T; --T){
		scanf("%lld",&n);
		if(is_prime(n)) puts("Prime");
		else ans = ((n&1)^1)<<1, pollard(n), printf("%lld\n",ans);
	}
	return 0;
}