浅谈几种筛法

Algor_pro_king_John

于 2019-05-03 09:00:33 发布

阅读量658

点赞数 1

本文链接：https://blog.csdn.net/Algor_pro_king_John/article/details/89784763

版权

杜教筛

问题一般是求 $\sum_{i=1}^{n}f(i)$ 这样的式子。

然后我们有一种很妙的想法，那就是构造两个积性函数 $h, g$ ，使得 $h = f * g$

然后尝试推一下 $h$ 的前缀和，发现：
$\sum_{i=1}^{n}h(i)=\sum_{i=1}^{n}\sum_{d|i}g(d)\cdot f(\frac{i}{d})\\=\sum_{d=1}^{n}g(d)\cdot\sum_{i=1}^{\lfloor\frac{n}{d}\rfloor}f({i})$

如果记 $S(n)=\sum_{i=1}^nf(i)$ ，那么就可以写成： $\sum_{i=1}^{n}h(i)=\sum_{d=1}^{n}g(d)\cdot S(\lfloor\frac{n}{d}\rfloor)$

把 $d = 1$ 提出来，得到我们想要的式子： $g(1)S(n)=\sum_{i=1}^{n}h(i)-\sum_{d=2}^{n}g(d)\cdot S(\lfloor\frac{n}{d}\rfloor)$

那么如果 $h$ 的前缀和我们可以在一个比较优秀的时间复杂度内求出，则后面那一部分的时间复杂度在进行整除分块后，可以达到一个优秀的时间复杂度 $O(n^{\frac{2}{3}})$ 。

一些题目

51nod1244

题意：

求 $\sum _{i=a}^{b} \mu(i)（a\le b\le 10^{10}）$
解法一：

直接套用上面的卷积。我们知道 $\mu * 1=\epsilon$

且 $1,\epsilon$ 都是积性函数，所以可以直接得到： $\displaystyle S(n)=1 - \sum _{i=2}^{n} S(\lfloor \frac{n}{i} \rfloor)$
解法二：

还是来推一下式子，其实思路是一样的，都用到了 $\sum_{d|n}\mu(d)=[n=1]$ 这条式子。
$\sum_{i=1}^{n} \sum_{d|i}\mu(d)=1\\ \therefore \displaystyle \sum_{i=1}^{n} \sum_{j=1}^{\lfloor \frac{n}{i}\rfloor}\mu(j)=1$
提出 $j = 1$ ，得到 $\sum _{i=1}^{n} \mu(i)=1-\sum_{i=2}^{n} \sum_{j=1}^{\lfloor \frac{n}{i}\rfloor}\mu(j)$ 然后就转化为上面的式子。
Trick

这里有几个需要注意的地方。一个是我们可以预处理一些 $\le \sqrt{n}$ 的前缀和来加快速度（事实上，根据唐教的分析，如果不预处理，时间复杂度是 $O(n^{\frac{3}{4}})$ 的），另外一个就是hash时候的技巧，与 $M i n 25$ 筛不同，杜教筛的筛法需要用map记录一些已经算过的前缀和来保证时间复杂度，而这个时候我们其实可以不用map，而是用两个表去处理，按照是否 $\le \sqrt{n}$ 来处理。时间复杂度优秀许多，具体参见代码。

$code_1$ （用map储存）：

注意map <LL,LL> :: iterator iter，iter = h.find(x)，iter->second的用法

#include <iostream>
#include <cstring>
#include <cstdio>
#include <map>

#define F(i, a, b) for (int i = a; i <= b; i ++)

const int B = 5000000;

typedef long long LL;

using namespace std;

LL a, b, S[B + 1], mu[B]; int z[B];
bool bz[B + 1];

map <LL, LL> H;
map <LL, LL> :: iterator iter;

void Init() {
	mu[1] = 1;
	F(i, 2, B) {
		if (!bz[i])
			z[++ z[0]] = i, mu[i] = - 1;
		F(j, 1, z[0]) {
			LL x = 1ll * z[j] * i;
			if (x > B) break;
			bz[x] = 1;
			if (i % z[j] == 0) {
				mu[x] = 0;
				break;
			}
			mu[x] = mu[z[j]] * mu[i];
		}
	}
	F(i, 2, B)
		mu[i] += mu[i - 1];
}

