【数学】杜教筛

前置知识:莫比乌斯反演,狄利克雷卷积等亿点数论知识。

不会不要紧,这里顺便一提。

积性函数

对于一个数论函数 \(\mathbf f\),若对于任意的 \(a\perp b\),都有 \(\mathbf f\left(ab\right) = \mathbf f\left(a\right)\mathbf f\left(b\right)\),则称 \(f\) 为一个积性函数。即使不满足 \(a\perp b\) 时仍有 \(\mathbf f\left(ab\right) = \mathbf f\left(a\right)\mathbf f\left(b\right)\) 的数论函数被称为完全积性函数

一些积性函数的例子:

  • \(1\),1 函数。对于任意的 \(n\),有 \(1\left(n\right) = 1\)。(完全积性函数)
  • \(Id_k\),幂函数。\(Id_k\left(n\right) = n^k\),当 \(k = 1\) 时下标可省略不写。(完全记性函数)
  • \(\varphi\),欧拉函数。
  • \(\mu\) 莫比乌斯函数。
  • \(\varepsilon\)。单位元。\(\varepsilon\left(n\right) = [n = 1]\)。(完全积性函数)
  • \(\gcd\left(k, n\right)\)\(k\) 为常数)。

狄利克雷卷积

对于两个数论函数 \(\mathbf f\)\(\mathbf g\),定义其狄利克雷卷积为:

\[\mathbf {f *g}\left(n\right) = \sum_{d|n}\mathbf f\left(d\right)\mathbf g\left(\frac{n}{d}\right) \]

可以证明,\(*\) 运算符合交换律和结合律。

其单位元为 \(\varepsilon\)

两个积性函数的狄利克雷卷积仍旧是一个积性函数。

一些例子

  • \(\mathbf f * \varepsilon = \mathbf f\)
  • \(\mu * 1 = \varepsilon\)
  • \(\varphi * 1 = Id\)
  • \(Id * \mu = \varphi\)
  • \(1 * Id = \sigma\)
  • \(1 * 1 = d\)

后面两个分别是约数和函数和约数个数函数。

莫比乌斯反演

若两个数论函数 \(\mathbf f\)\(\mathbf g\) 满足如下关系

\[\mathbf f\left(n\right) = \sum_{d|n}\mathbf g\left(d\right) \]

那么就有

\[\mathbf g\left(n\right) = \sum_{d|n}\mathbf f\left(d\right)\mu\left(\frac{n}{d}\right) \]

证明:把上面的式子写成 \(\mathbf f = \mathbf g * 1\),然后两边同时卷上 \(\mu\) 即可。

杜教筛

杜教筛可以用于在非线性时间内求数论函数前缀和。

直接从具体例子入手,可能会好理解一些。

思考一个问题:\(\varphi\) 函数的前缀和怎么求?求出其 \(1\)\(n\) 的前缀和,\(n\)\(10^9\) 级别。

当数据范围很小的时候,我们的做法是通过线性筛筛出所有的 \(\varphi\) 函数值,然后求个前缀和。

然而在这种数据范围下,这种做法显然行不通了。

那我们来思考几个简单的问题:

  • \(1\) 的前缀和是什么?显然是 \(n\)
  • \(Id\) 的前缀和是什么?显然是 \(\frac{n\left(n + 1\right)}{2}\)
  • \(Id\left(1\right)\) 的值是多少?显然是 \(1\)

而我们考虑到 \(\varphi * 1 = Id\)

考虑能不能用这些很好求前缀和的数论函数来求 \(\mu\) 的前缀和。

\(\sum_{i = 1}^n\varphi\left(i\right) = s\left(n\right)\)

那么就有

\[s * 1 = \sum_{d|n}\sum_{i = 1}^d\varphi\left(i\right) \]

改变枚举顺序,考虑每个 \(1\left(d\right) =1\) 的贡献,那么原式等于

\[\sum_{d = 1}^n\sum_{i = 1}^{n / d}\varphi\left(i\right) \]

\[=\sum_{d = 1}^ns\left(\lfloor\frac{n}{d}\rfloor\right) \]

那么 \(s\left(n\right)\) 值是什么?答案是:

\[s\left(n\right) = \sum_{d = 1}^ns\left(\lfloor\frac{n}{d}\rfloor\right) - \sum_{d = 2}^ns\left(\lfloor\frac{n}{d}\rfloor\right) \]

\[= \sum_{i = 2}^n\left(\varphi * 1\right)\left(i\right) - \sum_{d = 1}^ns\left(\lfloor\frac{n}{d}\rfloor\right) \]

\[= \sum_{i = 2}^nId\left(i\right) - \sum_{d = 1}^ns\left(\lfloor\frac{n}{d}\rfloor\right) \]

