libreoj 6235 区间素数个数 min25筛

缘起

学习min_25 筛, 遂找到此模板题 libreoj 6235. 区间素数个数

分析

题是中文的

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
给定一个正整数n， 求[1,n]中的素数的个数.

【输入】
一个正整数n（2<=n<=1e11）

【输出】
[1,n]范围内正整数的个数

【样例输入】
10

【样例输出】
4

【限制】
时间: 1s
内存 512M

首先，先说一种完全基于线性筛的方法——分段打表. 即将 [1,1e11]分成5000个2000w的区间. 然后就可以分段统计，每段上的质数的个数可以使用线性筛搞定（这种区间长度有上限的做法在【1】中已经讲过了）.事实上，打比赛的当天也有人这么a掉的，只是吃相太难看了. 下面给出基于dp的方法（素数个数用dp? 一开始我也觉得脑洞略大~）.

其实思想还是基于线性筛法. 你想想看线性筛我们做了什么? 每次都使用一个质数去筛掉它的倍数（即设置它的倍数为合数）.那么使用质数{p1,.p2,….,pi}筛不掉[1,n]中哪个合数呢（即使用p1,…,pi筛过之后, 剩下某个数——虽然它是合数，但依旧被认定为质数）? 答案是筛不掉最小质因子>pi的合数! 所以如果我们使用的pi足够大, 大到[1,…,,n]中不存在最小素因子>pi的时候, 则我们只需要使用{p1,.,,,.pi}去筛，则剩下的就必然全部是素数了. 而我们只需要令p_cnt为 <=sqrt(n)的最大质数，然后使用{p1,…,p_cnt}去筛[1,…,n]. 剩下的就是[2,n]范围内的素数. 所以我们自然的, 可以令

1
2
g(x,i) 为使用 {p1,...,pi}去筛 [2,..,x]剩下的数的个数（这些数显然由两部分组成,第一部分就是质数, 第
二部分是最小质因子>pi的数）. 其中 2<=x<=n

则g(n,cnt)就是答案. 于是, dp的味道就出来了. 自然要考虑

1
g(x,i)是如何通过g(x,i-1)得到的. 2<=x<=n

则我们有（基于该dp公式的算法就称为min_25筛，由一位名为min_25的菊苣发明，事实上，严格讲，只是min_25算法的一部分.）

1
2
3
4
g(x,i-1) - g(x,i) = g(x/pi,i-1)-(i-1), if pi^2<=x<=n
g(x,i-1) - g(x,i) = 0, if 2<=x<pi^2
i=1,...,cnt
g(x, 0)=x-1, 2<=x<=n

第一行是因为等式的左边=[1,n]中的素数个数+最小质因子>p_{i-1}的合数个数-[1,n]中的素数个数-最小质因子>pi的合数的个数=最小质因子=pi的合数的个数

而等式的右边为什么要减去i-1呢? 因为g(x/pi,i-1)中包含[1,..,x/pi]中的所有素数, 其中 p1,…,p_{i-1}是不能被算进去的（注意, x/pi因为 pi不超过sqrt(x),所以 x/pi 一定大于p_{i-1}）. 为啥不能算进去? 因为乘以pi还原回去的时候是p_{i-1}*pi, 但是这个数的最小质因子!=pi的. 所以要去掉.

第二行是因为使用{p1,…,pi}和{p1,…,p_{i-1}}进行筛法的时候没区别，因为使用pi已经不能筛掉任何一个数了（因为pi>sqrt(x), 所以[2,x]中不存在最小质因子>pi的数）

第四行根据定义是成立的. 其中p0=1——则筛法不能筛掉[2,x]中任何一个数——即认为[2,x]中所有数都是质数.

是不是觉得dp在手，天下我有? 一开始我也这么认为，天真了~ 先给出如下超时解法.

基于非记忆化递归（TLE）

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
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
//#include "stdafx.h"
#include <stdio.h>
#include <math.h>
#include <string.h>
//#define LOCAL
#define LL long long

const int maxn = 1e6+5;

int primes[maxn], num;
bool isprime[maxn];

void sieve()
{
	memset(isprime, true, sizeof(isprime));
	for (int i = 2; i<maxn; i++)
	{
		if (isprime[i])
		{
			for (int j = (i<<1); j<maxn; j+=i)
			{
				isprime[j] = false;
			}
		}
	}
	for (int i = 2; i<maxn; i++)
	{
		if (isprime[i])
		{
			primes[++num] = i;
		}
	}
}

