關于素數：求不超過n的素數，素數的判定(Miller Rabin 測試)

關于素數的基本介紹請參考百度百科here和維基百科here的介紹

首先介紹幾條關于素數的基本定理：

定理1：如果n不是素數，則n至少有一個( 1, sqrt(n) ]範圍内的的因子

定理2：如果n不是素數，則n至少有一個(1, sqrt(n) ]範圍内的素數因子

定理3：定義f(n)為不大于n的素數的個數，則 f(n) 近似等于 n/ln(n) (ln為自然對數) ，具體請參考here

求不超過n的素數本文位址

算法1：埃拉托斯特尼篩法，該算法的核心思想是：首先标記2~n的數都為素數，然後周遊2~n的數組，如果它是素數，就把它的倍數标記為非素數（即把所有素數的倍數都标記為非素數）代碼如下：

1 unsigned int findPrime(const unsigned int n, unsigned int prime[])
 2 {
 3     if(n < 2)return 0;
 4     unsigned int primeNum = 0;
 5     bool * isPrime = NULL;
 6     //如果n太大，下面可能會申請記憶體失敗
 7     try{
 8         isPrime = new bool[n+1];//0表示是素數，1表示非素數
 9     }catch(const std::bad_alloc &e){
10         std::cout<<"bad_alloc: new failed\n";
11         return -1;
12     }
13     memset(isPrime, 0, sizeof(bool)*(n+1));
14     for(long long i = 2; i <= n; i++)
15         if(isPrime[i] == 0)
16         {//所有素數的倍數都不是素數
17             prime[primeNum++] = i;
18             long long t ;
19             for(unsigned int k = 2; (t = i*k) <= n; k++)
20                 isPrime[t] = 1;
21         }
22     delete []isPrime;
23     return primeNum;//傳回素數的個數
24 }

算法2：查表法，該算法主要根據定理2，如果一個數是素數，那麼它不能被(1, sqrt(n) ]範圍内的素數整除，稱之為查表法，主要是因為當我們判斷n是否為素數時，(1, sqrt(n) ]範圍内的素數已經都計算出來了，可以利用已經計算出來的素數表，代碼如下：

1 unsigned int findPrime_table(const unsigned int n, unsigned int prime[])
 2 {
 3     if(n < 2)return 0;
 4     unsigned int primeNum = 0;
 5     prime[primeNum++] = 2;
 6     for(long long i = 3; i <= n; i++)
 7     {
 8         unsigned int tmp  = sqrt(i);
 9         bool flag = true;
10         for(unsigned int j = 0; prime[j] <= tmp; j++)
11         {
12             if(i % prime[j] == 0)
13             {
14                 flag = false;
15                 break;
16             }
17         }
18         if(flag)
19             prime[primeNum++] = i;
20     }
21     return primeNum;//傳回素數的個數
22 }

對于兩個算法的效率，我們分别計算不超過100,500,1000,10000,100000,1000000,10000000的素數，用時如下，左邊是算法1的耗時，右邊是算法2的耗時，機關微妙（測試環境，win7 x64 7G記憶體，i5處理器，codeblock with gcc）：

1.65536，4.635

2.64857，16.2225

4.30393，26.1546

65.5521，313.524

607.847，5474.92

5904.66，74654.9

212453，1.38515e+006

通過上面的資料可知，算法1的性能優于算法2

在記憶體和時間允許的情況下，上面提供的代碼（在int是32位的機器上）理論上能計算2^32-1以内的所有素數

素數的判定本文位址

算法1：最naive的方法，根據定理一來判斷，代碼如下：

1 bool isPrime_naive(const unsigned int n)
2 {
3     if(n < 2)return false;
4     unsigned int k = sqrt(n);
5     for(unsigned int i = 2; i <= k; i++)
6         if(n%i == 0)
7             return false;
8     return true;
9 }

算法2：根據定理2，先調用上面的素數生成算法生成sqrt(n)以内的素數（由于上面已經證明篩選法較快，是以下面判定算法中調用篩選法來生成素數表），然後再看他們是否都不能被n整除，代碼如下：

