大数的质数判定 && 质因子拆分：MR && PR
Jμdge
2018-11-13 21:56:17
最近搞这两个算法焦头烂额...但也是有一丝丝收获的。
首先声明，本博客不提供大多数证明。~~仅仅是根据结论来解释如何构造代码（emmm会用结论就好了嘛）~~
# Miller-Rabin 判定大数
大家听说过费马小定理吧？ 就是说如果一个数 p 是质数，那么就有：
>$ n^{p-1}≡1 (mod p) $
于是我们大胆猜测：如果一个数 p 对于 任意一个数 n 满足条件：
>$ n^{p-1}≡1 (mod p) $
那么这个数是质数。
...
然鹅这个猜想是错的，并且可以十分轻松的证明（打个代码试试就知道了QvQ）。
但是假如 n 是任意正整数，那么这个猜想会不会成立呢？
然鹅事实告诉我们，如此美好的性质并不存在。
比如？ 561 。
没错它通过了任意正整数 n 的检测，然鹅它并不是质数。
于是 $Miller - Rabin$ （注意是两个人）提出了另一种检测质数的方法:二次探测。
对于一个数 p，当有 X（X<p） 满足： $X^{2}≡1 (mod\ p) $ 且 $X=1 or X=p-1$ 时， p为质数。
~~证明？我不会啊。。。~~
看看某大佬的证明：
----
>∵ $X^{2}≡1\ (mod\ p) $
>∴ $X^{2}-1≡0\ (mod\ p) $
>∴ $(X+1)(X-1)≡0\ (mod\ p) $
>又 ∵ $X<p$
>∴ $X=1\ or\ X=p-1$
----
那么这个检测方法是否恒正确？
答案依旧是：很遗憾，存在反例。
比如 $46856248255981$
但是这个特别的数字是 $1e16$ 范围内唯一的反例。所以可以说 $Miller-Rabin$ 检测法的正确率已经相当高了。
~~（毕竟这个特别的数字可以特判掉）~~
那么这个算法如何实现？ 步骤如下：
>>1.将 $p-1$ 提取出所有 $2$ 的因数，假设有 $tm$ 个。设剩下的部分为 $ba$ 。
>>2.枚举(或随机)一个底数 $a$ 并计算 $r=a^{ba}(mod p)$。
>>3.将 $r$ 连续平方 $tm$ 次。如果没有出现过 $p-1$ ，那么 $p$ 未通过二次探测，不是质数。
>>4.否则，若底数已经足够，则跳出。否则回到步骤 2。
那么这样的检测正确率有保障么？ 有！ 正确率 75% （好低...）
但是我们多取几个随机数检测一下，错误的概率就渐进 $0$ 了。
（考虑每次判断错误率为 $\dfrac{1}{4} $ ，那么我们判个20 次大概就是 $\dfrac{1}{4^{20}} $ 了 ）
那么我们就可以大胆使用这个算法了（但其实这里有优化版的，不过效率不会高太多，所以还是老实写常规版的好）。
```
namespace Miller_Rabin{
inline ll qmul(ll a,ll p,ll c){ ll res=0;
for(a%=c,p%=c;p;a<<=1,a%=c,p>>=1)
if(p&1) res=(res+a)%c; return res;
} inline ll qpow(ll x,ll p,ll c){ ll s=1;
for(x%=c;p;x=qmul(x,x,c),p>>=1)
if(p&1) s=qmul(s,x,c); return s;
} inline bool check(ll x){
if(x==2) return 1;
if(x==1||(x&1^1)) return 0;
ll ba=x-1,r,ti=0,j;
while(ba&1^1) ba>>=1,++ti;
for(int i=0;i<22;++i){ //普通的话嘛，其实10来次差不多了
r=rand()%(x-1)+1;
r=qpow(r,ba,x);
if(r==1||r==x-1)
continue;
for(j=1;j<=ti;++j){
r=qmul(r,r,x);
if(r==x-1)break;
} if(j>ti)return 0;
} return 1;
}
} using namespace Miller_Rabin; ll ans;
```
```
namespace Miller_Rabin{
const ll bace[5]={2,3,7,61,24251};
//这个优化非常有趣，强制转精度的骚操作（关键没有用int128）
inline ll qmul(ll a,ll p,ll c){ static ll d;
d=(long double)a/c*p+1e-8,d=a*p-d*c; return d<0?d+c:d;
} inline ll qpow(ll x,ll p,ll c){ ll res=1;
for(x%=c;p;x=qmul(x,x,c),p>>=1)
if(p&1) res=qmul(res,x,c); return res;
} inline bool check(ll n){
if(n==1||(n&1^1)) return false;
for(int i=0;i<4;++i)
if(bace[i]==n) return 1;
static ll ba,r,j,ti; ba=n-1,ti=0;
for(int i=0;i<5;++i){ //每次判定失败的概率为 1/4 ,（配合二次探测）判个五六次其实就差不多了
ll a=rand()%(n-1)+1;
if(qpow(a,n-1,n)!=1)
return false;
} while(ba&1^1) ba>>=1,++ti;
for(int i=0;i<4;++i){ //然后是二次探测, 只用四个质数就好了
r=qpow(bace[i],ba,n);
if(r==1||r==n-1)continue;
for(j=1;j<=ti;++j){
r=r*r%n; if(r==n-1)break;
} if(j>ti)return 0;
} return 1; //经过了两步检测之后, 基本可以确定是质数
}
}
```
如果你看不懂的话...附上一些不错的文章:
>[1](https://www.xuebuyuan.com/552593.html)
>[2](http://www.cnblogs.com/TenosDoIt/p/3398112.html)
>[3](http://blog.chinaunix.net/uid-21712186-id-1818141.html)
# Pollard-Rho 分解大数
首先 PR 是建立在 MR 基础上的算法（因为分解过程中要判定质数）。
所以，确保你对 MR 算法的基本框架已经熟悉，再往下观看。
这个算法大概讲的就是：
>>1.判断当前处理的数 n 是否是质数，是的话就记录下来直接返回。
>>2.随机（又是随机）一个数 p（ <=n ），然后经过一系列的运算得到一个新的数字 q（仍不大于 n），然后判断 q 是不是与 n 互质:
>>>(1) 如果非互质的话就可以得到 n 的一个质因子，即 $gcd(n,q)$。
>>>(2) 如果互质的话就继续经过一系列的运算得到新的 q ，反复迭代
>>3.若返回值 $gcd(n,q)$ 不为 n ，则对 n 分解成功（但并不意味着 $gcd(n,q)$ 为质数），否则重复步骤 2 直至分解成功
>>4.此时我们将 $n$ 分成了两份，设一份为 $p$ ，则另一份为 $n/p$，但如步骤 3 所述， p 并非一定是质数，因此我们对这两份数分别递归分解，即将 $p$ 与 $n/p$ 置入步骤 1
然后这样我们就可以将 n 分解质因子了。
但是这里最大的疑问就是那个一系列运算究竟是怎样~~（玄学）~~ 的运算，又是如何保证正确性的。。。
相关证明看 [这里](https://blog.csdn.net/maxichu/article/details/45459533) 吧。
其次，计算 gcd 之前还有一个判断，具体与 生日悖论有关， [这篇blog](https://blog.csdn.net/doyouseeman/article/details/51204612) 里面有讲
总之会用就好了呗。
so？ 复杂度？ $O(n^{1/4})$ ，就是 n 开两次根。
附上代码
```
//namespace Miller_Rabin 这玩意儿,就是上面的
using namespace Miller_Rabin; ll ans;
ll gcd(ll x,ll y){ return y?gcd(y,x%y):x; }
const int M=(1<<7)-1;
inline ll Abs(ll x){ return x>=0?x:-x; }
inline ll g(ll x,ll n,ll a) { //对,这就是“一系列运算”
ll t=qmul(x,x,n)+a;
return t<n?t:t-n;
} ll Pollard_Rho(ll n){
if(!(n%2)) return 2;
if(!(n%3)) return 3;
ll x=0,y=1,t=1,q=1;
ll c=rand()%(n-1)+1;
for(int k=2;;k<<=1,y=x,q=1){
for(int i=1;i<=k;++i){
x=g(x,n,c),
q=qmul(q,Abs(x-y),n);
if(!(i&M)){ //这里用了生日悖论
t=gcd(q,n);
if(t>1) return t;
}
} if(t>1||(t=gcd(q,n))>1) return t;
}
} void find(ll n){
if(n==1||ans>=n)
return ;
if(check(n)){
ans=n;
return ;
} ll p=n;
while(p>=n)
p=Pollard_Rho(n);
while(n%p==0) n/=p;
find(p),find(n);
}
```
没了？没了。 ヾBye~Bye~