与素数有关的知识---模版

转自:素数算法

一、引言

在平时做题目或者进行预算的时候，素数的出现次数总是十分频繁。今天我们就来一点一点的说一说关于素数的一些算法。

二、朴素判断素数算法

就判断素数而言，事实上是非常简单的了。根据定义，判断一个整数n是否是素数，只需要去判断在整数区间[2, n-1]之内，是否具有某个数m，使得n % m == 0。代码可以这么写：

[cpp] view plain copy

 

1. int isPrime(int n) {  
2.     int i;  
3.     for (i = 2; i < n; ++i) {  
4.         if (n % i == 0) return 0;  
5.     }  
6.     return 1;  
7. }  

事实上，这个算法是O(n)的，感觉是很快了，但是依旧无法满足需求。所以有一个算法是O(sqrt(n))的算法。代码可以这么写：

[cpp] view plain copy

 

1. int isPrime(long long n) {  
2.     long long i;  
3.     for (i = 2; i * i <= n; ++i) {  
4.         if (n % i == 0) return 0;  
5.     }  
6.     return 1;  
7. }  

原理很巧妙，也仅仅是把代码的i < n变成了i * i <= n而已，但是优化却是极其高的。可能你会注意到，在上一份代码里面，我定义的n为int类型，而后面一份代码，我却定义成了long long，这是因为在1s内，上一份代码能判断出来的数量级为1e8，而后面一份代码判断出来的数量级却几乎可以达到1e16。而原因仅仅是因为a * b = c的话一旦a是c的约数，那么b也是，如果a不是，那么b也不是。

这个方法也可以称作试除法。

三、Miller_Rabin素性测试

尽管上面的O(sqrt(n))的算法已经非常优秀了，但是面对更高数量级的“大数”却会显得力不从心。而这个时候，Miller_Rabin的优越性就显示出来了。

Miller_Rabin的理论基础来源于费马小定理。值得一提的是费马小定理是数论四大定理之一。

费马小定理：n是一个奇素数，a是任何整数(1≤ a≤n-1) ，则 a^(n-1)≡1(mod n)

要测试n是否是一个素数，首先将n-1分解为(2^s) * d，在每次测试开始时，先随机选择一个介于[1, N - 1]的整数a，之后如果对于所有的r∈[0, s - 1]，若a^d mod N ≠ 1且a^((2^r) * d) mod N ≠ -1，那么n就是一个合数，否则n有3/4的几率是素数。

这也是为什么说这个算法只是素性测试了，他不能完全保证结果正确，但是当测试次数大约为20的时候，基本可以确保正确率达到97%以上。

它的代码可以这么写：

[cpp] view plain copy

 

1. const int S = 20;  
2. LL mod_mul(LL a, LL b, LL n) {  
3.     LL res = 0;  
4.     while (b) {  
5.         if (b & 1) res = (res + a) % n;  
6.         a = (a + a) % n;  
7.         b >>= 1;  
8.     }  
9.     return res;  
10. }  
11. LL mod_exp(LL a, LL b, LL n) {  
12.     LL res = 1;  
13.     while (b) {  
14.         if (b & 1) res = mod_mul(res, a, n);  
15.         a = mod_mul(a, a, n);  
16.         b >>= 1;  
17.     }  
18.     return res;  
19. }  
20. bool miller_rabin(LL n) {  
21.     if (n == 2 || n == 3 || n == 5 || n == 7 || n == 11) return true;  
22.     if (n == 1 || !(n % 2) || !(n % 3) || !(n % 5) || !(n % 7) || !(n % 11)) return false;  
23.   
24.     LL x, pre, u = n - 1, k = 0;  
25.   
26.     while (!(u & 1)) {  
27.         ++k;  
28.         u >>= 1;  
29.     }  
30.   
31.     srand((LL)time(NULL));  
32.     for (int i = 0; i < S; ++i) {                     //进行S次测试，S越大，结果越准确  
33.         x = rand() % (n - 2) + 2;                   //在[2, n)中取随机数  
34.         if (x % n == 0) continue;  
35.   
36.         x = mod_exp(x, u, n);                       //计算x^u % n  
37.         pre = x;  
38.         for (int j = 0; j < k; ++j) {  
39.             x = mod_mul(x, x, n);  
40.             if (x == 1 && pre != 1 && pre != n - 1)  
41.                 return false;  
42.             pre = x;  
43.         }  
44.         if (x != 1) return false;  
45.     }  
46.     return true;  
47. }  

