质数计数

前言

很久以前,我出了区间质数,当然求 之间的质数等于π(r)π(l1)\pi(r)-\pi(l-1),其中π(n)\pi(n)为质数计数函数,表示 中质数的个数。

我最初给出的方法是分段打表+Miller Rabin素性判定,可以解决n109n\le10^9规模,但是更大的规模无法打表,因为内存不够……最近看到洛谷P3912后,zyy提供了LibreOJ #6235,求π(n),n1011\pi(n),n\le10^{11}。对此我想大致了解一下质数计数算法。

筛法

计算π(n)\pi(n)最简单的想法就是筛法,这里同时来比较一下埃氏筛和欧拉筛的实际效率,以及bitset优化的效率。与算法竞赛不同的是,代码避免使用静态数组,而是根据输入规模动态分配内存。因此与静态数组的性能有一些差异。

埃氏筛

原始版本(era.cpp)

n\sqrt n内的质数的倍数筛掉,时间复杂度O(nloglogn)O(n\log\log n),空间复杂度O(n)O(n)。一个小的优化是直接忽略偶数(除了2)。这里假设n2n\ge2

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
#include <iostream>
#include <algorithm>
#include <ctime>
using namespace std;
typedef unsigned int_t;
int main()
{
int_t n;
cin >> n;
clock_t start = clock();
bool *p = new bool[n + 1];
fill(p + 2, p + n + 1, true);
for (int_t i = 3; i * i <= n; i += 2)
if (p[i])
for (int_t j = i * i; j <= n; j += i + i)
p[j] = false;
int_t cnt = 1;
for (int_t i = 3; i <= n; i += 2)
cnt += p[i];
cout << cnt << endl;
delete[] p;
cout << clock() - start << endl;
return 0;
}

bitset(era_bitset.cpp)

值得注意的是,std::bitset是静态的模板,因此我实际上用了std::vector\。时间复杂度不变,但是常数小,因为空间局部性吧。空间复杂度O(n8)O(\frac n8)

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
#include <iostream>
#include <vector>
#include <ctime>
using namespace std;
typedef unsigned int_t;
int main()
{
int_t n;
cin >> n;
clock_t start = clock();
vector<bool> p(n + 1, true);
for (int_t i = 3; i * i <= n; i += 2)
if (p[i])
for (int_t j = i * i; j <= n; j += i + i)
p[j] = false;
int_t cnt = 1;
for (int_t i = 3; i <= n; i += 2)
cnt += p[i];
cout << cnt << endl;
cout << clock() - start << endl;
return 0;
}

分段筛(era_segmented.cpp)

我最早是在某篇论文上看到的,但是http://primesieve.org/segmented_sieve.html的介绍和代码很容易理解。先生成n\sqrt n内的质数表,然后分段筛,段长度接近L1数据缓存。能很大程度提高效率,因为增加了缓存命中率。空间复杂度O(π(n)+B)O(\pi(\sqrt n)+B),BB为段长度。

另外,这可以计算较小的区间内的质数,如 。这样取段大小10810^8,应该可以通过LOJ那题。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
#include <iostream>
#include <vector>
#include <algorithm>
#include <cmath>
#include <ctime>
using namespace std;
typedef unsigned int_t;
const int_t seg_size = 32768;
int main()
{
int_t n;
cin >> n;
clock_t start = clock();
int_t m = sqrt(n);
bool *p = new bool[m + 1];
fill(p + 2, p + m + 1, true);
for (int_t i = 2; i * i <= m; i++)
if (p[i])
for (int_t j = i * i; j <= m; j += i)
p[j] = false;
bool *buf = new bool[seg_size];
vector<int_t> primes;
vector<int_t> nxt;
int_t cnt = 1;
for (int_t low = 0, s = 3, now = 3; low <= n; low += seg_size)
{
fill(buf, buf + seg_size, true);
int_t high = min(low + seg_size - 1, n);
for (; s * s <= high; s += 2)
if (p[s])
{
primes.push_back(s);
nxt.push_back(s * s - low);
}
for (size_t i = 0; i < primes.size(); i++)
{
int_t j = nxt[i];
for (int_t k = primes[i] * 2; j < seg_size; j += k)
buf[j] = false;
nxt[i] = j - seg_size;
}
for (; now <= high; now += 2)
cnt += buf[now - low];
}
cout << cnt << endl;
delete[] p;
delete[] buf;
cout << clock() - start << endl;
return 0;
}

欧拉筛

原始版本(euler.cpp)

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
#include <iostream>
#include <vector>
#include <algorithm>
#include <ctime>
using namespace std;
typedef unsigned int_t;
int main()
{
int_t n;
cin >> n;
clock_t start = clock();
bool *p = new bool[n + 1];
fill(p + 2, p + n + 1, true);
vector<int_t> primes;
for (int_t i = 3; i <= n; i += 2)
{
if (p[i])
primes.push_back(i);
for (int_t j : primes)
{
if (i * j > n)
break;
p[i * j] = false;
if (i % j == 0)
break;
}
}
cout << primes.size() + 1 << endl;
delete[] p;
cout << clock() - start << endl;
return 0;
}

bitset(euler_bitset.cpp)

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
#include <iostream>
#include <vector>
#include <ctime>
using namespace std;
typedef unsigned int_t;
int main()
{
int_t n;
cin >> n;
clock_t start = clock();
vector<bool> p(n + 1, true);
vector<int_t> primes;
for (int_t i = 3; i <= n; i += 2)
{
if (p[i])
primes.push_back(i);
for (int_t j : primes)
{
if (i * j > n)
break;
p[i * j] = false;
if (i % j == 0)
break;
}
}
cout << primes.size() + 1 << endl;
cout << clock() - start << endl;
return 0;
}

效率

开O2优化,用clock测量,取最小值,单位:毫秒。

nn π(n)\pi(n) era era_bitset era_seg euler euler_bitset
10710^7 664,579664,579 7979 2020 1212 5353 4444
10810^8 5,761,4555,761,455 1,0141,014 538538 9191 621621 492492
10910^9 50,847,53450,847,534 11,75211,752 7,1317,131 1,1941,194 6,5336,533 5,3095,309