根據費馬小定理，若p為素數，則必有a^(p-1) mod p=1 對和p互質的a成立。

根據二次探測定理：如果p是素數，且0<x<p，則方程x^2 mod p=1的解為1或p-1。

所以若p為素數，則必有a^(p-1) mod p 的平方根為1或-1

分解p-1為d*2^s，其中d為奇數

從i=0逐次計算a^(d*2^(s-i))，相當於“開方”，若得到-1或追查到a^d=1 (mod p)，則p通過測試，否則不通過

時間復雜度O(k*(logn)^3) （其中k為選的a的個數（the more the better?））

 1 #include<cstdio>
 2
 
 #include<cstring>
 3 #include<iostream>
 4 using namespace std;
 5 
 6 const int prime[9]={2,3,5,7,11,13,17,19,23};
 7 
 8 int n;
 9 
10 int quick(int a,int b,int mod){
11     int sum=1;
12     for(;b;b>>=1,a=a*a%mod)
13         if(b&1)  sum=sum*a%mod;
14     return
 
 sum;
15 }
16 
17 bool Rabin_Miller(int p,int a){
18     if(p==2)  return 1;
19     if((p&1)==0||p==1)  return 0;
20     int d=p-1;
21     while((d&1)==0)  d>>=1;
22     int m=quick(a,d,p);
23     if(m==1)  return 1;
24     for(;d<p;d<<=1,m=m*m%p)
25         if(m==p-1)  return
 
 1;
26     return 0;
27 }
28 
29 bool isprime(int x){
30     for(int i=0;i<9;i++){
31         if(x==prime[i])  return 1;
32         if(!Rabin_Miller(x,prime[i]))  return 0;
33     }
34     return 1;
35 }
36 
37 int main(){
38     scanf("%d",&n);
39     if(isprime(n))  puts("Yes!");
40     else  puts("No!");
41     return 0;
42 }

Miller-Rabin 素性測試