四、筛选法

上面介绍的一些素数判断的算法，在某些问题是基本可以适用的了。但是对于另外一类问题却十分尴尬。比如问你整数区间[1, n]内的素数个数是多少。这个时候如果一个个枚举，一个个判断，对于朴素方法来说，最优也是O(nsqrt(n))，即使是Miller_Rabin算法也无法在O(n)的时间内得到结果。于是就有了埃氏筛选法（埃拉托斯特尼筛法）。

对于筛选整数n以内的素数，算法是这么描述的：先把素数2的倍数全部删除，剩下的数第一个为3，再把素数3的倍数全部删除，剩下的第一个数为5，再把素数5的倍数全部删除······直到把n以内最后一个素数的倍数删除完毕，得到n以内的所有素数。代码可以这么写：

[cpp] view plain copy

 

1. const int maxn = 1e7 + 5;  
2. int pri[maxn];  
3. void getPrime(int n) {  
4.     for (int i = 0; i <= n; ++i) pri[i] = i;  
5.     pri[1] = 0;  
6.     for (int i = 2; i <= n; ++i) {  
7.         if (!pri[i]) continue;  
8.         pri[++pri[0]] = i;  
9.         for (int j = 2; i * j <= n && j < n; ++j) {  
10.             pri[i * j] = 0;  
11.         }  
12.     }  
13. }  

而事实上，观察可以发现，上面的这种写法有很多次重复计算，这样显然无法做到线性筛选，而另外一种写法却可以得到线性筛选，达到时间复杂度为O(n)，代码可以这么写：

[cpp] view plain copy

 

1. const int MAXN = 1e7 + 5;  
2. void getPrime(){  
3.     memset(prime, 0, sizeof(prime));  
4.     for (int i = 2;i <= MAXN;i++) {  
5.         if (!prime[i])prime[++prime[0]] = i;  
6.         for (int j = 1;j <= prime[0] && prime[j] <= MAXN / i;j++) {  
7.             prime[prime[j] * i] = 1;  
8.             if (i%prime[j] == 0) break;  
9.         }  
10.     }  
11. }  

来自kuangbin的模板。

五、容斥原理

从上面的代码可以发现，显然这种筛法只能应付达到1e7这种数量级的运算，即使是线性的筛选法，也无法满足，因为在ACM竞赛中，1e8的内存是极有可能获得Memery Limit Exceed的。

于是可以考虑容斥原理。

以AHUOJ 557为例，1e8的情况是筛选法完全无法满足的，但是还是考虑a * b = c的情况，1e8只需要考虑10000以内的素数p[10000]，然后每次先减去n / p[i]，再加上n / (p[i] * p[j])再减去n / (p[i] * p[j] * p[k])以此类推...于是就可以得到正确结果了。

代码如下：

[cpp] view plain copy

 

