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

谁借莪１个温暖的怀抱￠ 2023-06-01 12:58 5阅读 0赞

*  @1 - 素性测试：Miller-Rabin算法@
    
     *  @1.1 - 算法来源@
     *  @1.2 - 算法描述@
     *  @1.3 - 算法实现@
 *  @2 - 因数分解：Pollard-Rho算法@
    
     *  @2.0 - 参考资料@
     *  @2.1 - 算法来源@
     *  @2.2 - 算法描述@
     *  @2.3 - 算法实现@

--------------------

## @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 - 参考资料@ ###

[本节内容参考了这一篇博客][Link 1]。

### @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)); }

转载于:https://www.cnblogs.com/Tiw-Air-OAO/p/11526893.html

[Link 1]: https://blog.csdn.net/qq_39972971/article/details/82346390