引言：

如果你遇見很多個非常大的數n，有10的18次方這麼大，讓你來判斷是不是素數，你有什麼快速的方法呢？

用普通的判斷會超時，用篩法又篩不到，這時候就可以用上大素數判斷演算法。

（其實這個演算法理解原理後很簡單，我也就wa了11次）

原理部分：

首先我們知道，對於一個質數P，他一定滿足費馬小定理，

即 a^(p-1)≡1(mod p)

那麼我們反過來想，對於任意一個數n，若它滿足 a^(n-1)≡1(mod n)，那麼它是不是就是一個質數呢?

答案是 不一定，但是 在絕大多數的狀況下，n就是個質數。因此，對於一個很大的數n，我們可以隨便取一個a,然後用

快速冪計算a^(n-1)mod n,看最後的結果是不是1，以此來判斷是不是質數。

但這樣還是有很大的誤差，為此我們引入一個新定理——二次探測原理。

若P是一個質數，且滿足0<x<p、x^2 mod p =1，則可以推斷出x一定等於1或p-1。

因為除了2，所有的質數都是奇數，所以n若是偶數且不為2，那他一定不是質數，若他是奇數，則（n-1）一定是偶數,

我們就可以把a^(n-1)改寫成(a^{^(n-1)/2})²,就滿足了二次探測原理的x^2的形式，我們就可以用二次探測原理來驗證n是不是質

數，即若x^2%n=1,但x不等於1或n-1，則n一定不是一個質數。

更關鍵的是，若（n-1）/2仍是一個偶數，我們還可以把a^{^(n-1)/2}改寫乘某個數的平方，以此類推。

即把n-1化為2^s*d , a^(n-1)=

我們就可以在求a^(n-1)mod n的過程中，多次驗證二次探測原理。

通過上面兩個原理，我們再多用幾個a驗證，判斷n是不是質數準確率就可以達到99.9999999999999....%。

另外有一個結論：如果我們選擇2,3,7,61,24251,9875321這6個質數作為a的話，在10^14內不存在誤判。

好了，如果沒有看懂上面的的理論就看下面的程式碼吧。

程式碼部分：

inline LL quick_mul(LL a,LL b,LL m) //快速乘 
{
    ll ans = 0;
    a %= m;
    b %= m;
    
 
while (b) {
        if (b & 1) {
            ans = (ans + a) % m;
        }
        a = (a + a) % m;
        b >>= 1;
    }
    return ans;
}

LL quick_pow(LL a,LL b,LL m) //快速冪 
{
    LL res=1;
    a%=m;
    while(b)
    {
        if(b&1) res=quick_mul(res,a,m);
        a=quick_mul(a,a,m);
        b>>=1;
    }
    return res;
}
bool Miller_rabin(LL n,int num)
{
//先將特殊情況判斷一下
    if(n==2||n==3) return true;
    if(n%2==0||n==1) return false;
    srand((unsigned)time(NULL)); //為接下來a的隨機取值用 
    LL d=n-1;
    int s=0;
    while(!(d&1)) s++,d>>=1;//若d的二進位制的最後一位不是1，則說明d還是偶數 
    for(int i=0;i<num;i++)
    {
        LL a=rand()%(n-2)+2;//2~n-1;
        LL x=quick_pow(a,d,n), y=0;
        for(int j=0;j<s;j++)//一共平方s次 
        {
            y=quick_mul(x,x,n);//先平方 
            if(y==1&&x!=1&&x!=(n-1)) return false;//驗證二次探測原理 
            x=y;
        }
        if(y!=1) return false;//不滿足費馬小定理，那就肯定不是質數
    }
    return true;
}

另外補充一點就是，對於上面的快速乘，還有一種更快速的寫法

inline LL quick_mul(LL a,LL b,LL m)
{
    LL t=(ld)a/m*b;
    LL res=(ull)a*b-(ull)t*m;
    return (res+m)%m;
}

利用的是 相乘溢位的是高位，但由於取模我們只需要的是低位。

還有就是a*b%m=a*b-(a*b/m)*m