求素数大有学问
最近看了几道历年来找工作的笔试题目,很有几道是和素数相关的,本来也没有怎么上心,就觉得求素数么,不就弄个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 }
好了,大功告成,全部的代码奉上,希望各位多拍砖头。
![](https://images.cnblogs.com/OutliningIndicators/ContractedBlock.gif)
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 }