bool isPrime_checkList(const unsigned int n)
{
    if(n < 2)return false;
    if(n == 2 || n == 3) return true;
    unsigned int tmp = sqrt(n);
    unsigned int *prime = new unsigned int[tmp+1];
    unsigned int primeNum = findPrime(tmp, prime);//生成sqrt(n)以内的素數
    for(unsigned int i = 0; prime[i] <= tmp && i < primeNum; i++)
        if(n%prime[i] == 0)
        {
            delete []prime;
            return false;
        }
    delete []prime;
    return true;
}

算法3：機率算法，Miller Rabin 素數測試（參考資料：資料一，資料二，資料三），先講幾個有關的定理

費馬小定理：假如p是素數，且(a,p)=1，那麼 a^(p-1) ≡1（mod p）。即：假如p是素數，且a,p互質，那麼a的(p-1)次方除以p的餘數恒等于1。

二次探測定理：若 p 為素數，且 0 < x < p，則方程 x * x = 1 ( mod p ) 的解為 x = 1, p - 1 ( mod p )

在以上兩個定理的基礎上得出Miller Rabin素數測試的理論基礎：如果n是一個奇素數，将n - 1 分解成 ( 2^s ) * d，其中 s 是非負整數，d 是正奇數，a 是和n互素的任何整數，那麼a^d≡1(mod n) 成立或者對某個r(0≤r ≤s －1) 等式 a^(2^r*d) ≡－1(mod n)成立

那麼我們如何根據上述理論來測試某個數是否是素數呢，方法如下：對于一個數n, 将n - 1 分解成 ( 2^s ) * d，其中 s 是非負整數，d 是正奇數, 随機選取[1, n-1]内的整數a, 如果a^d≡1(mod n) 并且對所有的r(0≤ r ≤s －1) a^(2^r*d) ≡－1(mod n),那麼n是合數，否則n有3/4的機率是素數（關于3/4這個機率的證明見here）

按照上述方法一次判斷，把某個數誤判成素數的機率1/4，但是不會把某個素數誤判成合數，如果我們運作t次上述判斷，那麼誤判的機率是 (1/4)^t，當t = 10時這個機率是百萬分之一，是以我們總結出miller rabin素數測試的算法步驟如下：

----------------------------------------------------------

算法輸入：某個大于3的奇數n, 測試輪數t

算法循環t次下面的操作，隻有當t次的輸出結果都是素數時，該數就判定為素數，否則判定為合數

1、首先将n-1表示成(2^s)*d，其中s是非負整數，d是正奇數

2、在 [2, n-2] 内随機選擇整數a，計算x = a^d mod n, 如果x = 1或n-1 ，傳回 “是素數”，否則繼續

3、i = [1, s-1]循環如下操作

　　3.1，x = x^2 mod n

　　3.2、如果x = 1,傳回“是合數”

　　3.3、如果x = n-1,傳回“是素數”

4、傳回“是合數”

注意：其中3.2的判斷是基于以下定理：設x、y和n是整數，如果x^2=y^2 (mod n),但x ≠ ±y(mod n)，則(x-y)和n的公約數中有n的非平凡因子. 3.2中如果 x =1, 即 a^(2^i * d) = 1(mod n) 由上述定理(y = 1)可知：若式子（ttt） a^(2^(i-1) * d) ≠±1(mod n) 成立，則a^(2^(i-1) * d) - 1 和n有非平凡因子，進而可判斷n為合數。而上述式子（ttt）中a^(2^(i-1) * d)恰好是上一輪循環中的x，每次循環要繼續的條件包括x不等于1和 n-1（mod n 情況下n-1 和 -1等價。x=1或n-1時函數傳回了）。

算法3代碼如下，調用Miller_Rabin函數判斷，預設是測試40次

1 unsigned int exp_mod(unsigned int a, unsigned int b, unsigned int n)
 2 {   //a^b mod n
 3     unsigned long long ret = 1;
 4     unsigned long long a_copy = a;//防止大整數相乘溢出
 5     for (; b; b>>=1, a_copy = a_copy*a_copy%n)
 6         if (b&1)
 7             ret = ret*a_copy%n;
 8     return (unsigned int)ret;
 9 }