LL Solve(LL x) {
	if (x <= B) return mu[x];
	iter = H.find(x);
	LL ans = iter -> second; //我竟然第一次知道要这样用。。。
	if (iter == H.end()) {
		ans = 0;
		for (LL i = 2, j; i <= x; i = j + 1) {
			j = x / (x / i);
			ans += (j - i + 1) * Solve(x / i);
		}
		ans = 1 - ans;
		H[x] = ans;
	}
	return ans;
}

int main() {
	scanf("%lld%lld", &a, &b);

	Init();
	
	printf("%lld\n", Solve(b) - Solve(a - 1));
}

$code_2$ （推荐储存方式）：

#include <iostream>
#include <cstring>
#include <cstdio>
#include <map>

#define F(i, a, b) for (int i = a; i <= b; i ++)

const int B = 5000000;

typedef long long LL;

using namespace std;

LL a, b, n, S[B + 1], vis[B / 100], mu[B];
int z[B]; bool bz[B + 1];

map <LL, LL> H;
map <LL, LL> :: iterator iter;

void Init() {
	mu[1] = 1;
	F(i, 2, B) {
		if (!bz[i])
			z[++ z[0]] = i, mu[i] = - 1;
		F(j, 1, z[0]) {
			LL x = 1ll * z[j] * i;
			if (x > B) break;
			bz[x] = 1;
			if (i % z[j] == 0) {
				mu[x] = 0;
				break;
			}
			mu[x] = mu[z[j]] * mu[i];
		}
	}
	F(i, 2, B)
		mu[i] += mu[i - 1];
}

LL Solve(LL x) {
	if (x <= B) return mu[x];
	LL y = n / x;
	if (!vis[y]) {
		LL ans = 0;
		for (LL i = 2, j; i <= x; i = j + 1) {
			j = x / (x / i);
			ans += (j - i + 1) * Solve(x / i);
		}
		vis[y] = ans = 1 - ans;
	}
	return vis[y];
}

LL Doit(LL x) {
	memset(vis, 0, sizeof vis), n = x;
	return Solve(x);
}

int main() {
	scanf("%lld%lld", &a, &b);

	Init();
	
	printf("%lld\n", Doit(b) - Doit(a - 1));
}

Min25筛

也是一种很优秀的求积性函数前缀和的筛法。~~据zzq大神说它的时间复杂度是 $O(\frac{n}{log_nlog_n})$ ？？~~

为了下面叙述方便，记 $P$ 表示素数集合，记 $P_i$ 表示第 $i$ 个质数，其中 $|P|=\sqrt{n}$ 。

此类问题的 $F (i)$ 一般满足：

$\begin{cases} f(p)是一个关于p的多项式\ \ \ \ p\in P\\ f(p^c)可以快速计算出\ \ \ \ \ \ \ \ \ \ \ \ otherwise\end{cases}$

则min25筛的想法是把所有的 $F (i)$ 的和分成质数和合数去考虑，即总的贡献 = 质数的贡献 + 合数的贡献，下面讲解两种方法：

递归版本

那么先考虑质数的贡献，即 $i$ 为质数，且保证括号内条件成立时的情况下的贡献。

若记 $m i n (p)$ 表示 $i$ 的最小质因子以及 $g(n,j)=\sum_{i=1}^n[i\in P, min(p)\gt P_j]F(i)$

则由定义，不难推出 $g(n,j)=\begin{cases} g(n,j-1)&{P_j}^2\gt n\\ g(n,j-1)-F(P_j)[g(\frac{n}{P_j},j-1)-\sum_{i=1}^{j-1}F(P_i)]&{P_j}^2\le n\end{cases}$

其中 $g (n, 0)$ 初值比较难理解。事实上，我们观察这整一个递推过程，实际上很类似埃氏筛，每一次都把一个 $P_j$ 的贡献给删去，最终筛到 $\sqrt{n}$ 时，保证了还没有被质数筛去的数在 $n$ 范围内一定不是合数。

那么继续考虑埃氏筛，实际上一开始我们是把所有数都看成质数的。那么 $g (n, 0)$ 的初值也就不难想了，可能是某个与 $n$ 有关的多项式，例如在求质数个数中 $f (i) = 1$ ，那么 $g (n, 0) = n - 1$ .

现在回到原问题 $\sum_{i=1}^nF(i)$ 上，我们要求质数的贡献，实质上就是 $g(n,|P|)-g(P_{j-1},j-1)$

