Miller-Rabin素数快速检测

【update 2017-03-26】http://www.cnblogs.com/candy99/p/6624643.html

 

---

 

满足费马小定理 a^(n-1) === 1(mod n)

--->伪素数      

对于所有a belong Zn*，总存在满足的合数n，称为Carmichael数

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

【Miller-Rabin】:

1.随机找多个s个a

2.二次探测定理：　如果p是奇素数，则 x² === 1(mod p)的解为 x = 1 || x = p - 1(mod p)    ｛如：5的话，1或4｝

//Miller-Rabin
//n  prime a -->a^(n-1)===1(mod n) -->fastPowMod(a,n-1,n)==1
//warn: Carmichael/lucky

ll mulModhaoxiangmeiyonghenman(ll a,ll b,ll n){
    ll ans=0;
    for(;b;a=(a<<1)%n,b>>=1)
        if(b&1)
            ans=(ans+a)%n;
    return ans;
}

ll mulMod(ll a,ll b,ll n){    //黑科技
    ll ans=(a*b-(ll)((long double)a/n*b+0.5)*n);
    return ans<0?ans+n:ans;
}

ll powMod(ll a,ll b,ll n){
    ll ans=1;
    for(;b;a=mulMod(a,a,n),b>>=1)
        if(b&1)
            ans=(ans*a)%n;
    return ans;
}

bool witness(ll a,ll n,ll u,int t){
    ll now=powMod(a,u,n),pre=now;
    
    for(int i=1;i<=t;i++){
        now=mulMod(now,now,n);
        if(now==1&&pre!=1&&pre!=n-1)
            return true;
        pre=now;
    }
    if(now!=1) return true;
    return false;
}

bool mrP(ll n){
    if(n<=1) return false;
    if(n==2) return true;
    if((n&1)==0) return false;
    
    ll u=n-1;
    int t=0;
    while((u&1)==0) u>>=1,t++;  //n-1=2^t *u
    
    int a[6]={2,3,5,7,11,13};   //or random
    for(int i=0;i<6;i++){
        if(n==a[i]) return true;
        else if(witness(a[i],n,u,t)) return false;
    }
    return true;
}