Miller-Rabin隨機性素數測試算法(Miller_Rabin模板)

普通的素數測試咱們有O(√ n)的試除算法。事實上，咱們有O(slog³n)的算法。

定理一：假如p是質數，且(a,p)=1，那麼a^(p-1)≡1(mod p)。即假如p是質數，且a,p互質，那麼a的(p-1)次方除以p的餘數恆等於1。（費馬小定理）

該定理的逆命題是不必定成立的，可是使人可喜的是大多數狀況是成立的。

因而咱們就獲得了一個定理的直接應用，對於待驗證的數p，咱們不斷取a∈［1，p-1]且a∈Z，驗證a^(p-1) mod p是否等於1，不是則p果斷不是素數，共取s次。其中a^(p-1) mod p能夠經過把p-1寫成二進制，由(a*b)mod c=(a mod c)*b mod c，能夠在t=log(p-1)的時間內計算出解，如考慮整數相乘的複雜度，則一次計算的總複雜度爲log³(p-1)。這個方法叫快速冪取模。

爲了提升算法的準確性，咱們又有一個能夠利用的定理。
定理二：對於0<x<p，x^2 mod p =1 => x=1或p-1。

咱們令p-1=(2^t)*u，即p-1爲u二進制表示後面跟t個0。咱們先計算出x[0]=a^u mod p ，再平方t次並在每一次模p，每一次的結果記爲x[i]，最後也能夠計算出a^(p-1) mod p。若發現x[i]=1而x[i-1]不等於1也不等於p-1，則發現p果斷不是素數。

能夠證實，使用以上兩個定理之後，檢驗s次出錯的機率至多爲2^(-s)，因此這個算法是很可靠的。

須要注意的是，爲了防止溢出（特別大的數據），a*b mod c 也應用相似快速冪取模的方法計算。固然，數據不是很大就能夠免了。

下面是個人程序。

typedef unsigned long long LL;  
   
LL modular_multi(LL x,LL y,LL mo)  
{  
    LL t;  
    x%=mo;  
    for(t=0;y;x=(x<<1)%mo,y>>=1)  
        if (y&1)  
            t=(t+x)%mo;  
    return t;  
}  
   
LL modular_exp(LL num,LL t,LL mo)  
{  
    LL ret=1,temp=num%mo;  
    for(;t;t>>=1,temp=modular_multi(temp,temp,mo))  
        if (t&1)  
            ret=modular_multi(ret,temp,mo);  
    return ret;  
}  
   
bool miller_rabbin(LL n)  
{  
    if (n==2)return true;  
    if (n<2||!(n&1))return false;  
    int t=0;  
    LL a,x,y,u=n-1;  
    while((u&1)==0) t++,u>>=1;  
    for(int i=0;i<S;i++)  
    {  
        a=rand()%(n-1)+1;  
        x=modular_exp(a,u,n);  
        for(int j=0;j<t;j++)  
        {  
            y=modular_multi(x,x,n);  
            if (y==1&&x!=1&&x!=n-1)  
                return false;  
            ///其中用到定理，若是對模n存在1的非平凡平方根，則n是合數。  
            ///若是一個數x知足方程x^2≡1 (mod n),但x不等於對模n來講1的兩個‘平凡’平方根：1或-1，則x是對模n來講1的非平凡平方根  
            x=y;  
        }  
        if (x!=1)///根據費馬小定理,若n是素數，有a^(n-1)≡1(mod n).所以n不多是素數  
            return false;  
    }  
    return true;  
}