LL g(LL n, int lim) // 非记忆化递归是超时的根源
{
	if (!lim)
	{
		return n-1;
	}
	if (primes[lim]*primes[lim]>n)
	{
		return g(n, lim-1);
	}
	return g(n, lim-1)-(g(n/primes[lim], lim-1)-(lim-1));
}

int getlim(int x) // 找 primes 中<=x的最大素数， 使用二分查找
{
	int lo = 1, hi = num;
	while(lo<hi)
	{
		int mi = (lo+hi)>>1;
		x<primes[mi]?hi=mi:lo=mi+1;
	}
	return --lo;
}

int main(){
#ifdef LOCAL
	freopen("d:\\data.in", "r", stdin);
#endif
	sieve();
	LL n;
	scanf("%lld",&n);
	int sq = sqrt(n*1.0);
	printf("%lld", g(n, getlim(sq)));
	return 0;
}

结果是吃T（只能骗小数据分, 大的一定超时）

编号 题目 状态 分数 总时间 内存 代码 / 答案文件 提交者 提交时间
#569661 #6235. 区间素数个数 Time Limit Exceeded 55 9787 ms 2944 K C++ / 1.0 K yfs 2019-08-13 13:48:56

不用说也知道——递归可能导致状态被大量重复计算. 但是我们又没办法打表进行记忆化. 因为空间复杂度不够（1000亿~）.所以我们必须要进行空间复杂度的优化. 首先，dp公式是二维的，滚动数组显然可以滚掉一维. 但是仍然不够——因为一维仍旧是1e11的——依旧是无法承受之痛. 怎么办? 这时候我们要扪心自问——g(x,i)（2<=x<=n，n可达1e11）都是对于求 g(n,cnt)所必须的么? 其实不然， 我们只需找到一个[2,n]的子集S（它其实可以视作dp的状态集）.然后可以证明 S在dp公式逐步转移过程中跳不出这个范围就好了.

我们取

1
S={n/l,n/2,...,n/lim,lim,lim-1,...,1}

其中, 如果n是完全平方数的话, n/lim和lim只保留一个（这是显然的，不然n/lim=lim都要进行一次dp，那岂不是一轮dp中n/lim=lim发生了2次dp? 这显然不是我们的本意）.

只需要证明

1
S/pj 是 S 的子集即可, 1<=j<=cnt

这是不难使用数学归纳法证明的. 就此略过.

而|S|=O(sqrt(n)) 这就是我们可以承受的了. 综上所述，其实我们就是在 S上而不是 {1,…,n}上进行dp. 而且我们需要的答案正好包含在S进行dp的结果中. 所以给出如下ac算法

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
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
//#include "stdafx.h"

#include <stdio.h>
#include <math.h>
#include <string.h>

//#define LOCAL

typedef long long LL;

LL n, lim, top;

const LL maxn = 2*1e6+5, maxnn = 1e6+5;

bool isprime[maxn], flag;

LL g[maxn], h[maxn], primes[maxn];

void sieve() // 线性筛出 100万以内的质数
{
	memset(isprime, true, sizeof(isprime));
	for (LL i = 2; i<maxnn; i++)
	{
		if (isprime[i])
		{
			for (LL j = (i<<1); j<maxnn; j+=i)
			{
				isprime[j] = false;
			}
		}
	}
	for (LL i = 2; i<maxnn; i++)
	{
		if (isprime[i])
		{
			primes[++top] = i;
		}
	}
}

LL sz(LL x)
{
	if (flag) // 如果n是完全平方数
	{
		if (x>=lim)
		{
			return n/x;
		}
		return (lim<<1)-x;
	}
	if (x>lim)
	{
		return n/x;
	}
	return (lim<<1)+1-x;
}

void init() // N/1,N/2,...,N/lim,lim,...,2,1 （如果N是完全平方数的话, 则只保留N/lim和lim中的一个）
{
	lim = sqrt(1.0*n);
	for (LL i = 1; i<=lim; i++)
	{
		h[i] = n/i;
		// 使用 map 存储 <n/i, i> 这个键值对，但是使用map这个stl会超时,所以使用sz来优化（sz就是用来干这个的）
		g[i] = n/i-1;
	}
	if (lim*lim<n)
	{
		for (LL i = lim+1; i<=(lim<<1); i++)
		{
			h[i] = (lim<<1)+1-i;
			//使用 map 存储 <(lim<<1)+1-i, i> 这个键值对,但是使用map这个stl会超时,所以使用sz来优化（sz就是用来干这个的）
			g[i] = (lim<<1)-i;
		}
	}
	else
	{
		for (LL i = lim+1; i<=(lim<<1)-1; i++)
		{
			h[i] = (lim<<1)-i;
			//使用 map 存储 <(lim<<1)-i, i> 这个键值对,但是使用map这个stl会超时,所以使用sz来优化（sz就是用来干这个的）
			g[i] = (lim<<1)-i-1;
		}
		flag = true; // n是完全平方数
	}
}

