@1 - 素性测试：Miller-Rabin算法@
@1.1 - 算法来源@
假如我们需要检测一个数 x 是否为素数，我们应该怎么做？
最显然的就是从 2~n-1 尝试去找到 x 的另一个因数。
当然可以稍微优化一点：只需从 2~√x 中寻找，这个算法复杂度是 O(√x) 的。
那么是否有更优秀的算法？在密码学中所使用的数论，如果仅用 O(√x) 的算法是远远不足的。
更优的确定性算法的确存在，但过于复杂，在算法竞赛中不大适用。
我们不妨转换方向，设计一个更为简便的随机算法。
考虑费马小定理：
它的逆定理为：
这个定理不一定成立，但是很大概率上成立。
我们不妨选择一些 a，检测定理是否成立，以此判断 p 是否为素数。
于是这个算法正确的概率非常高。
看起来很有道理，但是——存在一些合数 x，对于所有的 a 都满足 (a^{x-1} = 1 mod x)。
我们需要对这个算法进一步改进，而这个改进需要用到下面的二次探测定理：
它的逆定理为：
这个逆定理一样很大概率上成立。
可以证明，对于每一个合数 x，必然存在一个 a 不会同时满足上面两个逆定理。
@1.2 - 算法描述@
假如我们要检测的数为 x。
我们选取一些质数 a1, a2, ..., ak。
先排除掉 x = ai 以及 x 是 ai 的倍数的情况，以加快速度。
令 x - 1 = 2^k * m，且 m 是奇数（即提取 p - 1 中的 2 的幂）。
因为 x 为偶数的情况上面已经排掉了，所以 k >= 1。
接下来，对于每一个 ai 做一遍下面的过程：
求出 n0 = ai^m，如果此时 n0 = 1 或 p - 1 就视为检测成功。
依次求出 n1 = 2*n0，n2 = 2*n1 = 2^2*n0，...，nk = 2^k*n0。
依次检验 n1~nk。对于 ni，如果 ni = p - 1 视为检测成功；否则如果 ni = 1 视为检测失败；如果最终 nk ≠ 1 则视为检测失败。
可以发现上面的检测既要求满足二次探测，又要求满足费马小定理才会检测成功。
质数 a1, a2, ..., ak 视情况而定。
如果 x <= 10^12，a 取 {2, 13, 23, 1662803} 即可。
如果 x <= 10^18，a 取 {2, 3, 5, 7, 11, 13, 17, 19, 23} 即可。
取得太多虽然可以进一步保证正确性，但太慢。
这个算法，检测失败一定失败，检测成功不一定成功。这就是它的随机性所在。
但在算法竞赛的数据范围中是可以保证正确性的。
@1.3 - 算法实现@
算法实现还是比较巧妙的。
typedef long long ll;
const int prm[] = {2, 3, 5, 7, 11, 13, 17, 19, 23, 29};
ll mul_mod(ll a, ll b, ll mod) {
ll ret = 0;
while( b ) {
if( b & 1 ) ret = (ret + a)%mod;
a = (a + a)%mod;
b >>= 1;
}
return ret;
}
ll pow_mod(ll b, ll p, ll mod) {
ll ret = 1;
while( p ) {
if( p & 1 ) ret = mul_mod(ret, b, mod);
b = mul_mod(b, b, mod);
p >>= 1;
}
return ret;
}
bool Miller_Rabin(ll x, ll y, ll z) {
ll pw = pow_mod(y, z, x);
if( pw == 1 || pw == x-1 ) return true;
for(;z!=x-1;z<<=1) {
pw = mul_mod(pw, pw, x);
if( pw == x-1 ) return true;
else if( pw == 1 ) return false;
}
return false;
}
bool Prime_Test(ll x) {
for(int i=0;i<10;i++)
if( x == prm[i] ) return true;
else if( x % prm[i] == 0 ) return false;
ll k = x-1;
while( k % 2 == 0 ) k /= 2;
for(int i=0;i<10;i++)
if( !Miller_Rabin(x, prm[i], k) ) return false;
return true;
}
@2 - 因数分解：Pollard-Rho算法@
@2.0 - 参考资料@
@2.1 - 算法来源@
同素性测试一样，pollard-rho 算法是结合代码简便+时间相对优秀为一体的算法，方便用于算法竞赛。
假如 N 为合数（此处需要先用素性测试测试 N 是否为合数）。
则存在 p <= q 且 p ≠ 1，满足 N = p*q。
我们生成一个随机数列 {xi}，使得 x1 = rand(), xi = f(xi-1) mod N。
如果存在 i, j 满足 xi = xj mod p 且 xi ≠ xj mod N，则我们就算是找到了这个 p。
而取 d = gcd(xi - xj, N) 就一定可以找到 N 的一个非平凡因子 d。
注意按照我们定义出来的随机数列 {xi}，它必然存在一个循环节。且根据生日悖论，该循环节 r1 的期望长度为 O(√N)。
令数列 {yi} 为 yi = xi mod p，则 {yi} 的循环节 r2 的期望长度为 O(√p)，且是 r1 的一个因子。
当 r1 ≠ r2 时，我们就可以按照上面的方法找到非平凡因子 d。
@2.2 - 算法描述@
一般来说，我们可以取 f(x) = x^2 + a，其中 a 是一个随机参数。
我们取 x[i] 与 x[2*i]（具体可以令 z[i] = x[2*i]，则 z[i+1] = f(f(z[i]))）。
如果 gcd(x[i] - x[2*i], N) = 1，则 2*i - i = i 并不是个循环节，继续迭代。
如果 gcd(x[i] - x[2*i], N) = N，则 x[i] = x[2*i]，2*i - i = i 是 {xi} 的循环节；此时停止迭代，变更随机参数，再来一遍。
否则，gcd(x[i] - x[2*i], N) 就是我们要找的那个因子。
我们期望枚举 O(√p) 次以得到我们的答案，而 p = O(√N)。
所以最终算法时间复杂度期望为 (O(N^{frac{1}{4}}*alpha(N)))，其中 (alpha(N)) 为 gcd 的复杂度。
这个算法也是一个随机算法，不过和上面的不同，这个算法一定给出正确答案，但时间复杂度是期望的。
@2.3 - 算法实现@
下面的代码展示的是求一个数的最小因数的过程。
实现上依然有一些值得注意的小细节。
typedef long long ll;
const int INF = (1<<30);
const int prm[] = {2, 3, 5, 7, 11, 13, 17, 19, 23, 29};
ll mul_mod(ll a, ll b, ll mod) {
ll ret = 0;
while( b ) {
if( b & 1 ) ret = (ret + a)%mod;
a = (a + a)%mod;
b >>= 1;
}
return ret;
}
ll pow_mod(ll b, ll p, ll mod) {
ll ret = 1;
while( p ) {
if( p & 1 ) ret = mul_mod(ret, b, mod);
b = mul_mod(b, b, mod);
p >>= 1;
}
return ret;
}
bool Miller_Rabin(ll x, ll y, ll z) {
ll pw = pow_mod(y, z, x);
if( pw == 1 || pw == x-1 ) return true;
for(;z!=x-1;z<<=1) {
pw = mul_mod(pw, pw, x);
if( pw == x-1 ) return true;
else if( pw == 1 ) return false;
}
return false;
}
bool Prime_Test(ll x) {
for(int i=0;i<10;i++)
if( x == prm[i] ) return true;
else if( x % prm[i] == 0 ) return false;
ll k = x-1;
while( k % 2 == 0 ) k /= 2;
for(int i=0;i<10;i++)
if( !Miller_Rabin(x, prm[i], k) ) return false;
return true;
}
ll gcd(ll x, ll y) {
return y == 0 ? x : gcd(y, x%y);
}
ll abs(ll x) {
return x < 0 ? -x : x;
}
ll min(ll x, ll y) {
return x < y ? x : y;
}
ll Pollard_Rho(ll x) {
for(int i=0;i<10;i++)
if( x % prm[i] == 0 ) return prm[i];
ll a = rand() % (x-2) + 2, b = a;
ll c = rand() % (x-1) + 1, d = 1;
while( d == 1 ) {
a = (mul_mod(a, a, x) + c)%x;
b = (mul_mod(b, b, x) + c)%x;
b = (mul_mod(b, b, x) + c)%x;
d = gcd(abs(a-b), x);
}
return d;
}
ll Find_Min_Div(ll x) {
if( x == 1 ) return INF;
if( Prime_Test(x) ) return x;
ll y = Pollard_Rho(x);
return min(Find_Min_Div(y), Find_Min_Div(x/y));
}