时间戳

首页 > 技术文章 > 求素数大有学问

求素数大有学问

havePassed 原文

最近看了几道历年来找工作的笔试题目，很有几道是和素数相关的，本来也没有怎么上心，就觉得求素数么，不就弄个for循环，判断到当前要判断的数的开方即可，可是linFen的博客让我看的是一愣一愣的，所以在此做个笔记。

定理1：如果n不是素数, 则n有满足1<d<=sqrt(n)的一个因子d。则有代码如下：

 1 //版本1
 2 //定理：如果n不是素数, 则n有满足1<d<=sqrt(n)的一个因子d.
 3 int isPrime(unsigned long long n)
 4 {
 5     if (n<2)
 6     {
 7         return 0;
 8     }
 9     if (n==2)
10     {
11         return 1;
12     }
13     for (unsigned long long i=3;i*i<=n;i++)
14     {
15         if(n%i==0) {
16             return 0;
17         }
18     }
19     return 1;
20 }

然后改进，去掉大于2的偶数的判断，可以节省一半的时间：

 1 //版本2，去掉偶数的判断
 2 int isPrime_without_even(unsigned long long n)
 3 {
 4     if (n<2)
 5     {
 6         return 0;
 7     }
 8     if (n==2)
 9     {
10         return 1;
11     }
12     if (n%2==0)
13     {
14         return 0;
15     }
16     for (unsigned long long i=3;i*i<=n;i+=2)
17     {
18         if(n%i==0) {
19             return 0;
20         }
21     }
22     return 1;
23 }

然后还有一个关于素数的定理如下：

定理: 如果n不是素数, 则n有满足1<d<=sqrt(n)的一个"素数"因子d.
证明:

I1. 如果n不是素数, 则n有满足1<d<=sqrt(n)的一个因子d.
I2. 如果d是素数, 则定理得证, 算法终止.
I3. 令n=d, 并转到步骤I1.

上述过程需要我们知道已有的素数的集合，则计算素数和改进的版本如下：

 1 //版本3方法补充
 2 //此方法甚好，可是有一点，需要构造一个素数保存序列,下列是构造素数序列函数
 3 //求前n个素数的素数序列。
 4 void makePrimes(unsigned long long primes[],unsigned int n)
 5 {
 6 //    FILE* zxh=fopen("zxh.txt","w");
 7 //    fprintf(zxh,"{2,3");
 8     primes[0]=2;
 9     primes[1]=3;
10     int flag=true;
11     unsigned long long i=5;
12     unsigned int pointer=2,j;
13     for(;pointer<n;i+=2)
14     {
15         flag=true;
16         for (j=0;primes[j]*primes[j]<=i;j++)
17         {
18             if (i%primes[j]==0)
19             {
20                 flag=false;
21             }
22         }
23         if (flag==true)
24         {
25             primes[pointer++]=i;
26 //            fprintf(zxh,",%lld",i);
27         }
28     }
29 //    fprintf(zxh,"}");
30 //    fclose(zxh);
31 }
32 
33 //版本3
34 int isPrime_with_primes(unsigned long long primes[],unsigned long long n)
35 {
36     if (n<2)
37     {
38         return 0;
39     }
40     for (unsigned long long i=0;primes[i]*primes[i]<=n;i++)
41     {
42         if (n%primes[i]==0)
43         {
44             return 0;
45         }
46     }
47     return 1;
48 }

其实如果我们经常求做对素数的操作的时候，我们还可以把以前生成的素数都保存下来，以便于以后使用。正如上面makePrimes函数中我注释掉的，我们可以把计算好的素数保存下来到一个文件中，以后每次使用的时候只需要拷贝即可。

根据素数定理：

x/(lnx-0.5)<π(x)<x/(lnx-1.5)，其中π(x)表示小于x的素数的个数。

左边不等式对于x>=67成立,右边不等式对于x>√e3≈4.48169...成立.

  

则我们可以事先算好小于unsigned long long内的所有的素数即可。

假设我们要判断unsigned long long (2^64)以内的素数，我们就需要2^32内的素数即可判断，即有：

