本文共 4494 字，大约阅读时间需要 14 分钟。

一、问题描述

判定一个给定的整数是否为素数

二、Monte Carlo算法

1.预备知识 (1)费马小定理 其逆否命题： 利用该逆否命题进行素数判定(代码如下)：

//利用费马小定理(正确的概率较低)void Fermat(int n) {
       bool flag = true;    srand((unsigned)time(0));    int a = uniform(1, n); //[1, n-1]    if(Pow(a, n - 1) % n != 1)        flag = false;    Print(flag, n);  //flag = true未必是素数}

缺陷：

即使任意的a属于[1, n-1]，使得Pow(a, n - 1) % n = =1 也不能保证n为素数。(费马小定理的逆命题不成立)

结论：

不能通过验证Pow(a, n - 1) % n = =1来判定n是否为素数。

2.对Fermat测试改进

(1)理论知识与对应代码 设n是一个大于4的奇整数，s和t是使得$n-1=2^s\ast t$的正整数，其中t为奇数。 当且仅当$a\in [2,n-2]$且满足下述2个条件之一： ①$a^t\%n==1$ ②$\exists 整数i\in [0,s-1]$，满足$a^{2^i\ast t}\%n==n-1$ 则称n为一个以a为底的强伪素数，称a为n素性的强伪证据，即$a\in B(n)$

代码如下：

//对Fermat的改进//判定a是否为n素性的强伪证据bool Btest(int a, int n) {
     //n为奇数    int s = 0, t = n - 1;    do    {
           s++;        t /= 2;    } while (t % 2 != 1); //满足n - 1 = 2^s * t的最大奇数t;    //int x = Pow(a, t) % n;  Pow(a, t)溢出了。    int x = 1;    while(t--) {
           x = (x * (a % n)) % n;  //(a * b) % p = (a % p * b % p) % p    }    if(x == 1 || x == n - 1)        return true;    for (int i = 1; i < s; i++) {
           x = (x * x) % n;        if(x == n - 1)            return true;    }    return false;  //一定是合数}bool MillRab(int n) {
       //srand((unsigned)time(0));   //放这误差很大    int a = uniform(2, n - 1); //[2, n-2]    return Btest(a, n);  //测试n是否为强伪素数}

分析：1次调用MillRab函数，若返回false，一定为合数；若返回true，则判定为合数的概率<$\frac{1}{4}$。为缩小误判概率，调用k次MillRab函数，若每次都返回true，则判定为合数的概率<${\frac{1}{4}}^k$，当k取10时，误判概率低于百万分之一。

代码如下：

bool RepeatMillRab(int n, int k) {
       for (int i = 1; i <= k; i++) {
           if(MillRab(n) == false)            return false;    }    return true;}

三、完整代码及运行结果

(1)完整代码

//素数的判定#include 
   
    #include 
    
     #include 
     
      #include 
      
       using namespace std;typedef long long LL;void Print(bool flag, int n) {
       if(flag)        printf("%d is a prime.\n", n);    else        printf("%d is not a prime.\n", n);}//传统方法void TraditionalPrime(int n) {
       bool flag = true; //假设一开始为素数    int upper = floor(sqrt(n));    for (int i = 2; i <= upper; i++) {
           if(n % i == 0)            flag = false;    }    Print(flag, n);}bool TraditionalMethod(int n) {
       int upper = floor(sqrt(n));    for (int i = 2; i <= upper; i++) {
           if(n % i == 0)            return false;    }    return true;}//[a, b-1] (rand() % (b-a))+ aint uniform(int a, int b) {
       return (rand() % (b - a)) + a;}//整数幂LL Pow(int a, int n) {
       LL sum = 1;    for (int i = 1; i <= n; i++) {
           sum *= a;    }    return sum;}//利用费马小定理(正确的概率较低)void Fermat(int n) {
       bool flag = true;    srand((unsigned)time(0));    int a = uniform(1, n); //[1, n-1]    if(Pow(a, n - 1) % n != 1)        flag = false;    Print(flag, n);  //flag = true未必是素数}//对Fermat的改进//判定a是否为n素性的强伪证据bool Btest(int a, int n) {
     //n为奇数    int s = 0, t = n - 1;    do    {
           s++;        t /= 2;    } while (t % 2 != 1); //满足n - 1 = 2^s * t的最大奇数t;    //int x = Pow(a, t) % n;  Pow(a, t)溢出了。    int x = 1;    while(t--) {
           x = (x * (a % n)) % n;  //(a * b) % p = (a % p * b % p) % p    }    if(x == 1 || x == n - 1)        return true;    for (int i = 1; i < s; i++) {
           x = (x * x) % n;        if(x == n - 1)            return true;    }    return false;  //一定是合数}bool MillRab(int n) {
       //srand((unsigned)time(0));   //放这误差很大    int a = uniform(2, n - 1); //[2, n-2]    return Btest(a, n);  //测试n是否为强伪素数}bool RepeatMillRab(int n, int k) {
       for (int i = 1; i <= k; i++) {
           if(MillRab(n) == false)            return false;    }    return true;}//打印1万以内的素数void PrintPrimes() {
       int pn = 2, count = 2; //pn表示真素数的个数, count表示判定为素数的个数(含假阳性)；    printf("2 3 ");    int n = 5;    do{
           if(RepeatMillRab(n, floor(log10(n)))) {
               if(TraditionalMethod(n))                pn++;            count++;            printf("%d", n);            if(count % 20 == 0)                printf("\n");            else                printf(" ");        }        n += 2;    } while (n < 10000);    printf("误判个数：%d", count - pn);    double correct = (double)(pn) / count * 100;    printf("\n正确比例：%f%%\n", correct);}int main() {
       //传统方法    //int n = 2643;    //TraditionalPrime(n);    //Fermat(n);    srand((unsigned)time(0));          PrintPrimes();    return 0;}
      
     
    
   

(2)运行结果

转载地址：http://cthhz.baihongyu.com/