普通的素数测试咱们有O(√ n)的试除算法。事实上,咱们有O(slog³n)的算法。算法
定理一:假如p是质数,且(a,p)=1,那么a^(p-1)≡1(mod p)。即假如p是质数,且a,p互质,那么a的(p-1)次方除以p的余数恒等于1。(费马小定理)测试
该定理的逆命题是不必定成立的,可是使人可喜的是大多数状况是成立的。spa
因而咱们就获得了一个定理的直接应用,对于待验证的数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)。这个方法叫快速幂取模。code
为了提升算法的准确性,咱们又有一个能够利用的定理。
定理二:对于0<x<p,x^2 mod p =1 => x=1或p-1。orm
咱们令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果断不是素数。it
能够证实,使用以上两个定理之后,检验s次出错的几率至多为2^(-s),因此这个算法是很可靠的。class
须要注意的是,为了防止溢出(特别大的数据),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; }