問題
如題
思路
首先太暴力的就不談饺蔑，會折壽掠河，有一個強偽證的算法（strong liar）Miller Rabin Test蜈垮，思路是這樣（以下p代表某一個大于2的素數(shù)），找到判斷素數(shù)的必要條件菊碟，將多個條件組合起來就能得到近似充分的判斷素數(shù)條件
兩個引理說一下节芥，首先是二次探測定理：如果x2≡1(mod p)，那么x≡±1(mod p)逆害，或者說不存在模p等于1下的非平凡平方根（證明：p只有1和p兩個因子的情況下头镊，(x+1)(x-1)≡0(mod p)只可能有x-1=0或x+1=p）
然后是費馬小定理，a和p互質(zhì)（gcd(a,p)=1）時魄幕，ap-1≡1(mod p)
好了開始相艇，對于要判斷的數(shù)P，任取2-p內(nèi)的數(shù)a纯陨，這里我們方便起見取a=2
根據(jù)上述條件坛芽，p-1一定是偶數(shù)，可以分解成2r*d的形式翼抠，所以對于a2^r*d ≡1(mod p)咙轩，也就是(a2^(r-1)*d )2≡1(mod p)，設(shè)x(n) = (a(2^n*d)) mod p阴颖，那么目前的情況是活喊，我們僅有平凡平方根x(r-1)≡±1(mod p)
這里如果x(r-1)≡1(mod p)，問題就變得可遞歸了量愧，因為接下來的方程是a2^(r-1)*d≡1(mod p)钾菊，然后我們可以繼續(xù)往下解xn-1,xn-2...
如果是x(r-1)≡-1(mod p)帅矗，那么問題就更簡單了，因為接下來的解不確定煞烫，不用往下推了
所以從上往下推我們可以看到浑此，對于函數(shù)x的序列，x(r)=1是一定的滞详，而且如果x(t)=-1凛俱，那么x(t+1)到x(r)一定都是1，也就是說茵宪，對于序列中的數(shù)只要判斷兩個條件即可：
1.1之前的數(shù)一定是1或-1
2.最后一個數(shù)一定是1
解決
public boolean sPrime_RobinMiller(long p) {
long[] ar = new long[]{2, 3, 5, 7, 11};//if n < 2,152,302,898,747, it is enough to test a = 2, 3, 5, 7, and 11.
if (p < 4) return p > 1;
if (p == 5 || p == 7 || p == 11) return true;//a<p-1
for (long a : ar) {
long p1 = p - 1, u = p1 & -p1, d = p1 / u, cur = quickPow(a, d, p), prev = cur;
while (u > 0) {
cur = prev * prev % p;
if (cur == 1 && prev != 1 && prev != p1)
return false;//序列中最冰，1之前的數(shù)一定是1或者-1
prev = cur;
u >>= 1;
}
if (cur != 1)
return false;//序列的最后一個數(shù)一定是1
}
return true;//all passed
}
private long quickPow(long a, long b, long c) {
long res = 1;
for (; b > 0; b >>= 1, a = a*a%c)
if ((b & 1) == 1) {
res = res * a % c;
}
return res;
}
Tips
- 取a=2,3,5,7,11，對于int范圍內(nèi)所有數(shù)（小于2147483648的所有數(shù)）都可以判斷出結(jié)果稀火，注意這時需要對3,5,7,11單獨判斷一遍暖哨，因為前提條件是gcd(a,p)=1
- 并不是所有數(shù)據(jù)類型取long就能判斷出長整型范圍內(nèi)的數(shù)據(jù)了，因為快速冪中需要計算2d mod p凰狞，其中a*a是有可能溢出的
Ref
https://oi-wiki.org/math/prime/
https://zh.wikipedia.org/wiki/%E7%B1%B3%E5%8B%92-%E6%8B%89%E5%AE%BE%E6%A3%80%E9%AA%8C
