题解 P4718 【【模板】Pollard-Rho算法】
Piwry
2020-07-23 20:07:37
## 前言
半年前我曾写了一次 PR 的代码，不过当时是照着 AC 代码玄学调试（~~还是概率~~）过的，注释打的也大部分都是调试坑点，和算法正确性丝毫没有关系...（~~毕竟是道黑题怎么也得想办法过掉~~）
半年后~~为了做 lxl 的某道[题](https://www.luogu.com.cn/problem/P5071 "lxl毒瘤！（划）")~~，我再次打开了这道题目。不过发现自己以前写的代码太丑了，打算重写一份。
于是和着算法导论上的代码，差不多地打了一份；交上去 T 的很惨。接着开始试着寻找关于 Pollard-Rho 算法的资料，根据对复杂度的理解再打了一份，总算是稳过了。（当然中途也少不了被数据各种卡的调试部分）
（话说过程中因看不懂算导在 ~~bb~~ 写什么就去网上搜了一圈，结果最后还是看算导理解的qaq）
写完又翻了下题解，发现几乎没人涉及到算法复杂度的部分，尤其是**期望复杂度 $\sqrt p$ 是怎么来的**。于是就打算自己另写一篇题解了qwq
以上全是废话，由于我的问题语文读着还不太通，只是纪念性地写在这。看到这句话请无视以上所有文字接着看正文（
## 算法解析
首先我们思考一个下一项仅由上一项决定，且值域为模 $n$ 的剩余系的序列 $\text{<}x\text{>}$
例如 $x_i=(x_{i-1}^2+1)\pmod n$
由于它某项的值仅由前一项决定，且每一项可能的取值是有限的，因此该数列一定会在经历一定次数的迭代后陷入循环
更具体的来说，只要该序列中**出现一对相同的数**，就会在这对数间产生环
联想到生日悖论，若这个数列是期望随机的，则它的 "环" 及 "尾" 的期望长度就是 $\sqrt n$ 的
为了说明方便，下面我们就先假设由该递推式导出的序列都是**足够随机**的，并以此继续讨论。
现在我们假设 $n$ 不是质数，或者说它至少有一对非平凡因子（指该因子不为 $1$ 或 $n$）；则 $n$ 可以表示为 $n=p_1^{r_1}\cdot p_2^{r_2}\cdots p_m^{r_m}$
这里我们随意取一个 $p_i$ 并设为 $p$
再考虑一个序列 $\text{<}x'\text{>}$，其中 $x'_i=x_i\pmod p$
我们可以发现，由于数列 $\text{<}x\text{>}$ 的递推式仅有加减乘除运算，因此数列 $\text{<}x'\text{>}$ 其实和它有着同样的递推式：
${x'}_i=x_{i-1}\pmod p$
${x'}_i=((x_{i-1}^2+1)\pmod n)\pmod p$
${x'}_i=(x_{i-1}^2+1)\pmod p$
${x'}_i=((x_{i-1}\pmod p)^2+1)\pmod p$
${x'}_i=({x'}_{i-1}^2+1)\pmod p$
由此可以预知序列 $\text{<}x'\text{>}$ 也会产生 "环"。而根据生日悖论，序列 $\text{<}x'\text{>}$ 的环的长度就是期望 $\sqrt p$ 的。
注意到当 ${x'}_i\equiv {x'}_j\pmod p$ 时，一定有 $p | (|x_i-x_j|)$；此时若$x_i\neq x_j$，通过计算 $\gcd(|x_i-x_j|, n)$，一定就可以得到 $n$ 的因子 $p$
若我们沿着 "$\rho$" 的路径遍历数列，不难想到，是可以在期望记录 $n^{\frac 1 4}$ 个元素后就能找出一对 $x_i, x_j$，并以此计算出 $n$ 的某个非平凡因子
## 算法实现
### 随机数
首先我们要考虑就是序列 $\text{<}x\text{>}$ 的数是否足够随机
这里我没找到可行的证明（qaq）；不过根据算法实际效率来看，这种生成方式是足够优秀的
不过为了保证数列在更多情况下依然足够随机，我们使用更一般的公式 $x_i=x_{i-1}^2+c$，并在初始时随机带入一个 $x$ 和 $c$ 并进行迭代
不过注意当 $c=0$ 或 $c=2$ 的情况应当避免（来自算法导论，书上也没有具体解释原因）
### 初始的实现方式
设最大的不在 "环" 上的序列编号 $s$，对于每一个序号 $t, t>s$，一定能找到一个序号 $t'$ 使得 $x_t\equiv x_{t'}\pmod p$
因此我们可以考虑设一个定点 $y$ 和一个动点 $x$；$x$ 根据序列的迭代公式一步步计算，期望会在 $n^{\frac 1 4}$ 步内找到 $n$ 的一个非平凡因子
但我们很难确定 $y$ 是否在序列的 "环" 上。具体实现时我们倍增地设置 $y$ 的序号即可
```cpp
/*其中 rand(l, r) 代表取这个闭区间的随机数*/
ll pr(ll n){
ll x =rand(0, n-1), y =x;
ll c =rand(3, n-1), d =1;
for(int i =1, k =2; d == 1; ++i){
x =(x*x+c)%n;
d =gcd(abs(y-x), n);
if(i == k){
y =x;
k <<=1;
}
}
return d;
}
```
### Floyd 判环
上面的算法在大多数情况下已经足以应付了
但我们发现，在极少数情况下，该算法甚至在遍历整个模 $n$ 的 "$\rho$" 数列后仍然无法找到解
这时候我们就要判环并重新尝试
具体来说，我们另设一个点 $y$ 和 $x$ 在同一起点出发，但该点以 $x$ 两倍的速度按公式迭代，当 $x\equiv y\pmod n$ 时我们就找到了环。期望环的长度是 $\sqrt n$ 的，我们可以证明该算法能在 $O(\sqrt n)$ 的时间内判断出环
至于算法的证明，设起点距离进入环还有 $s$ 步，环长 $c$，点 $x$ 走了 $k$ 步；则想要令 $x, y$ 相遇就需要使 $k-s\equiv 2k-s\pmod c$ 且 $k> s$，这显然可以做到
同时注意在算法结束时，点 $y$ 会至少遍历环一次，但 $x$ 不一定会
### 基于 Floyd 的实现
由于序列 $\text{<}x'\text{>}$ 可以看做是 $\text{<}x\text{>}$ 的较小版本，或者说这两个序列存在一一对应的关系，我们的 Floyd 判环其实对序列 $\text{<}x'\text{>}$ **也是适用**的。
具体来说，当我们发现 Floyd 算法内的两动点满足 $x\equiv y\pmod p$ 时（尽管我们可能不会检查这个式子），$\gcd(|x-y|, n)$ 就会给出 $n$ 的一个非平凡因子。而 Floyd 算法显然会比我们上面的倍增法要快许多
```cpp
ll pr(ll n){
/*因为初始跳两步的原因，该写法没法分解 4*/
if(n == 4) return 2;
ll x =rand(0, n-1), y =x;
ll c =rand(3, n-1);
ll d =1;
x =(x*x+c)%n;
y =(y*y+c)%n, y =(y*y+c)%n;
while(x != y){
x =(x*x+c)%n;
y =(y*y+c)%n, y =(y*y+c)%n;
d =gcd(abs(x-y), n);
if(d != 1)
return d;
}
return n;
}
```
### 初始实现的判环改进
当然，如果加上判环部分，一开始的初始实现也是可行的
```cpp
ll pr(ll n){
/*因为初始跳两步的原因，该写法没法分解 4*/
if(n == 4) return 2;
ll x =rand(0, n-1), y =x, y2 =x;
ll c =rand(3, n-1);
x =(x*x+c)%n;
y2 =(y2*y2+c)%n, y2 =(y2*y2+c)%n;
for(int i =1, k =2; x != y2; ++i){
x =(x*x+c)%n;
y2 =(y2*y2+c)%n, y2 =(y2*y2+c)%n;
ll d =gcd(abs(x, -y), n);
if(d != 1)
return d;
if(i == k){
y =x;
k <<=1;
}
}
return n;
}
```
### 倍增积累 gcd
我们发现在模 $n$ 意义下将每次的 `abs(x-y)` 相乘，若对于某一个 $|x_i-y_i|$ 有 $d=\gcd(|x_i-y_i|, n)$ 且 $d$ 为 $n$ 的非平凡因子；则它们的积，这里设为 $cnt$，也一定有 $d'=\gcd(|x_i-y_i|, n)$ 且 $d'$ 为 $n$ 的非平凡因子
具体来说，我们可以积累一定的样本再做 $\gcd$，这样能节省一些时间，且不会影响每个样本的答案
这个阈值可以设为倍增的，且可以设一个上限。以本题数据来说设为 $128$ 会有不错的效果
这里仅给出以 "基于 Floyd 的实现" 为模板改进的代码
```cpp
ll pr(ll n){
if(n == 4) return 2;
ll x =rand(0, n-1), y =x;
ll c =rand(3, n-1);
ll d =1;
x =(x*x+c)%n;
y =(y*y+c)%n, y =(y*y+c)%n;
for(int lim =1; x != y; lim =min(128, lim<<1)){
ll cnt =1;
for(int i =0; i < lim; ++i){
ll tmp =cnt*abs(x-y)%n;
/*这时要么原先的 cnt 以及此时的 |x-y| 含 n 质因数 ，要么 x-y == 0*/
/*为了避免之前的样本丢失，直接推出循环并对之前积累的样本做一次计算*/
if(!tmp)
break;
cnt =tmp;
x =(x*x+c)%n;
y =(y*y+c)%n, y =(y*y+c)%n;
}
d =gcd(cnt, n);
if(d != 1)
return d;
}
return n;
}
```
## CODE
代码里有详尽的注译
其中 `pr()` 函数中，我将 "初始实现的判环改进" 注译起来了，并标为 "算法2"
同时代码中使用的随机数生成函数是 `std::mt19937`；不过在本题的数据表现中，其和一般的 `rand()` 函数差别不大
```cpp
#include <cstdio>
#include <cstdlib>
#include <ctime>
#include <cmath>
#include <random>
#define ll long long
#define ull unsigned long long
#define lll __int128
//#pragma GCC target("avx")
//#pragma GCC optimize("Ofast")
//#pragma GCC optimize(2)
/*------------------------------Base------------------------------*/
inline ll mul(ll a, ll b, const ll M){
ll d =a*(long double)b/M;
ll ret =a*b-d*M;
if(ret < 0)
ret +=M;
if(ret >= M)
ret -=M;
return ret;
}
inline ll plus2(ll a, ll b, const ll M){
ll d =(a+(long double)b)/M+0.5;
ll ret =a+b-d*M;
if(ret < 0)
ret +=M;
return ret;
}
inline ll plus(const lll a, const lll b, const lll M){
lll c =a+b;
if(c >= M)
return c-M;
else
return c;
}
ll Pow(ll x, ll k, const ll M){
ll ret =1;
for(; k; x =mul(x, x, M), k >>=1) if(k&1) ret =mul(ret, x, M);
return ret;
}
ll gcd(ll a, ll b){
while(b) b ^=a ^=b ^=a %=b;
return a;
}
inline ll Abs(ll x){ return (x < 0) ? -x : x; }
/*------------------------------Rand------------------------------*/
static std::mt19937 engine;
/*------------------------------Miller Robin------------------------------*/
bool mr(ll p){
if(p < 2) return 0;
if(p == 2) return 1;
if(p == 3) return 1;
std::uniform_int_distribution<ll> Rand(2, p-2);
ll d =p-1, r =0;
while(!(d&1)) ++r, d >>=1;
for(ll k =0; k < 10; ++k){
/*[2, p-2]*/
ll a =Rand(engine);
ll x =Pow(a, d, p);
if(x == 1 || x == p-1) continue;
for(int i =0; i < r-1; ++i){
x =mul(x, x, p);
if(x == p-1) break;
}
if(x != p-1) return 0;
}
return 1;
}
/*------------------------------Pollard Rho------------------------------*/
using std::min;
/*貌似 -c 在本题的数据的效率高些 (0.2s)，原因尚不确定*/
inline ll getnext(ll x, ll c, ll n){ return plus(mul(x, x, n), -c, n); }
ll pr(ll n){
/*因为初始跳两步的原因，下面写法均没法分解 4 (即使下面的 Rand 范围设置为 [0, n-1] )*/
if(n == 4) return 2;
std::uniform_int_distribution<ll> Rand(3, n-1);
ll x =Rand(engine), y =x;
ll c =Rand(engine);
ll d =1;
/*以下两种写法的期望复杂度都是正确的，但写法 1 的表现更好*/
/*----------写法 1----------*/
x =getnext(x, c, n);
y =getnext(y, c, n), y =getnext(y, c, n);
for(int lim =1; x != y; lim =min(128, lim<<1)){/*提升约 0.1s */
// for(int lim =1; x != y; lim =lim<<1){
ll cnt =1;
for(int i =0; i < lim; ++i){
ll tmp =mul(cnt, Abs(plus(x, -y, n)), n);
if(!tmp)/*提升约 0.6s；这时要么原先的 cnt 含 n 质因数 (这时 x-y 也含 )，要么 x-y == 0*/
break;
cnt =tmp;
x =getnext(x, c, n);
y =getnext(y, c, n), y =getnext(y, c, n);
}
d =gcd(cnt, n);
if(d != 1)
return d;
}
return n;
/*----------写法 2----------*/
/*这里还加了倍增 gcd 优化，可以作为参考*/
/*
x =getnext(x, c, n);
y =getnext(y, c, n), y =getnext(y, c, n);
ll x2 =x, cnt =1;
for(int i =1, i2 =1, k =2, lim =1; x != y; ++i, ++i2){
x =getnext(x, c, n);
y =getnext(y, c, n), y =getnext(y, c, n);
// d =gcd(Abs(plus(x, -x2, n)), n);
ll tmp =mul(cnt, Abs(plus(x, -x2, n)), n);
if(tmp)
cnt =tmp;
if(i2 == lim || !tmp || x == y){
i2 =1;
lim =min(128, lim<<1);
d =gcd(cnt, n);
if(d != 1)
return d;
cnt =1;
}
// if(d != 1)
// return d;
if(i == k){
x2 =x;
k <<=1;
}
}
return n;*/
}
ll mxp;
inline void push(ll p){
if(p > mxp)
mxp =p;
}
/*函数要求保证 n 可分解*/
void dfs(ll n){
srand(time(0));
ll d =pr(n), d2;
while(d == n)
d =pr(n);
d2 =n/d;
if(mr(d))
push(d);
else
dfs(d);
if(mr(d2))
push(d2);
else
dfs(d2);
}
ll getfact(ll n){
mxp =0;
if(mr(n))
return n;
else
dfs(n);
return mxp;
}
/*------------------------------Main------------------------------*/
inline ll read(){
ll x =0; char c =getchar();
while(c < '0' || c > '9') c =getchar();
while(c >= '0' && c <= '9') x = (x<<3) + (x<<1) + (48^c), c =getchar();
return x;
}
int main(){
srand(time(0));
for(int t =0, T =read(); t < T; ++t){
ll n =read();
ll fact =getfact(n);
if(fact == n)
puts("Prime");
else
printf("%lld\n", fact);
}
}
```