问题
如题
思路
首先太暴力的就不谈,会折寿,有一个强伪证的算法(strong liar)Miller Rabin Test,思路是这样(以下p代表某一个大于2的素数),找到判断素数的必要条件,将多个条件组合起来就能得到近似充分的判断素数条件
两个引理说一下,首先是二次探测定理:如果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互质(gcd(a,p)=1)时,ap-1≡1(mod p)
好了开始,对于要判断的数P,任取2-p内的数a,这里我们方便起见取a=2
根据上述条件,p-1一定是偶数,可以分解成2r*d的形式,所以对于a2^r*d ≡1(mod p),也就是(a2^(r-1)*d )2≡1(mod p),设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),然后我们可以继续往下解xn-1,xn-2...
如果是x(r-1)≡-1(mod p),那么问题就更简单了,因为接下来的解不确定,不用往下推了
所以从上往下推我们可以看到,对于函数x的序列,x(r)=1是一定的,而且如果x(t)=-1,那么x(t+1)到x(r)一定都是1,也就是说,对于序列中的数只要判断两个条件即可:
1.1之前的数一定是1或-1
2.最后一个数一定是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之前的数一定是1或者-1
prev = cur;
u >>= 1;
}
if (cur != 1)
return false;//序列的最后一个数一定是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范围内所有数(小于2147483648的所有数)都可以判断出结果,注意这时需要对3,5,7,11单独判断一遍,因为前提条件是gcd(a,p)=1
- 并不是所有数据类型取long就能判断出长整型范围内的数据了,因为快速幂中需要计算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