转载自：http://www.dxmtb.com/blog/miller-rabbin/
普通的素数测试我们有O(√ n)的试除算法。事实上，我们有O(slog³n)的算法。
定理一：假如p是质数，且(a,p)=1，那么a^(p-1)≡1(mod p)。即假如p是质数，且a,p互质，那么a的(p-1)次方除以p的余数恒等于1。（费马小定理）
该定理的逆命题是不一定成立的，但是令人可喜的是大多数情况是成立的。
于是我们就得到了一个定理的直接应用，对于待验证的数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)。这个方法叫快速幂取模。
为了提高算法的准确性，我们又有一个可以利用的定理。
定理二：对于0<x<p，x^2 mod p =1 => x=1或p-1。
我们令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果断不是素数。
可以证明，使用以上两个定理以后，检验s次出错的概率至多为2^(-s)，所以这个算法是很可靠的。
需要注意的是，为了防止溢出（特别大的数据），a*b mod c 也应用类似快速幂取模的方法计算。当然，数据不是很大就可以免了
codevs 1702 素数判定 2
一个数，他是素数么？
设他为P满足(P<=2^63-1)
P
Yes|No
2
Yes
算法导论——数论那一节
注意Carmichael Number
分类标签 Tags 点此展开
1 /*数据确实没有错，我的错误点是，在两个long long相乘的时候，直接进行了乘法运算，很有可能溢出，所以错了后面的点，long long数据要再写一个相乘%的子函数*/ 2 #include<iostream> 3 using namespace std; 4 #include<cstdio> 5 #define S 10 6 #include<cstdlib> 7 #include<ctime> 8 #define ll long long 9 ll cas, maxz; 10 ll read() 11 { 12 ll ans=0;char c; 13 c=getchar(); 14 while(c<'0'||c>'9') c=getchar(); 15 while(c>='0'&&c<='9') 16 { 17 ans=ans*10+c-'0'; 18 c=getchar(); 19 } 20 return ans; 21 } 22 ll quick_mul_mod(ll a,ll b,ll c)//a*b%c 23 { 24 ll ret=0; 25 a%=c;b%=c; 26 while(b) 27 { 28 if(b&1) 29 { 30 ret+=a;/*相当于加了b次a*/ 31 ret%=c; 32 b--; 33 } 34 a<<=1; 35 a%=c; 36 b>>=1; 37 } 38 return ret; 39 } 40 ll quick_mod(ll a,ll b,ll c)//ji suan a^b%c 41 { 42 ll ans=1; 43 a%=c; 44 while(b) 45 { 46 if(b&1) 47 { 48 b--; 49 ans=quick_mul_mod(ans,a,c);/*注意只要是long long的乘法，很有可能溢出的*/ 50 } 51 b>>=1; 52 a=quick_mul_mod(a,a,c); 53 } 54 return ans; 55 } 56 bool Miller_rabin(ll n) 57 { 58 if(n==2) return true; 59 if(n<=1||!(n&1)) return false; 60 ll u=n-1,t=0; 61 while(!(u&1)) 62 { 63 u>>=1; 64 t++; 65 } 66 for(int i=0;i<S;++i) 67 { 68 ll x=rand()%(n-1)+1; 69 x=quick_mod(x,u,n); 70 for(int i=1;i<=t;++i) 71 { 72 ll y=quick_mul_mod(x,x,n); 73 if(y==1&&x!=1&&x!=n-1) 74 return false; 75 x=y; 76 } 77 if(x!=1) return false; 78 } 79 return true; 80 } 81 82 int main() 83 { 84 ll n=read(); 85 if(Miller_rabin(n)) 86 printf("Yes "); 87 else printf("No "); 88 return 0; 89 }