本博客rss订阅地址: http://feed.cnblogs.com/blog/u/147990/rss

关于素数：求不超过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






            

            
posted @ 
2013-10-30 22:41 
tenos 
阅读(3962) 
评论(2) 
编辑 
收藏 
举报


        

	    
	    
    









    
    

    
    

    

        

    

        

            
                
                
            
        

    

        

            
                
            
        

    

    

    



	

	

	
	

	

	

	





    

        
 本博客rss订阅地址: http://feed.cnblogs.com/blog/u/147990/rss 



公益页面-寻找遗失儿童