用Monte Carlo算法实现素数判定

一、问题描述
判定一个给定的整数是否为素数

二、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 <cstdio>
#include <cmath>
#include <random>
#include <ctime>
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))+ a
int 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)运行结果