@总结 - 10@ Miller-Rabin素性测试与Pollard-Rho因数分解

@1 - 素性测试：Miller-Rabin算法@
@2 - 因数分解：Pollard-Rho算法@

@1 - 素性测试：Miller-Rabin算法@

@1.1 - 算法来源@

假如我们需要检测一个数 x 是否为素数，我们应该怎么做？

最显然的就是从 2~n-1 尝试去找到 x 的另一个因数。

当然可以稍微优化一点：只需从 2~√x 中寻找，这个算法复杂度是 O(√x) 的。

那么是否有更优秀的算法？在密码学中所使用的数论，如果仅用 O(√x) 的算法是远远不足的。

更优的确定性算法的确存在，但过于复杂，在算法竞赛中不大适用。

我们不妨转换方向，设计一个更为简便的随机算法。

考虑费马小定理：

\[当 p 为素数时，a^{p-1} = 1\mod p
\]

它的逆定理为：

\[若 a^{p-1} = 1 \mod p，则 p 为素数
\]

这个定理不一定成立，但是很大概率上成立。

我们不妨选择一些 a，检测定理是否成立，以此判断 p 是否为素数。

于是这个算法正确的概率非常高。

看起来很有道理，但是——存在一些合数 x，对于所有的 a 都满足 \(a^{x-1} = 1 \mod x\)。

我们需要对这个算法进一步改进，而这个改进需要用到下面的二次探测定理：

\[当 p 为素数时，若 x^2 = 1 \mod p，则 x = 1 或 p-1 \mod p
\]

它的逆定理为：

\[若 x^2 = 1 \mod p，且 x = 1 或 p-1 \mod p，则 p 为素数
\]

这个逆定理一样很大概率上成立。

可以证明，对于每一个合数 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));

}

巴特西