10 
11 //n是否通過以a為基的測試M-R測試
12 //傳回true時是素數，false是合數
13 bool witness(unsigned int a, unsigned int n)
14 {
15     unsigned int d = n-1,s = 0;
16     while((d&1) == 0)
17     {//n-1 轉換成 2^s * d
18         d >>= 1;
19         s++;
20     }
21     unsigned int x = exp_mod(a, d, n);//a^d mod n
22     if(x == 1 || x == n-1)return true;
23     for(unsigned int i = 1; i <= s-1; i++)
24     {
25         x = exp_mod(x, 2, n);
26         if(x == 1)return false;
27         else if(x == n-1)return true;
28     }
29     return false;
30 }
31 
32 bool Miller_Rabin(unsigned int n, int testTimes = 40)
33 {
34     if(n < 2)return false;
35     if(n == 2 || n == 3)return true;//後面要生成[2,n-2]的随機數，是以2,3提前判斷
36     if(n%2 == 0 || n%3 == 0)return false;
37     srand(time(NULL));
38     for(int i = 1; i <= testTimes; i++)
39     {
40         //産生[2,n-2]的随機數
41         unsigned int a = ((unsigned long long)rand())*(n-4)/RAND_MAX + 2;
42         if(witness(a, n) == false)
43             return false;
44     }
45     return true;
46 }

上述三個算法都可判斷unsigned int 可表示的整數範圍内的數，關于三個算法的效率：算法3效率是最高的，一般一個數在100~200微妙左右就能夠判斷，數越大，其優勢越明顯，總體而言速度是：算法3 > 算法2 > 算法1

下面提供一個計算某個區間[m,n]内的素數的函數，并把相應的素數寫進一個txt檔案，每行1000個，檔案的最後一行寫入該檔案内素數的個數，檔案以輸入區間命名，如果隻是想單純計算個數，把相應的檔案操作注釋即可，經過計算2^32-1内的素數總共203280221個，如果有需要可以下載下傳：不超過2^32-1的所有素數下載下傳本文位址

1 //生成[m,n]之間的素數，寫進檔案
 2 unsigned int GeneratePrime(unsigned int m, unsigned int n)
 3 {
 4     if(n<2 || m>n)return 0;
 5     unsigned int primeNum = 0;
 6     unsigned int tmp = sqrt(n);
 7     unsigned int *prime = new unsigned int[tmp+1];
 8     unsigned int k = findPrime(tmp, prime);//生成sqrt(n)以内的素數
 9     char filename[50];
10     sprintf(filename, "prime%u-%u.txt", m, n);
11     FILE *fp = fopen(filename, "w");
12     if(m < 2)m = 2;
13     for(long long i = m; i <= n; i++)
14     {
15         unsigned int tmp  = sqrt(i);
16         bool flag = true;
17         for(unsigned int j = 0; prime[j] <= tmp && j < k; j++)
18         {
19             if(i % prime[j] == 0)
20             {
21                 flag = false;
22                 break;
23             }
24         }
25         if(flag)
26         {
27             fprintf(fp, "%lld ", i);
28             primeNum++;
29             if(primeNum % 1000 == 0)fprintf(fp, "\n");
30         }
31     }
32     fprintf(fp, "\nprime number:%d", primeNum);
33     fclose(fp);
34     delete []prime;
35     return primeNum;//傳回素數的個數
36 }

本文所有代碼下載下傳（備用位址）

關于素數：求不超過n的素數，素數的判定(Miller Rabin 測試)

繼續閱讀

C語言第四章自述2第四章選擇結構程式設計

面試題:vector和map的差別，異同。空間分布，100萬資料存哪個比較合适。一、疊代器差別二、vector三、Map、Set四、vector_map 為什麼比map效率高五、如何選擇六、容器選擇原則七、效率對比

C++ 多線程用條件變量确定線程的執行順序而不是使用 sleep(1)

POJ 1284 Primitive Roots (歐拉函數&原根定理)

CQ V1.0分詞bates(基于雙數組tire樹)—應該是目前最快的中文分詞算法

成員函數初始化清單

2021-08-13c++——類之操作符重載

swmm與lisflood-fp源碼如何一起編譯 CMake指令

Windows下VS開發環境環境安裝工程項目設定關于Debug和Release的提示

一文看懂字元串的加減乘除

C++ 第十五周報告1--《冒泡法排序》

C++實作簡單順序表

C經典書籍筆記——C陷阱與缺陷②(文法陷阱之優先級)一、錯誤案列二、優先級規律

線性表之順序表的實作

C++判斷素數、求最大公約數代碼判斷一個數是否為素數求兩個數的最大公約數

SequoiaDB巨杉資料庫C++驅動概述