LL gettot() // 二分查找 <=lim的最大素数
{
	LL lo = 1, hi = top;
	while(lo < hi)
	{
		LL mid = (lo+hi)>>1;
		primes[mid]> lim?hi = mid:lo = mid+1;
	}
	return --lo;
}

LL solve()
{
	init();
	LL tot = gettot(); // <=lim的素数的个数
	for (LL i = 1; i<=tot; i++) // dp进行tot次, 因为最后想求的是 g(n,tot)
	{
		for (LL j = 1; h[j]>=primes[i]*primes[i]; j++) // 这是dp公式
		{
			int t = sz(h[j]/primes[i]); // sz 是根据值,反求在h中的索引（注意, 不能使用map/hash_map/unordered_map, 否则map存取将导致超时）
			g[j] -= g[t]-(i-1); // 这是dp公式
		}
	}
	return g[1];
}

int main() {

#ifdef LOCAL
	freopen("d:\\data.in", "r", stdin);
#endif
	sieve();
	scanf("%lld", &n);
	printf("%lld", solve());
	return 0;
}

ac情况

#571072 #6235. 区间素数个数 Accepted 100 5759 ms 12384 K C++ / 1.5 K yfs 2019-08-14 12:42:48

其实上面的代码可以有个优化, 就是不需要先筛出100w内的所有质数——只需要筛出[1,lim]中所有质数就行了(这其实仰仗的是只有一组测试用例). 而且这还避免了使用二分查找去查找<=lim的最大质数（因为它在筛[1,lim]的过程中就算出来了）.

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
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
//#include "stdafx.h"

#include <stdio.h>
#include <math.h>
#include <string.h>

//#define LOCAL

typedef long long LL;

LL n, lim, top;

const LL maxn = 2*1e6+5, maxnn = 1e6+5;

bool isprime[maxnn], flag;

LL g[maxn], h[maxn], primes[maxnn];

void sieve() // 筛出 [1,lim]中的所有质数,而不是筛出[1,100w]中的,因为后者我还要使用二分查找去查找<=lim的最大质数
{
	memset(isprime, true, sizeof(isprime));
	for (LL i = 2; i<=lim; i++)
	{
		if (isprime[i])
		{
			for (LL j = (i<<1); j<=lim; j+=i)
			{
				isprime[j] = false;
			}
		}
	}
	for (LL i = 2; i<=lim; i++)
	{
		if (isprime[i])
		{
			primes[++top] = i;
		}
	}
}

LL sz(LL x)
{
	if (flag) // 如果n是完全平方数
	{
		if (x>=lim)
		{
			return n/x;
		}
		return (lim<<1)-x;
	}
	if (x>lim)
	{
		return n/x;
	}
	return (lim<<1)+1-x;
}

void init()
{
	lim = sqrt(1.0*n);
	sieve();
	for (LL i = 1; i<=lim; i++)
	{
		h[i] = n/i;
		g[i] = n/i-1;
	}
	if (lim*lim<n)
	{
		for (LL i = lim+1; i<=(lim<<1); i++)
		{
			h[i] = (lim<<1)+1-i;
			g[i] = (lim<<1)-i;
		}
	}
	else
	{
		for (LL i = lim+1; i<=(lim<<1)-1; i++)
		{
			h[i] = (lim<<1)-i;
			g[i] = (lim<<1)-i-1;
		}
		flag = true; // n是完全平方数
	}
}

LL solve()
{
	init();
	for (LL i = 1; i<=top; i++)
	{
		for (LL j = 1; h[j]>=primes[i]*primes[i]; j++) // 这是dp公式
		{
			int t = sz(h[j]/primes[i]);
			g[j] -= g[t]-(i-1); // 这是dp公式
		}
	}
	return g[1];
}

int main() {

#ifdef LOCAL
	freopen("d:\\data.in", "r", stdin);
#endif
	scanf("%lld", &n);
	printf("%lld", solve());
	return 0;
}