如果记 $S(n,j)=\sum_{i=2}^nF(i)[i的最小质因子大于等于P_j]$

那么可以推出 $S(n,j)=g(n,|P|)-\sum_{i=1}^{j-1}f(P_i)+\sum_{k\ge j}^{|P|}\sum_{e}(f({P_k}^e)S(\frac{n}{{P_k}^e},k+1)+f({P_k}^{e+1}))$

事实上，后面那一部分式子就是合数的贡献，枚举质数 $k$ ，次数 $e$ ，然后与求 $g$ 的思想是完全类似的，注意要加上一个 $f({P_k}^{e+1})$ 就好了.

注意，这里我们所说的递归版本，实际上指的是 $s$ 的求法是用递归，实际实现中， $g$ 的求法可以不用递归，具体请看代码，并且求 $s$ 时可以不用记忆化，时间复杂度依旧优秀。

递推版本

实际上递推版本比递归更慢。因为它求的是 $\sqrt{n}$ 个 $\frac{n}{i}$ 的前缀和。

具体细节与递归版本类似。我们直接设 $S'(n,i)=\sum_{i=2}^nF(i)[i的最小质因子不小于p_i或i是质数]$ .

可以得到： $S'(n,i)=\left\{ \begin{array}{} S'(n,i+1)+\sum_{e\geq 1且\ p_i^{e+1}\ \leq n}\ \ F(p_i^e)[S(\lfloor\frac{n}{p_i^e}\rfloor,i+1)-g(p_i,i)]+F(p_i^{e+1})& &{p_i^2\leq n} \\ S'(n,i+1)& &{p_i^2\geq n} \end{array}\right.$

