Pollard Rho 因数分解
Pollard Rho 是一种较高效的、基于随机化的算法。它可以在期望 \(O(n^\frac{1}{4})\) 的时间内，对一个整数 \(n\) 分解出它的一个非平凡因子。
很可惜，这还是一个伪多项式算法。
试除法
对于一个正整数 \(n\)，它的因子分布是对称的——对于 \(n\) 一对因子 \((x,y)\)，满足 \(x<y,xy=n\)，则 \(x\) 必然落在 \([1,\lfloor\sqrt n\rfloor]\) 中，而 \(y\) 则必然落在 \([\lceil\sqrt n\rceil,n]\) 中。
因此，我们可以暴力扫描 \([1,\lfloor\sqrt n\rfloor]\) 范围内的数，并进行“试除”。
这个算法的复杂度是 \(O(n^\frac 1 2)\) 的。和我们的目标差别也不太大，就是次幂还需要折半而已。
生日“悖论”
这就是大名鼎鼎的生日悖论！
生日悖论有很多种描述，这里给出的描述是：
在一个房间里，有 \(k\) 个人，不妨假设他们每个人的生日都是在一年中独立均匀随机的。如果一年有 \(n\) 天，那么当房间里存在两个人的生日相同的概率不小于 \(p\) 时，\(k\) 至少是多少？
问题本身并不难解决。我们设事件 \(A\) 为“存在两个人生日相同”，那么 \(P(A)\ge p\)。反过来，\(\overline{A}\) 表示“所有人的生日互不相同”，这样 \(P(\overline A)\le 1-p\)。
容易计算出：
而根据一个稀奇古怪的不等式：
我们可以得到：
解这个不等式可以得到：
近似地得到，较接近下界的 \(k\) 为 \(\sqrt{2n\ln\frac{1}{1-p}}\)，因此当 \(1-p\) 不太小的时候，我们都可以认为 \(k\) 是 \(O(n^{\frac 1 2})\) 的。
实际上，当 \(p=99\%\) 时，也可以算出 \(\sqrt{2n\ln \frac 1{1-p}}\approx \sqrt{8n}\)。
Pollard Rho 因数分解
基础思想
假如我们得到了一个大整数 \(n\)，需要找到 \(n\) 的一个非平凡因子。我们可以用随机化的方法来“优化”试除法——随便在 \([2,\lfloor\sqrt n\rfloor]\) 范围内选一个数，看它是不是 \(n\) 的一个因子，如果是，我们就终于成功了。
这样一来，选中一个因子的概率大约是 \(\frac{d(n)}{2\sqrt n}\)，因此该算法的期望运行次数为 \(O(\frac{\sqrt n}{d(n)})\)。考虑到 \(d(n)\) 增长比较缓慢，这样优化的效果不是很明显。
但是，我们并没有必要正正好选出一个因子。实际上，如果我们可以在 \([2,n-1]\) 中找到一个合适的 \(k\)，使得 \(\gcd(k,n)>1\)，我们实际上也成功找到了一个 \(n\) 的非平凡因子。虽然因子比较少，但是满足 \(\gcd(k,n)>1\) 的 \(k\) 会多不少，因此成功的概率也会高不少。
随机数
假如我们有一个随机整数的序列，为 \(a_1,a_2,\dots,a_m\)。再假设 \(n\) 有一个非平凡的因子 \(p\)。
我们考察这样两个序列：\(b_i=a_i\bmod n,c_i=a_i\bmod p\)。如果我们发现了 \(c_i=0\) 且 \(b_i\not=0\)，那么我们直接取 \(\gcd(n,c_i)\) 就找到了一个 \(n\) 的非平凡因子。
但是，这样一个暴力的做法，就好像上述的“随机因子”的方法，太直接了。我们可以套用生日悖论：如果存在 \(i,j\)，满足 \(c_i=c_j\) 且 \(b_i\not=b_j\)，我们也算找到了一个非平凡因子 \(\gcd(n,|c_i-c_j|)\)。由于 \(c\) 本质上也是一个随机序列，因此我们知道 \(c\) 中不同的数的数量是 \(O(\sqrt p)\) 的。如果我们认为 \(p\) 是 \(n\) 的最小非平凡因子，则应有 \(p\le \sqrt n\)。因此，该算法的期望时间复杂度（也就是 \(m\)）是 \(O(n^{\frac 1 4})\) 的。
我们假设基础运算、生成随机数、取 \(\gcd\) 都是可以 \(O(1)\) 完成的。
虽然实际实现时基本上不可能做到。
“随机数”生成和 Rho 环
一般来说，在 Pollard Rho 的实现中，我们一般会采取 \(x_{i+1}=x_i^2+c\) 来迭代生成下一个“随机数”，其中 \(c\) 是一个随机的常数。
不要问为什么要用这个公式迭代，问就是可能和分形有关，利用它本身的复杂性吧
但是，所谓“随机”也只是在一定范围内才能当作随机的“伪随机”罢了。oi-wiki上面举了 \(x_1=1,c=2,n=50\) 的例子，这个情况下 \(\{x_n\}\) 迅速地掉入了一个循环。
顺便一说，正因为循环的存在，通过排列之后我们可以得到一个类似于 \(\rho\) 的形状，这正是 Pollard Rho 中 Rho 的来历吧。
虽然我觉得所有不幸的伪随机数序列最后都会掉到这样一个环里
处理环
如果我们发现了 \(b\) 中出现了循环，那么之后的所有搜索都是无效的，因为我们是在做重复的无用功。这个时候我们应当及时地退出，重新选定种子再搜索因子。
问题就出现了：如何发现我们有没有进入一个循环？Trivial 的方法是，直接开一个数组用于标记，但在 \(n\) 相当大的时候这个方法就不管用了。不过，有两种巧妙的方法可以帮助我们解决这个问题：
Floyd 判环法
这是应当被收藏的巧妙方法之一。
我们假想，有一个人和他的宠物犬在环形操场上跑步。假设人的速度为 \(v_h\)，狗狗的速度为 \(v_d\)，并且满足 \(v_d=2v_h\)。那么，在出发后，人和犬再相遇时，我们就知道此时两者的路程之差一定是环长的正整数倍。
放到我们的问题上来，我们可以利用类似的技巧。维护指针 \(i,j\)，并始终让 \(j=2i\)。这样，当我们发现 \(x_i=x_j\) 时，就说明我们一定走到环上了。
倍增判环法
由于环长是有限的，因此我们大可对环长进行倍增。
具体地说，我们可以依次走 \(1,2,4,8,\dots,2^k,\dots\) 步。维护指针 \(i\)，指向起点；维护另一个指针 \(j\)，指向当前走到的位置。如果我们发现 \(x_i=x_j\)，则一定出现环了，退出整个流程；否则直到 \(j=i+2^k\) 时这个阶段才会终止。移动起点（否则 \(i\) 很有可能从来就没有走到过环上），并进行下一阶段的移动。
算法实现
样本累积
频繁调用计算 \(\gcd\) 的函数会导致运行效率较低下。利用 \(\gcd\) 的性质：如果 \(\gcd(n,a)>1\)，那么对于 \(c\in N_+,\gcd(ac\bmod n,n)=\gcd(ac,n)>1\)，我们可以将样本累积多次，一次求出 \(\gcd\) 且不会导致我们找不到因子。
一般来说，样本累积数大概会取 \(2^k-1\)，也就是每累积 \(2^k-1\) 份就做一次 \(\gcd\)。实际实现的时候，通常会取 \(k=7\)。
原因我也不太清楚，可能就是在平衡复杂度吧
参考实现
template<typename _T>
_T ABS( const _T a ) {
return a < 0 ? -a : a;
}
unsigned GetSeed() { char *tmp = new char; return time( 0 ) * ( unsigned long long ) tmp; }
LL PollardRho( const LL n ) {
static std :: mt19937 rng( GetSeed() );
std :: uniform_int_distribution<LL> gen( 1, n - 1 );
bool cir; int acc;
LL x, y, c, d, smp;
while( true ) {
cir = false;
acc = 0, smp = 1;
y = x = gen( rng ), c = gen( rng );
for( int stp = 1 ; ; stp <<= 1, x = y ) {
for( int i = 0 ; i < stp ; i ++ ) {
y = ( ( __int128 ) y * y % n + c ) % n;
smp = ( __int128 ) smp * ABS( x - y ) % n, acc ++;
if( x == y || smp == 0 ) { cir = true; break; }
if( acc == 127 ) {
d = Gcd( n, smp );
if( d > 1 ) return d;
acc = 0, smp = 1;
}
}
if( cir ) break;
d = Gcd( n, smp );
if( d > 1 ) return d;
}
}
}