ac情况

#571079 #6235. 区间素数个数 Accepted 100 5331 ms 10440 K C++ / 1.7 K yfs 2019-08-14 13:02:19

你看，妥妥的节约了 400ms之多~

完全类似的问题是 hdu 5901 Count primes (唯一的区别就是这里是读取多组数据，而loj是读取一组数据)

简单改一下，就a掉了

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
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
#include "stdafx.h"

#include <stdio.h>
#include <math.h>
#include <string.h>

#define LOCAL

typedef long long LL;

LL n, lim, top;

const LL maxn = 2*1e6+5, maxnn = 1e6+5;

bool isprime[maxn], flag;

LL g[maxn], h[maxn], primes[maxn];

void sieve() // 线性筛出 100万以内的质数
{
	memset(isprime, true, sizeof(isprime));
	for (LL i = 2; i<maxnn; i++)
	{
		if (isprime[i])
		{
			for (LL j = (i<<1); j<maxnn; j+=i)
			{
				isprime[j] = false;
			}
		}
	}
	for (LL i = 2; i<maxnn; i++)
	{
		if (isprime[i])
		{
			primes[++top] = i;
		}
	}
}

LL sz(LL x)
{
	if (flag) // 如果n是完全平方数
	{
		if (x>=lim)
		{
			return n/x;
		}
		return (lim<<1)-x;
	}
	if (x>lim)
	{
		return n/x;
	}
	return (lim<<1)+1-x;
}

void init() // N/1,N/2,...,N/lim,lim,...,2,1 （如果N是完全平方数的话, 则只保留N/lim和lim中的一个）
{
	flag = false;
	lim = sqrt(1.0*n);
	for (LL i = 1; i<=lim; i++)
	{
		h[i] = n/i;
		// 使用 map 存储 <n/i, i> 这个键值对，但是使用map这个stl会超时,所以使用sz来优化（sz就是用来干这个的）
		g[i] = n/i-1;
	}
	if (lim*lim<n)
	{
		for (LL i = lim+1; i<=(lim<<1); i++)
		{
			h[i] = (lim<<1)+1-i;
			//使用 map 存储 <(lim<<1)+1-i, i> 这个键值对,但是使用map这个stl会超时,所以使用sz来优化（sz就是用来干这个的）
			g[i] = (lim<<1)-i;
		}
	}
	else
	{
		for (LL i = lim+1; i<=(lim<<1)-1; i++)
		{
			h[i] = (lim<<1)-i;
			//使用 map 存储 <(lim<<1)-i, i> 这个键值对,但是使用map这个stl会超时,所以使用sz来优化（sz就是用来干这个的）
			g[i] = (lim<<1)-i-1;
		}
		flag = true; // n是完全平方数
	}
}

LL gettot() // 二分查找 <=lim的最大素数
{
	LL lo = 1, hi = top;
	while(lo < hi)
	{
		LL mid = (lo+hi)>>1;
		primes[mid]> lim?hi = mid:lo = mid+1;
	}
	return --lo;
}

LL solve()
{
	init();
	LL tot = gettot(); // <=lim的素数的个数
	for (LL i = 1; i<=tot; i++) // dp进行tot次, 因为最后想求的是 g(n,tot)
	{
		for (LL j = 1; h[j]>=primes[i]*primes[i]; j++) // 这是dp公式
		{
			int t = sz(h[j]/primes[i]); // sz 是根据值,反求在h中的索引（注意, 不能使用map/hash_map/unordered_map, 否则map存取将导致超时）
			g[j] -= g[t]-(i-1); // 这是dp公式
		}
	}
	return g[1];
}

int main() {

#ifdef LOCAL
	freopen("d:\\data.in", "r", stdin);
#endif
	sieve();
	while(~scanf("%lld", &n))
	{
		printf("%lld\n", solve());
	}
	return 0;
}

ac情况

30266815 2019-08-14 13:24:12 Accepted 5901 1762MS 13712K 2338 B G++ yfsyfsyfs

其实解题过程中，无耻的看了很多解题报告~ 但是dp公式都能看懂，就是不晓得代码为什么像他们那么写~ 最后还是放弃读懂他们的代码的念头（本人太弱了），自己理解了思想，自己写了 QAQ

参考

【1】https://yfsyfs.github.io/2019/08/11/lightoj-1197-Help-Hanzo-%E7%B4%A0%E6%95%B0%E7%AD%9B%E6%B3%95/