关于素数：求不超过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 }
【版权声明】转载请注明出处：http://www.cnblogs.com/TenosDoIt/p/3398112.html
