BZOJ3667: Rabin-Miller算法 （Miller-Rabin&&pol_rho&&特技快速乘学习笔记）

学习了两个新算法，判素数和分解大合数。

Miller-Rabin(以下简称MR)算的时候大概就是这样子算。

---

那么怎么分解质因数呢？

这两篇博文应该已经说得比较清楚了。
http://www.cnblogs.com/thythy/p/5493624.html http://blog.csdn.net/maxichu/article/details/45459533

rho每次只找一个因数，要找n的因数，每次随机一个x出来，首先令y=x，然后令y不断做如下变化y=y*y+c(c是自己定的一个随机数)，然后用abs(x-y)去跟n求gcd，如果不为1，则找到一个，有一个步长扩二倍的优化（详见代码），不太明白为啥，估计是为了防止死循环。

找到一个因数之后递归求解（配合MR），找到所有质因数，期望复杂度$O(n^{1/4})$(口胡)

---

另外一个问题就是计算MOD为64位整数的a*b%MOD,代码给出了一种奇技淫巧，在上面链接的第一篇blog里面有讲解

//QWsin
#include<cmath>
#include<cstdio>
#include<cstring>
#include<cstdlib>
#include<iostream>
#include<algorithm>
using namespace std;
typedef long long ll;
const int a[]={2,3,5,7,11,13,17,19,23,29};

ll maxs=0;

ll mul(ll a,ll b,ll p){
    ll d=(long double)a/p*b+0.5;
    ll ret=a*b-d*p;
    if(ret<0) ret+=p;
    return ret;
}

ll Pow(ll x,ll p,ll MOD)
{
    ll ret=1;
    for(;p;p>>=1,x=mul(x,x,MOD)) if(p&1) ret=mul(ret,x,MOD);
    return ret; 
}

inline int MR(ll n)
{
    if(n==2) return 1;

    ll t=n-1;int cnt=0;
    while((t&1)==0) t>>=1,++cnt;
    for(int i=0;i<10;++i)
    {
        if(a[i]==n) return 1;
        ll fac=Pow(a[i],t,n),nxt;
        for(int i=1;i<=cnt;++i)
        {
            nxt=mul(fac,fac,n);
            if(nxt==1&&fac!=1&&fac!=n-1) return 0;
            fac=nxt;
        }
        if(fac!=1) return 0;
    }
    return 1;
}

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

ll pol_rho(ll n,ll c)
{
    ll k=2,x=rand()%n,y=x,p=1;
    for(ll i=1;p==1;++i){
        y=(mul(y,y,n)+c)%n;
        p=x>y?x-y:y-x;
        p=gcd(p,n);
        if(i==k) x=y,k<<=1;
    }
    return p;
}

inline void solve(ll n)
{
//  printf("%lld\n",n);
    if(MR(n)){maxs=max(maxs,n);return ;}

    ll t=n;
    while(t==n||t==1) 
        t=pol_rho(n,rand()%(n-1));
    solve(n/t);
    solve(t);
}

int main()
{
    srand(20010827);
    int T;cin>>T;
    while(T--) {
        ll n;cin>>n;
        maxs=0;
        solve(n);
        if(maxs==n) puts("Prime");
        else printf("%lld\n",maxs);
    }
    return 0;
}