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;
}