素数的测定

强伪素数的定义，如下：

则容易知道，满足条件的n有可能是素数也有可能是合数(强伪素数)，但是不满足条件的一定是合数。我们又有以下定理

所以当返回为真时，它为合数错误的概率小于1/4，当返回为假时，必为合数，那么重复调用k次返回为真时，该数为合数的概率应该小于$\frac{1}{4^k}$，只要k取10，错误概率就小于百万分之一。
整个代码如下：

#include<iostream>
#include<cstdlib>
#include<ctime>
#include<cmath>

using namespace std;

bool is_primes(int n){
    if(n==1)
        return false;       //1不是素数
    else{
        for(int i=2;i<sqrt(n)+1;i++)
            if(n%i==0)
                return false;            //合数
    }
    return true;    //素数
}

bool Btest(int a,int n){//返回true表示是强伪素数，返回false表示是合数
    int s=0;
    int t=n-1;
    while(t%2==0){
        t=t/2;
        s=s+1;
    }
    //计算a^t mod n
    int i=0;
    int x=1;
    while(i++<t){
        x=(x*a)%n;
    }
    if(x==1||x==n-1)
        return true;
    for(i=1;i<=s-1;i++){
        x=(x*x)%n;
        if(x==n-1)
            return true;
    }
    return false;
}

bool MillRab(int n){
    int a=rand()%(n-3)+2;  //a=2,3,...,n-4+2=n-2.
    return Btest(a,n);
}

bool RepeatMillRob(int n,int k){
    for(int i=0;i<k;i++)
        if(MillRab(n)==false)
             return false;          //一定是合数
    return true;
}

int PrintPrimes(){
   cout<<2<<endl<<3<<endl;
   int n=5;
   int err=0;
   while(n<=10000){
        int k=log10(n);
        if(RepeatMillRob(n,k)){
            if(is_primes(n))
                cout<<n<<endl;
            else
                err++;
        }
        n=n+2;
   }
   return err;
}

int main(){
    srand(time(0));
    cout<<PrintPrimes()<<endl;
    return 0;
}

运行结果：

在5-10000的素数检测中，用概率方法仅错了2个，增大k能进一步提高准确率。