Max unsigned long (MUL）= 2^32=4294967296

π(MUL)<MUL/(ln(MUL)-1.5)=207679878.59807090275597192690663≈207679879

即我们只要事先计算好207679879个素数的就可以很容易的判断2^64范围之内的所有的数是否为素数。

假设我们已经有计算好的素数的数组Primes，则当小于我们已经计算好的素数的最大值时候，只需要在素数数组中采用二分查找即可，当大于的时候需要上述方法来判断。代码如下：

 1 //  二分查找，看数组Primes里面是否存在一个elem
 2 //  若存在，返回索引值
 3 //  若不存在，返回-1
 4 int binary_search_prime(unsigned long long Primes[],unsigned long long elem,unsigned int low,unsigned int high)
 5 {
 6     if (low>PRIME_NUM-1||high<low||high>PRIME_NUM-1)
 7     {
 8         return 0;
 9     }
10     int mid=0;
11     while(low<=high){
12         mid=(low+high)/2;
13         if (elem==Primes[mid])
14         {
15             return mid;
16         }
17         else if(elem<Primes[mid])
18         {
19             high=mid-1;
20         }
21         else
22         {
23             low=mid+1;
24         }
25     }
26     return -1;
27 }
28 
29 //  总体来讲判断一个64位数是否是素数，
30 //  当小于我们已经计算好的素数的最大值时候，只需要在素数数组中查找即可，当大于的时候需要上述方法来判断。
31 //  假设我们所求的素数是在已经缓存素数的范围内，则在判断是否是素数的时候只需要二分查找即可
32 int is_prime_with_precompute(unsigned long long Primes[],unsigned long long elem)
33 {
34     if (elem<=Primes[PRIME_NUM-1])
35     {
36         return (-1!=binary_search_prime(Primes,elem,0,PRIME_NUM-1));
37     }
38     else
39     {
40         unsigned int j=0;
41         for (j=0;Primes[j]*Primes[j]<=elem;j++)
42         {
43             if (elem%Primes[j]==0)
44             {
45                 return 0;
46             }
47         }
48     }
49     return 1;
50 }

上面的这种方法已经省去了很多的时间，可以大大帮助判断64位整数是否是素数，但是根据素数定理：

x/(lnx-0.5)<π(x)<x/(lnx-1.5)，其中π(x)表示小于x的素数的个数。

左边不等式对于x>=67成立,右边不等式对于x>√e3≈4.48169...成立.

也就是说，我们已知了x之后，我们就可以判断小于x的素数大概有多少个，则我们在查找素数表的时候就不用在[0,PRIME_NUM-1]范围内查找了，新的查找范围是[x/(lnx-0.5)-1,x/(lnx-1.5)]，则is_prime_with_precompute算法可以改写如下：

 1 int is_prime_with_precompute_optimize(unsigned long long Primes[],unsigned long long elem)
 2 {
 3     if (elem>67&&elem<=Primes[PRIME_NUM-1])
 4     {
 5         unsigned int 
 6         low=(int)(elem/(log((double)elem)-0.5))-1,
 7         high=(int)(elem/(log((double)elem)-1.5));
 8 
 9         return (-1!=binary_search_prime(Primes,elem,low,high));
10     }
11     else
12     {
13         unsigned int j=0;
14         for (j=0;Primes[j]*Primes[j]<=elem;j++)
15         {
16             if (elem%Primes[j]==0)
17             {
18                 return 0;
19             }
20         }
21     }
22     return 1;
23 }

好了，大功告成，全部的代码奉上，希望各位多拍砖头。

  1 #include <math.h>
  2 #define PRIME_NUM 2076
  3 
  4 //版本1
  5 //定理：如果n不是素数, 则n有满足1<d<=sqrt(n)的一个因子d.
  6 int isPrime(unsigned long long n)
  7 {
  8     if (n<2)
  9     {
 10         return 0;
 11     }
 12     if (n==2)
 13     {
 14         return 1;
 15     }
 16     for (unsigned long long i=3;i*i<=n;i++)
 17     {
 18         if(n%i==0) {
 19             return 0;
 20         }
 21     }
 22     return 1;
 23 }
 24 //版本2，去掉偶数的判断
 25 int isPrime_without_even(unsigned long long n)
 26 {
 27     if (n<2)
 28     {
 29         return 0;
 30     }
 31     if (n==2)
 32     {
 33         return 1;
 34     }
 35     if (n%2==0)
 36     {
 37         return 0;
 38     }
 39     for (unsigned long long i=3;i*i<=n;i+=2)
 40     {
 41         if(n%i==0) {
 42             return 0;
 43         }
 44     }
 45     return 1;
 46 }
 47 
 48 //版本3方法补充
 49 //此方法甚好，可是有一点，需要构造一个素数保存序列,下列是构造素数序列函数
 50 //求前n个素数的素数序列。
 51 void makePrimes(unsigned long long primes[],unsigned int n)
 52 {
 53     FILE* zxh=fopen("zxh.txt","w");
 54     fprintf(zxh,"{2,3");
 55     primes[0]=2;
 56     primes[1]=3;
 57     int flag=true;
 58     unsigned long long i=5;
 59     unsigned int pointer=2,j;
 60     for(;pointer<n;i+=2)
 61     {
 62         flag=true;
 63         for (j=0;primes[j]*primes[j]<=i;j++)
 64         {
 65             if (i%primes[j]==0)
 66             {
 67                 flag=false;
 68             }
 69         }
 70         if (flag==true)
 71         {
 72             primes[pointer++]=i;
 73             fprintf(zxh,",%lld",i);
 74         }
 75     }
 76     fprintf(zxh,"}");
 77     fclose(zxh);
 78 }
 79 
 80 //版本3
 81 //定理: 如果n不是素数, 则n有满足1<d<=sqrt(n)的一个"素数"因子d.
 82 //证明: I1. 如果n不是素数, 则n有满足1<d<=sqrt(n)的一个因子d.
 83 //        I2. 如果d是素数, 则定理得证, 算法终止.
 84 //        I3. 令n=d, 并转到步骤I1.
 85 int isPrime_with_primes(unsigned long long primes[],unsigned long long n)
 86 {
 87     if (n<2)
 88     {
 89         return 0;
 90     }
 91     for (unsigned long long i=0;primes[i]*primes[i]<=n;i++)
 92     {
 93         if (n%primes[i]==0)
 94         {
 95             return 0;
 96         }
 97     }
 98     return 1;
 99 }
100 
101 //  根据素数定理：
102 //  这个公式的新近改进如下:
103 //  x/(lnx-0.5)<π(x)<x/(lnx-1.5)
104 //    左边不等式对于x>=67成立,右边不等式对于x>√e3≈4.48169...成立.
105 //  
106 //  则我们可以事先算好小于unsigned long long内的所有的素数即可。
107 //  假设我们要判断unsigned long long 2^64以内的素数，我们就需要2^32内的素数即可判断，即有：
108 
109 //  Max unsigned long (MUL）=2^32=4294967296
110 //  π(MUL)<MUL/(ln(MUL)-1.5)=207679878.59807090275597192690663≈207679879
111 //  即我们只要事先计算好207679879个素数的话 就可以很容易的判断2^64范围之内的所有素数
112 
113 
114 //  二分查找，看数组Primes里面是否存在一个elem
115 //  若存在，返回索引值
116 //  若不存在，返回-1
117 int binary_search_prime(unsigned long long Primes[],unsigned long long elem,unsigned int low,unsigned int high)
118 {
119     if (low>PRIME_NUM-1||high<low||high>PRIME_NUM-1)
120     {
121         return 0;
122     }
123     int mid=0;
124     while(low<=high){
125         mid=(low+high)/2;
126         if (elem==Primes[mid])
127         {
128             return mid;
129         }
130         else if(elem<Primes[mid])
131         {
132             high=mid-1;
133         }
134         else
135         {
136             low=mid+1;
137         }
138     }
139     return -1;
140 }
141 
142 //  总体来讲判断一个64位数是否是素数，
143 //  当小于我们已经计算好的素数的最大值时候，只需要在素数数组中查找即可，当大于的时候需要上述方法来判断。
144 //  假设我们所求的素数是在已经缓存素数的范围内，则在判断是否是素数的时候只需要二分查找即可
145 int is_prime_with_precompute(unsigned long long Primes[],unsigned long long elem)
146 {
147     if (elem<=Primes[PRIME_NUM-1])
148     {
149         return (-1!=binary_search_prime(Primes,elem,0,PRIME_NUM-1));
150     }
151     else
152     {
153         unsigned int j=0;
154         for (j=0;Primes[j]*Primes[j]<=elem;j++)
155         {
156             if (elem%Primes[j]==0)
157             {
158                 return 0;
159             }
160         }
161     }
162     return 1;
163 }
164 
165 //  上面的这种方法已经省去了很多的时间，可以大大帮助判断64位整数是否是素数
166 //  但是我们在上面的素数定理中也说了，有：
167 //  x/(lnx-0.5)<π(x)<x/(lnx-1.5)
168 //    左边不等式对于x>=67成立,右边不等式对于x>√e3≈4.48169...成立.
169 //  也就是说，我们已知了x之后，我们就可以判断小于x的素数大概有多少个，则我们在查找素数表的时候就不用在[0,PRIME_NUM-1]范围内查找了
170 //  新的查找范围是[x/(lnx-0.5)-1,x/(lnx-1.5)]，则上述算法可以改写如下：
171 int is_prime_with_precompute_optimize(unsigned long long Primes[],unsigned long long elem)
172 {
173     if (elem>67&&elem<=Primes[PRIME_NUM-1])
174     {
175         unsigned int 
176         low=(int)(elem/(log((double)elem)-0.5))-1,
177         high=(int)(elem/(log((double)elem)-1.5));
178 
179         return (-1!=binary_search_prime(Primes,elem,low,high));
180     }
181     else
182     {
183         unsigned int j=0;
184         for (j=0;Primes[j]*Primes[j]<=elem;j++)
185         {
186             if (elem%Primes[j]==0)
187             {
188                 return 0;
189             }
190         }
191     }
192     return 1;
193 }
View Code

推荐阅读