初值 $S^{'} (n, ∣ P ∣ + 1) = g (n, ∣ P ∣)$

一些题目

LOJ#6235. 区间素数个数

事实上就是求 $1\sim n$ 素数，这个可以直接通过求的 $g$ 来得到，连 $s$ 都不用求。

注意实现的时候只需要保存 $\frac{n}{i}$ 的 $\sqrt{n}$ 个值即可。

#include <bits/stdc++.h>

#define F(i, a, b) for (int i = a; i <= b; i ++)

const int N = 700000;

using namespace std;

long long n, m, id1[N], id2[N], w[N], g[N], z[N], sqr;
bool bz[N];

int main() {
	scanf("%lld", &n), sqr = int(sqrt(n));

	F(i, 2, sqr) {
		if (!bz[i])
			z[++ z[0]] = i, bz[i] = 1;
		F(j, 1, z[0]) {
			if (z[j] * i > sqr) break;
			bz[z[j] * i] = 1;
			if (i % z[j] == 0) break;
		}
	}

	for (long long i = 1, j; i <= n; i = j + 1) {
		j = n / (n / i); w[++ m] = n / i;
		if (w[m] <= sqr) id1[w[m]] = m; else id2[n / w[m]] = m;
		g[m] = w[m] - 1;
	}
	F(j, 1, z[0])
		for (int i = 1; i <= m && z[j] * z[j] <= w[i]; i ++) {
			int k = (w[i] / z[j] <= sqr) ? id1[w[i] / z[j]] : id2[n / (w[i] / z[j])];
			g[i] -= g[k] - (j - 1);
		}

	printf("%lld\n", g[1]);
}

#6053. 简单的函数

注意把 $f (p) = p - 1$ 拆成两个函数 $g (p) = p, h (p) = 1$ ，然后相减.

这样可以保证所求是积性函数.

#include <bits/stdc++.h>
#define F(i, a, b) for (int i = a; i <= b; i ++)

typedef long long LL;

const int N = 3e5;
const LL ny = 5e8 + 4, Mo = 1e9 + 7;

using namespace std;

bool bz[N]; int id1[N], id2[N], sqr;
LL g[N], h[N], z[N], w[N], s[N], n, m;

LL S(LL x, int y) {
	LL k = x <= sqr ? id1[x] : id2[n / x];
	LL sum = (g[k] - h[k] - (s[y - 1] - (y - 1)) + (y == 1) * 2) % Mo;
	for (LL i = y; i <= z[0] && z[i] * z[i] <= x; i ++)
		for (LL e = 1, p = z[i]; p * z[i] <= x; e ++, p *= z[i])
			sum = (sum + (z[i] ^ e) % Mo * S(x / p, i + 1) + (z[i] ^ (e + 1))) % Mo;
	return sum;
}

int main() {
	scanf("%lld", &n), sqr = int(sqrt(n));
	F(i, 2, sqr) {
		if (!bz[i])
			z[++ z[0]] = i, s[z[0]] = (s[z[0] - 1] + i) % Mo;
		F(j, 1, z[0]) {
			if (z[j] * i > sqr) break;
			bz[z[j] * i] = 1;
			if (i % z[j] == 0) break;
		}
	}
	for (LL i = 1, j; i <= n; i = j + 1) {
		j = n / (n / i); w[++ m] = n / i;
		if (w[m] <= sqr) id1[w[m]] = m; else id2[n / w[m]] = m;
		g[m] = (((w[m] + 2) % Mo * ((w[m] - 1) % Mo)) % Mo * ny) % Mo;
		h[m] = (w[m] - 1) % Mo;
	}
	F(j, 1, z[0])
		for (LL i = 1; i <= m && z[j] * z[j] <= w[i]; i ++) {
			LL k = (w[i] / z[j] <= sqr) ? id1[w[i] / z[j]] : id2[n / (w[i] / z[j])];
			h[i] = (h[i] - (h[k] - (j - 1))) % Mo;
			g[i] = (g[i] - z[j] * (g[k] - s[j - 1])) % Mo;
		}

	printf("%lld\n", n == 1 ? 1 : ((S(n, 1) + 1) % Mo + Mo) % Mo);
}

两种版本的时间复杂度都是 $O(\frac{n^{\frac{3}{4}}}{log(\sqrt n)})$ 。经实测，递归版会稍微快一些。

jzoj6027. 【GDOI2019模拟2019.2.23】签到

经过一系列的容斥后，可以计算。

比如说，要求 $1\sim n$ 中不是 $a_i(1\le i\le m)$ 的倍数，但要是 $b$ 的倍数。

我们就把 $n / b$ 之后进行计算。

再譬如，如果在此基础上，还要求 $f (x) = 1$ ，我们就要特殊考虑了，必须要考虑的是 $b$ 的质因数情况，因为我们还是在 $n / b$ 的基础上进行计算，经过一番特殊处理后还是可以做出来，只要透彻理解质数和合数的计算即可。

#include <bits/stdc++.h>

#define F(i, a, b) for (LL i = a; i <= b; i ++)
#define mem(a, b) memset(a, b, sizeof a)
#define get getchar()

typedef long long LL;

const LL M = 3e5 + 10;

using namespace std;

LL n, m, b, sqr, MAX, top, L, r[61], vis[M], s[M], S[M];
LL a[M], g[M], w[M], bz[M], z[M], id1[M], id2[M], R[M];
LL S0, S1, S2, S3; bool flag;

void Re(LL &x) {
	char c = get; x = 0; LL t = 1;
	for (; !isdigit(c); c = get) t = (c == '-' ? - 1 : t);
	for (; isdigit(c); x = (x << 3) + (x << 1) + c - '0', c = get); x *= t;
}

void Init() {
	Re(n), Re(m), Re(b);
	F(i, 1, m) Re(a[i]);
	F(i, 1, 60) Re(r[i]);
	r[0] = 1;
}

void GetPrime() {
	mem(bz, 0), z[0] = 0, mem(vis, 0), top = 0;
	F(i, 2, sqr) {
		if (!bz[i]) {
			z[++ z[0]] = i, s[z[0]] = s[z[0] - 1] + 1;
			vis[z[0]] = 1;
			while (top < m && a[top] < i) top ++;
			if (a[top] == i)
				s[z[0]] --, vis[z[0]] = 0;
		}
		F(j, 1, z[0]) {
			if (z[j] * i > sqr) break;
			bz[z[j] * i] = 1;
			if (i % z[j] == 0) break;
		}
	}
}

LL dg(LL x, LL y, LL n, LL t) {
	if (z[y] > x) return 0;
	LL k = (x <= sqr) ? id1[x] : id2[n / x];
	LL sum = (g[k] - s[y - 1]) * r[t];
	F(k, y, z[0]) {
		if (z[k] * z[k] > x) break;
		if (vis[k])
			for (LL s = z[k], cnt = 1; s * z[k] <= x; s *= z[k], cnt ++)
				if (!t) sum += dg(x / s, k + 1, n, t) + 1; else {
					if (r[cnt])
						sum += dg(x / s, k + 1, n, t);
					sum += r[cnt + 1];
				}
	}
	return sum;
}
LL DG(LL x, LL y, LL n,  LL res) {
	if (z[y] > x) return 0;
	LL k = (x <= sqr) ? id1[x] : id2[n / x], sum;
	if (res > 1) sum = 0; else
	if (res == 1) sum = (x >= MAX); else
		sum = (g[k] - s[y - 1]);
	F(k, y, z[0]) {
		if (z[k] * z[k] > x) break;
		if (vis[k]) {
			for (LL s = z[k], cnt = 1; s * z[k] <= x; s *= z[k], cnt ++) {
				if (r[cnt + R[k]])
					sum += DG(x / s, k + 1, n, res - !r[R[k]]);
				if ((res - !r[R[k]]) == 0)
					sum += r[cnt + R[k] + 1];
			}
			if (!r[R[k]]) break;
		}
	}
	return sum;
}

LL Solve(LL n, LL q) {
	sqr = q < 2 ? LL(sqrt(n)) : LL(sqrt(max(n, b))), L = 0, GetPrime();
	for (LL i = 1, j; i <= n; i = j + 1) {
		j = n / (n / i); w[++ L] = n / i;
		if (w[L] <= sqr) id1[w[L]] = L; else id2[n / w[L]] = L;
		g[L] = w[L] - 1;
	}
	F(j, 1, z[0])
		for (LL i = 1; i <= L && z[j] * z[j] <= w[i]; i ++) {
			LL k = (w[i] / z[j]) <= sqr ? id1[w[i] / z[j]] : id2[n / (w[i] / z[j])];
			g[i] -= (g[k] - (j - 1));
		}
	LL j = m;
	F(i, 1, L) {
		while (a[j] > w[i]) j --;
		g[i] -= j;
	}
	if (q <= 1)
		return dg(n, 1, n, q);
	else {
		LL x = b, res = 0;
		F(i, 1, z[0]) {
			while (x % z[i] == 0)
				R[i] ++, x /= z[i];
			s[i] = s[i - 1] + r[R[i] + 1] * vis[i];
			S[i] = S[i - 1] + vis[i];
			if (vis[i] * (!r[R[i]]))
				res ++, MAX = z[i];
		}
		if (x > 1) {
			z[++ z[0]] = x, vis[z[0]] = 1;
			while (top < m && a[top] < x) top ++;
			if (a[top] == x) vis[z[0]] = 0;
			s[z[0]] = s[z[0] - 1] + r[2] * vis[z[0]];
			S[z[0]] = S[z[0] - 1] + vis[z[0]];
			if (vis[z[0]] * (!r[1]))
				res ++, MAX = x;
		}
		LL j = z[0];
		F(i, 1, L) {
			while (j > 0 && z[j] > w[i]) j --;
			g[i] = (g[i] - S[j]) * (r[1] != 0) + s[j];
		}
		return DG(n, 1, n, res) + (res == 0);
	}
}

int main() {
	freopen("qiandao.in", "r", stdin);
	freopen("qiandao.out", "w", stdout);

	Init();
	S0 = Solve(n, 0);
	flag = 1;
	F(i, 1, m)
		if (b % a[i] == 0) { flag = 0; break; }
	S1 = flag == 0 ? 0 : Solve(n / b, 0) + (b > 1);

//	printf("%lld\n", S0);
//	printf("%lld\n", S1);

	S2 = Solve(n, 1);
	S3 = flag == 0 ? 0 : Solve(n / b, 2);
//	printf("%lld\n", S2);
//	printf("%lld\n", S3);

	printf("%lld\n", n - (n - S0 + S1 + S2 - S3));
}

Algor_pro_king_John

关注

1
点赞
踩
0

收藏

觉得还不错? 一键收藏
1
评论
浅谈几种筛法

杜教筛问题一般是求∑i=1nf(i)\sum_{i=1}^{n}f(i)i=1∑nf(i)这样的式子。然后我们有一种很妙的想法，那就是构造两个积性函数h,gh,gh,g，使得h=f∗gh=f*gh=f∗g然后尝试推一下hhh的前缀和，发现：∑i=1nh(i)=∑i=1n∑d∣ig(d)⋅f(id)=∑d=1ng(d)⋅∑i=1⌊nd⌋f(i)\sum_{i=1}^{n}h(i)=\su...
复制链接

扫一扫