用Monte Carlo算法实现素数判定

墨蓝 2022-11-22 00:15 193阅读 0赞

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

二、Monte Carlo算法  
1.预备知识  
(1)费马小定理  
![在这里插入图片描述][20201103142924945.png_pic_center]  
其逆否命题：  
![在这里插入图片描述][20201103143113780.png_pic_center]  
利用该逆否命题进行素数判定(代码如下)：

//利用费马小定理(正确的概率较低)
    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 ∗ t n - 1 = 2^s \* t n−1=2s∗t的正整数，其中t为奇数。  
当且仅当 a ∈ \[ 2 , n − 2 \] a\\in\[2, n - 2\] a∈\[2,n−2\]且满足下述2个条件之一：  
① a t    %    n = = 1 a^t \\;\\%\\; n == 1 at%n==1  
② ∃ 整 数 i ∈ \[ 0 , s − 1 \] \{\\exists\}整数i\\in\[0, s - 1\] ∃整数i∈\[0,s−1\]，满足 a 2 i ∗ t    %    n = = n − 1 a^\{2^i \* t\}\\;\\%\\;n==n-1 a2i∗t%n==n−1  
则称n为一个以a为底的**强伪素数**，称a为n素性的**强伪证据**，即 a ∈ B ( n ) a\\in B(n) a∈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是否为强伪素数
    }

![在这里插入图片描述][watermark_type_ZmFuZ3poZW5naGVpdGk_shadow_10_text_aHR0cHM6Ly9ibG9nLmNzZG4ubmV0L3dlaXhpbl8zNzQ3NzAwOQ_size_16_color_FFFFFF_t_70_pic_center]  
**分析**：1次调用MillRab函数，若返回false，一定为合数；若返回true，则判定为合数的概率< 1 4 \\frac\{1\}\{4\} 41。为缩小误判概率，调用k次MillRab函数，若每次都返回true，则判定为合数的概率< 1 4 k \\frac\{1\}\{4\}^k 41k，当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)运行结果  
![在这里插入图片描述][watermark_type_ZmFuZ3poZW5naGVpdGk_shadow_10_text_aHR0cHM6Ly9ibG9nLmNzZG4ubmV0L3dlaXhpbl8zNzQ3NzAwOQ_size_16_color_FFFFFF_t_70_pic_center 1]

[20201103142924945.png_pic_center]: /images/20221120/5259828a8c1041918b4d4c137d009012.png
[20201103143113780.png_pic_center]: https://img-blog.csdnimg.cn/20201103143113780.png#pic_center
[watermark_type_ZmFuZ3poZW5naGVpdGk_shadow_10_text_aHR0cHM6Ly9ibG9nLmNzZG4ubmV0L3dlaXhpbl8zNzQ3NzAwOQ_size_16_color_FFFFFF_t_70_pic_center]: https://img-blog.csdnimg.cn/20201103144741523.png?x-oss-process=image/watermark,type_ZmFuZ3poZW5naGVpdGk,shadow_10,text_aHR0cHM6Ly9ibG9nLmNzZG4ubmV0L3dlaXhpbl8zNzQ3NzAwOQ==,size_16,color_FFFFFF,t_70#pic_center
[watermark_type_ZmFuZ3poZW5naGVpdGk_shadow_10_text_aHR0cHM6Ly9ibG9nLmNzZG4ubmV0L3dlaXhpbl8zNzQ3NzAwOQ_size_16_color_FFFFFF_t_70_pic_center 1]: https://img-blog.csdnimg.cn/20201103145839305.png?x-oss-process=image/watermark,type_ZmFuZ3poZW5naGVpdGk,shadow_10,text_aHR0cHM6Ly9ibG9nLmNzZG4ubmV0L3dlaXhpbl8zNzQ3NzAwOQ==,size_16,color_FFFFFF,t_70#pic_center