Pollard的rho启发式因子分解算法 & [CodeVS 4939] 欧拉函数：Miller-Rabin + Pollard-rho 质因数分解

Pollard的rho启发式因子分解算法用于给出整数的一个因子。在一定的合理假设下，如果n有一个因子p，可在O(p√)的期望时间内可找出n的一个因子p。

关于其复杂度，Wikipedia是这样叙述的：

If the pseudo random number x = g(x) occurring in the Pollard ρ algorithm were an actual random number, it would follow that success would be achieved half the time, by the Birthday paradox in O(p^(1/2)) ≤ O (n^(1/4)) iterations. It is believed that the same analysis applies as well to the actual rho algorithm, but this is a heuristic claim, and rigorous analysis of the algorithm remains open.

以下绝大部分来自《算法导论》。

Pollard-rho算法的核心是用递推式xi+1=(x2i+c)modn生成一个最终会进入循环的“随机”序列（表示成有向图，看起来就像ρ，这就是算法名字的由来）。虽然只显式地生成了一个序列，实际上同时生成了许多形如ρ的序列（后面将会推证）；只要两个指针都进入某个ρ的圈圈里，把它们所指向的值作差，取绝对值，和n求gcd，就能得到n的一个因子。

伪代码如下：

Pollard-Rho(n, c)
i = 1
k = 2
y = x = a random integer in [0, n)
d = 1
while d == 1
i = i+1
x = (x*x+c) mod n
if x == y
return n
d = gcd(n, abs(x-y))
if i == k
k = k*2
y = x
return d

它的正确性是显然的。算法可能会失败地返回一个平凡因子n，也可能成功地返回一个n的某个非平凡因子。

设xi+1=(x2i+c)modn的循环大小为C，循环的第一项是xt，进入循环后的某一时刻，k会被赋予一个不小于C的值，此时的x被保存为y，再转一圈，y固定不动，x会回到y，算法以失败终止。

必有2的某次幂落在区间[t, 2t]内，因此在不超过第2t步，循环内的某个值将被保存。此时的k可能小于C，无妨，因为必有2的某次幂落在区间[C, 2C]内，2C步之内，有k不小于C成立。于是，至多走(2*min(t, C)+C)步，算法终止。这个算法叫Brent判圈算法（Brent’s cycle finding method），与Floyd判圈算法均为线性，但常数优于后者（至少3C次计算后继结点）。根据Wikipedia，Pollard的原始版本采用Floyd判圈算法，后来由Richard Brant用自己新发明的判圈算法加以改进。

设n有一个因子p，其实我们同时在计算x′i=ximodp，并且递推式具有相同的形式：

x′i+1=xi+1modp=(x2i+c)modnmodp=(x2i+c)modp=((ximodp)2+c)modp=(x′2i+c)modp

同样是ρ形。循环内两个数作差，是p的倍数，和n取gcd，因子p或p的某倍便被呈现出来。与上面的论证类似，进入循环后不超过t′步，循环内的某个值将被保存，下一步，p被呈现。显然t′≤t，C′≤C。

假设{xi}是随机的，则t′和C′都是O(p√)的。

生日悖论：从n个数中可重复地随机选择k个，当k≥2n‾‾‾√时，存在两数相等的概率大于1/2。

用数学期望描述这个命题：相等数对的数目不少于1。这样计算起来会简单一些。

对1≤i<j≤k定义指示器随机变量Xij=I{第i个数和第j个数相等}，则E[Xij]=1n，

E[∑i=1k∑j=i+1kXij]=∑i=1k∑j=i+1kE[Xij]=∑i=1k∑j=i+1k1n=k(k−1)2n

令上式≥1，得到一个充分条件k≥2n‾‾‾√。

把这个结论运用到我们的问题中来。t′=O(p√)，C′=O(p√)。

至此，一切都看起来很美好。如果n有某个因子p，则O(p√)左右次循环后，我们就能找到它。

然而，推导这一切的假设显然不成立。{xi}不是随机的。另外，存在这样一种可能：所有ρ的尺寸相同，这将导致只能找到平凡因子n；这种情况需要换一个参数c，《算法导论》告诉我们，除了0和2，其他数都是不错的选择。所以，本算法的名称有定语“启发式”。

CodeVS 4939 欧拉函数

这道题是mhb同学放到CodeVS上的Orz 据其本人所言，当时试除+卡常数，最终把数据范围开到这么变态……真是太神啦……题解页面中提到的qwertyu也很神，感谢Ta的代码。

结合Miller-Rabin、Pollard-rho两个算法可以进行质因数分解。如果n是素数，返回；否则，求n的一个非平凡因子d，递归。由于本题计算的是欧拉函数值，只需要知道有哪些质因子，而不需要知道每个质因子的指数，所以可以将n中的d全部除掉再递归。质因子的数目算上重复的也不超过ceiling(log2n)，所以数组不用开很大。

#include <cstdio>
#include <cstdlib>
#include <ctime>
#include <algorithm>
typedef long long ll;
const int s = 10, MAX_F = 70;
ll cnt, f[MAX_F];

inline ll mul_mod(ll a, ll b, ll m)
{
ll c = a*b-(ll)((long double)a*b/m+0.5)*m;
return c<0 ? c+m : c;
}

ll fast_exp(ll a, ll x, ll m)
{
ll b = 1;
while (x) {
if (x & 1)
b = mul_mod(b, a, m);
a = mul_mod(a, a, m);
x >>= 1;
}
return b;
}

bool MR(ll n)
{
if (!(n&1))
return n == 2;
ll t = 0, u;
for (u = n-1; !(u&1); u >>= 1)
++t;
for (int i = 0; i < s; ++i) {
ll a = rand()%(n-2)+2, x = fast_exp(a, u, n);
for (ll j = 0, y; x != 1 && j < t; ++j, x = y) {
y = mul_mod(x, x, n);
if (y == 1 && x != n-1)
return false;
}
if (x != 1)
return false;
}
return true;
}

inline ll abs(ll x)
{
return x<0 ? -x : x;
}

ll gcd(ll a, ll b)
{
return b ? gcd(b, a%b) : a;
}

ll PR(ll n, ll a)
{
ll x = rand()%n, y = x, k = 1, i = 0, d = 1;
while (d == 1) {
if ((x = (mul_mod(x, x, n)+a)%n) == y)
return n;
d = gcd(n, abs(y-x));
if (++i == k) {
k <<= 1;
y = x;
}
}
return d;
}

void decomp(ll n)
{
if (n == 1)
return;
if (MR(n)) {
f[cnt++] = n;
return;
}
ll d = n, c = n-1;
while (d == n)
d = PR(n, c--);
do {
n /= d;
} while (!(n%d));
decomp(d);
decomp(n);
}

int main()
{
srand(time(0));
ll n;
while (scanf("%lld", &n), n) {
cnt = 0;
decomp(n);
std::sort(f, f+cnt);
cnt = std::unique(f, f+cnt)-f;
ll ans = n;
for (int i = 0; i < cnt; ++i) {
printf("%lld\n", f[i]);
ans = ans/f[i]*(f[i]-1);
}
printf("%lld\n", ans);
}
return 0;
}