為了測試一個大整數是不是素數，我們不能夠使用傳統的測試是否有因子的方法，因為那樣的時間復雜度至少也是$O(n)$，空間復雜度是$O(n)$（使用線性篩數法），時間復雜度還好說，空間復雜度在n比較大的時候完全不可以接受，這就需要我們使用$Miller\_Rabin$算法進行解決。

算法的核心使用了費馬小定理和二次探測定理

費馬小定理：若$p$為質數，則$a^{p?1}\equiv 1(\mathrm{m}\mathrm{o}\mathrm{d}p)$

這是一個質數的必要條件，但是有一類特殊的數也滿足這個等式并且是合數（卡米歇爾數），為了剔除這些特殊的數，我們還需要二次探測定理幫助判斷

二次探測定理：如果$p$為素數，$x^2\equiv 1(\mathrm{m}\mathrm{o}\mathrm{d}p)$，則$x=1或x=p?1$

因此對于滿足費馬小定理的數$p$，即${a^{\frac{p?1}{2}}}^2\equiv 1(\mathrm{m}\mathrm{o}\mathrm{d}p)$，我們判斷$a^{\frac{p?1}{2}}$是否為1或者$p?1$，如果不是說明$p$不是素數，如果是我們可以繼續對$a^{\frac{p?1}{2}}$進行探測

每一次檢測都無法100%確定這是一個素數，但是判斷8次幾乎不可能出錯。

這里的mult函數是乘積函數，因為接近long long 的數字相乘會爆long long,這里按位相乘。

在研究的很多人的實現后我實現了以下，覺得我這個復雜度低一些(因為自己寫的緣故看起來很親切)。實際過題時間也比之前看到的快。

typedef long long ll;
ll mult(ll x,ll y,ll p)
{long double d=1; d=d*x/p*y;return ((x*y-((ll)d)*p)%p+p)%p;
}ll quick_pow(ll a,ll b,ll p)
{a%=p; ll ret=1;while(b){if(b&1) ret=mult(ret,a,p);a=mult(a,a,p); b>>=1;}return ret;
}bool Miller_Rabin(ll n)
{if(n<2) return false;const int times=8;const ll prime[times]={2,3,5,7,11,13,17,61};for(int i=0;i<times;i++)if(n==prime[i]) return true;else if(n%prime[i]==0) return false;ll x=n-1; while(~x&1) x>>=1;for(int i=0;i<times;i++){ll a=prime[i]; ll now=quick_pow(a,x,n); ll last;if(x==n-1){if(now!=1) return false;}else{while(x!=n-1){x<<=1; last=now; now=mult(now,now,n);if(now==1){if(last!=1 && last!=n-1) return false;break;}}}}return true;
}

之前看到的別人的板子

ll mult(ll a, ll b, ll p) {a %= p;b %= p;ll res = 0, tmp = a;while(b) {if (b & 1) {res += tmp;if (res > p) res -= p;}tmp <<= 1;if (tmp > p) tmp -= p;b >>= 1;}return res;
}ll quick_pow(ll a, ll b, ll p) {ll res = 1;ll tmp = a % p;while(b) {if (b & 1) res = mult(res, tmp, p);tmp = mult(tmp, tmp, p);b >>= 1;}return res;
}bool check(ll a, ll n, ll x, ll t) {ll res = quick_pow(a, x, n);ll last = res;for (ll i = 1; i <= t; i++) {res = mult(res, res, n);if (res == 1 && last != 1 && last != n-1) return true;last = res;}return res != 1;
}bool Miller_Rabin(ll n) {if (n < 2) return false;if (n == 2) return true;if ((n & 1) == 0) return false;ll x = n-1;ll t = 0;while ((x & 1) == 0) {x >>= 1;t++;}srand(time(NULL));const ll tims = 8;for (ll i = 0; i < tims; i++) {ll a = rand() % (n-1) + 1;if (check(a, n, x, t)) return false;}return true;
}