1. #include <cmath>  
2. #include <cstdio>  
3. using namespace std;  
4.    
5. const int maxn = 10005;  
6. int sqrn, n, ans = 0;  
7. bool vis[maxn];  
8. int pri[1500] = {0};  
9. void init(){  
10.     vis[1] = true;  
11.     int k = 0;  
12.     for(int i = 2; i < maxn; i++){  
13.         if(!vis[i]) pri[k++] = i;  
14.         for(int j = 0; j < k && pri[j] * i < maxn; j++){  
15.             vis[pri[j] * i] = true;  
16.             if(i % pri[j] == 0) break;  
17.         }  
18.     }  
19. }  
20. void dfs(int num, int res, int index){  
21.     for(int i = index; pri[i] <= sqrn; i++){  
22.         if(1LL * res * pri[i] > n){  
23.             return;  
24.         }  
25.         dfs(num + 1, res * pri[i], i+1);  
26.         if(num % 2 == 1){  
27.             ans -= n / (res * pri[i]);  
28.         }else{  
29.             ans += n / (res * pri[i]);  
30.         }  
31.    
32.         if(num == 1) ans++;  
33.     }  
34. }  
35. int main(){  
36.     init();  
37.     while(~scanf("%d",&n) && n){  
38.         ans = n;  
39.         sqrn = sqrt((double)n);  
40.         dfs(1,1,0);  
41.         printf("%d\n",ans-1);  
42.     }  
43.     return 0;  
44. }  

六、Meissel-Lehmer算法

最后介绍的这个算法可以说是黑科技级别的，它能够在3s内求出1e11之内的素数个数。据说还有在300ms内求出1e11的个数的。可以参考wiki里面原理。然后给出来自Codeforces 665F题目里面的代码。

[cpp] view plain copy

 

1. #define MAXN 100    // pre-calc max n for phi(m, n)  
2. #define MAXM 10010 // pre-calc max m for phi(m, n)  
3. #define MAXP 40000 // max primes counter  
4. #define MAX 400010    // max prime  
5. #define setbit(ar, i) (((ar[(i) >> 6]) |= (1 << (((i) >> 1) & 31))))   
6. #define chkbit(ar, i) (((ar[(i) >> 6]) & (1 << (((i) >> 1) & 31))))  
7. #define isprime(x) (( (x) && ((x)&1) && (!chkbit(ar, (x)))) || ((x) == 2))  
8.   
9. namespace pcf {  
10.     long long dp[MAXN][MAXM];  
11.     unsigned int ar[(MAX >> 6) + 5] = { 0 };  
12.     int len = 0, primes[MAXP], counter[MAX];  
13.   
14.     void Sieve() {  
15.         setbit(ar, 0), setbit(ar, 1);  
16.         for (int i = 3; (i * i) < MAX; i++, i++) {  
17.             if (!chkbit(ar, i)) {  
18.                 int k = i << 1;  
19.                 for (int j = (i * i); j < MAX; j += k) setbit(ar, j);  
20.             }  
21.         }  
22.   
23.         for (int i = 1; i < MAX; i++) {  
24.             counter[i] = counter[i - 1];  
25.             if (isprime(i)) primes[len++] = i, counter[i]++;  
26.         }  
27.     }  
28.   
29.     void init() {  
30.         Sieve();  
31.         for (int n = 0; n < MAXN; n++) {  
32.             for (int m = 0; m < MAXM; m++) {  
33.                 if (!n) dp[n][m] = m;  
34.                 else dp[n][m] = dp[n - 1][m] - dp[n - 1][m / primes[n - 1]];  
35.             }  
36.         }  
37.     }  
38.   
39.     long long phi(long long m, int n) {  
40.         if (n == 0) return m;  
41.         if (primes[n - 1] >= m) return 1;  
42.         if (m < MAXM && n < MAXN) return dp[n][m];  
43.         return phi(m, n - 1) - phi(m / primes[n - 1], n - 1);  
44.     }  
45.   
46.     long long Lehmer(long long m) {  
47.         if (m < MAX) return counter[m];  
48.   
49.         long long w, res = 0;  
50.         int i, a, s, c, x, y;  
51.         s = sqrt(0.9 + m), y = c = cbrt(0.9 + m);  
52.         a = counter[y], res = phi(m, a) + a - 1;  
53.         for (i = a; primes[i] <= s; i++) res = res - Lehmer(m / primes[i]) + Lehmer(primes[i]) - 1;  
54.         return res;  
55.     }  
56. }