\[= \dfrac{n\left(n + 1\right)}{2} - \sum_{d = 1}^ns\left(\lfloor\frac{n}{d}\rfloor\right) \]

哇,我们发现了什么?

后面的那个 \(s\) 我们可以在一边数论分块,一边递归的求。

这样我们就可以筛出一些比较小的时候的 \(s\) 的值,然后用个 map 记忆化一下计算。

LL sum_phi(LL n) {
	if(n <= N) return sphi0[(int)n];
	if(sphi.count(n)) return sphi[n];
	LL res = 1ll * n * (n + 1) / 2;
	for(LL l = 2, r = 1; l <= n; l = r + 1) {
		r = n / (n / l);
		res -= 1ll * (r - l + 1) * sum_phi(n / l);
	}
	return sphi[n] = res;
}

推广到一般情况,则可以得出下面的式子

对于一个数论函数 \(\mathbf f\),设 \(s\left(n\right) = \sum_{i = 1}^n\mathbf f\left(i\right)\),另外有数论函数 \(\mathbf g\)

\[s\left(n\right) = \dfrac{\sum_{i = 1}^n\left(\mathbf {f * g}\right)\left(i\right) - \sum_{i = 2}^n\mathbf g\left(i\right)s\left(\lfloor\frac{n}{i}\rfloor\right)}{\mathbf g\left(1\right)} \]

如果您能够快速求出 \(\left(\mathbf {f * g}\right)\) 的前缀和,\(\mathbf g\) 的前缀和,以及 \(\mathbf g\left(1\right)\) 的值的话,也就意味着您能够用类似于上面求 \(\sum_{i = 1}^n \varphi\left(i\right)\) 的方式,快速求出 \(\mathbf f\) 的前缀和!

算法的复杂度是 \(\mathcal O(n^{\frac{2}{3}})\) 的。

这就是杜教筛啦。

完整代码 洛谷 P4213 【模板】杜教筛

#include <bits/stdc++.h>
#define LL long long

template <typename Temp> inline void read(Temp & res) {
	Temp fh = 1; res = 0; char ch = getchar();
	for(; !isdigit(ch); ch = getchar()) if(ch == '-') fh = -1;
	for(; isdigit(ch); ch = getchar()) res = (res << 3) + (res << 1) + (ch ^ '0');
	res = res * fh;
}

namespace Math {
	#define N (10000000)
	const int Maxn = 1e7 + 5;
	bool isprime[Maxn]; int prime[Maxn], cnt_prime; LL phi0[Maxn], mu0[Maxn];
	void sieve() {
		phi0[1] = mu0[1] = 1;
		for(int i = 2; i <= N; ++i) {
			if(!isprime[i]) {
				prime[++cnt_prime] = i;
				phi0[i] = i - 1; mu0[i] = -1;
			}
			for(int j = 1; j <= cnt_prime && prime[j] <= N / i; ++j) {
				int cur = i * prime[j];
				isprime[cur] = 1;
				if(i % prime[j] == 0) {
					phi0[cur] = phi0[i] * prime[j];
					mu0[cur] = 0; break;
				} else {
					phi0[cur] = phi0[i] * (prime[j] - 1);
					mu0[cur] = -mu0[i];
				}
			}
		}
		for(int i = 1; i <= N; ++i) phi0[i] += phi0[i - 1], mu0[i] += mu0[i - 1];
	}
	
	std::map<LL, LL> mu, phi;
	//mu * 1 = epsilon
	LL sum_mu(LL n) {
		if(n <= N) return mu0[(int)n];
		if(mu.count(n)) return mu[n];
		LL res = 1ll;
		for(LL l = 2, r = 1; l <= n; l = r + 1) {
			r = n / (n / l);
			res -= 1ll * (r - l + 1) * sum_mu(n / l);
		}
		return mu[n] = res;
	}
	//phi * 1 = id 
	LL sum_phi(LL n) {
		if(n <= N) return phi0[(int)n];
		if(phi.count(n)) return phi[n];
		LL res = 1ll * n * (n + 1) / 2;
		for(LL l = 2, r = 1; l <= n; l = r + 1) {
			r = n / (n / l);
			res -= 1ll * (r - l + 1) * sum_phi(n / l);
		}
		return phi[n] = res;
	}
	#undef N
}

int T, n;

int main() {
	using namespace Math;
	sieve(); 
	read(T);
	while(T--) {
		read(n);
		printf("%lld %lld\n", sum_phi(n), sum_mu(n));
	}
	return 0;
}
posted @ 2021-03-21 18:24  zimujun  阅读(35)  评论(0编辑  收藏  举报