杜教筛
杜教筛被用于处理一类数论函数的前缀和问题。对于数论函数 \(f\),杜教筛可以在低于线性时间的复杂度内计算 \(S(n)=\sum_{i=1}^{n}f(i)\)。
算法思想¶
我们想办法构造一个 \(S(n)\) 关于 \(S\left(\left\lfloor\frac{n}{i}\right\rfloor\right)\) 的递推式。
对于任意一个数论函数 \(g\),必满足:
其中 \(f*g\) 为数论函数 \(f\) 和 \(g\) 的 狄利克雷卷积。
略证
\(g(d)f\left(\frac{i}{d}\right)\) 就是对所有 \(i\leq n\) 的做贡献,因此变换枚举顺序,枚举 \(d\),\(\frac{i}{d}\)(分别对应新的 \(i,j\))
那么可以得到递推式:
假如我们可以构造恰当的数论函数 \(g\) 使得:
- 可以快速计算 \(\sum_{i=1}^n(f * g)(i)\);
- 可以快速计算 \(g\) 的单点值,以用数论分块求解 \(\sum_{i=2}^ng(i)S\left(\left\lfloor\frac{n}{i}\right\rfloor\right)\)。
则我们可以在较短时间内求得 \(g(1)S(n)\)。
注意
无论数论函数 \(f\) 是否为积性函数,只要可以构造出恰当的数论函数 \(g\), 便都可以考虑用杜教筛求 \(f\) 的前缀和。
如考虑 \(f(n)=\mathrm{i}\varphi(n)\), 显然 \(f\) 不是积性函数,但可取 \(g(n)=1\), 从而:
计算 \(\sum_k (f*g)(k)\) 和 \(g\) 的时间复杂度均为 \(O(1)\), 故可以考虑使用杜教筛。
时间复杂度¶
我们认为计算 \(\sum_{i=1}^n(f * g)(i)\) 和 \(g(n)\) 的时间复杂度均为 \(O(1)\), 设计算 \(S(n)\) 的复杂度为 \(T(n)\), 此时我们不妨将 \(S(n)\) 简化为如下形式:
由 整除分块/数论分块 可知 \(\left\lfloor \dfrac n i \right\rfloor\) 共有 \(O(\sqrt{n})\) 种取值,故有:
Note
\(O\left(\sum_{j=2}^{\sqrt{\frac{n}{i}}} T\left(\left\lfloor\frac{n}{ij}\right\rfloor\right)\right)\) 视作高阶无穷小,从而可以舍去。
故:
如果可以通过线性筛预处理出 \(S(1)\) 到 \(S(k)\) 的值,此时的 \(\sum_i S\left(\left\lfloor \dfrac{n}{i}\right\rfloor\right)\) 中我们只需要计算 \(k<\left\lfloor \dfrac{n}{i}\right\rfloor\leq n\) 的部分。设计算这一部分的复杂度为 \(T'(n)\),则有:
从而:
由均值不等式可知,当 \(k=\Theta\left(n^{\frac{2}{3}}\right)\) 时,\(T(n)\) 取得最小值 \(O\left(n^{\frac{2}{3}}\right)\).
例题¶
问题一¶
P4213【模板】杜教筛(Sum)
求 \(S_1(n)= \sum_{i=1}^{n} \mu(i)\) 和 \(S_2(n)= \sum_{i=1}^{n} \varphi(i)\) 的值,\(1\leq n<2^{31}\).
我们知道:
时间复杂度的推导见 时间复杂度 一节。
对于较大的值,需要用 map
/ unordered_map
存下其对应的值,方便以后使用时直接使用之前计算的结果。
当然也可以用杜教筛求出 \(\varphi (x)\) 的前缀和,但是更好的方法是应用莫比乌斯反演。
由于题目所求的是 \(\sum_{i=1}^n \sum_{j=1}^i [\gcd(i,j)=1]\), 所以我们排除掉 \(i=1,j=1\) 的情况,并将结果除以 \(2\) 即可。
观察到,只需求出莫比乌斯函数的前缀和,就可以快速计算出欧拉函数的前缀和了。时间复杂度 \(O\left(n^{\frac 2 3}\right)\).
求 \(S(n)=\sum_{i=1}^n\varphi(i)\).
同样的,\(\varphi * 1=\operatorname{id}\), 从而:
代码实现
#include <algorithm>
#include <cstdio>
#include <cstring>
#include <map>
using namespace std;
const int maxn = 2000010;
long long T, n, pri[maxn], cur, mu[maxn], sum_mu[maxn];
bool vis[maxn];
map<long long, long long> mp_mu;
long long S_mu(long long x) { // 求mu的前缀和
if (x < maxn) return sum_mu[x];
if (mp_mu[x]) return mp_mu[x]; // 如果map中已有该大小的mu值,则可直接返回
long long ret = (long long)1;
for (long long i = 2, j; i <= x; i = j + 1) {
j = x / (x / i);
ret -= S_mu(x / i) * (j - i + 1);
}
return mp_mu[x] = ret; // 路径压缩,方便下次计算
}
long long S_phi(long long x) { // 求phi的前缀和
long long ret = (long long)0;
long long j;
for (long long i = 1; i <= x; i = j + 1) {
j = x / (x / i);
ret += (S_mu(j) - S_mu(i - 1)) * (x / i) * (x / i);
}
return (ret - 1) / 2 + 1;
}
int main() {
scanf("%lld", &T);
mu[1] = 1;
for (int i = 2; i < maxn; i++) { // 线性筛预处理mu数组
if (!vis[i]) {
pri[++cur] = i;
mu[i] = -1;
}
for (int j = 1; j <= cur && i * pri[j] < maxn; j++) {
vis[i * pri[j]] = true;
if (i % pri[j])
mu[i * pri[j]] = -mu[i];
else {
mu[i * pri[j]] = 0;
break;
}
}
}
for (int i = 1; i < maxn; i++)
sum_mu[i] = sum_mu[i - 1] + mu[i]; // 求mu数组前缀和
while (T--) {
scanf("%lld", &n);
printf("%lld %lld\n", S_phi(n), S_mu(n));
}
return 0;
}
问题二¶
「LuoguP3768」简单的数学题
大意:求
其中 \(n\leq 10^{10},5\times 10^8\leq p\leq 1.1\times 10^9\),\(p\) 是质数。
利用 \(\varphi * 1=\operatorname{id}\) 做莫比乌斯反演化为:
其中 \(F(n)=\frac{1}{2}n(n+1)\)
对 \(\sum_{d=1}^nF\left(\left\lfloor\frac{n}{d}\right\rfloor\right)^2\) 做数论分块,\(d^2\varphi(d)\) 的前缀和用杜教筛处理:
需要构造积性函数 \(g\),使得 \(f\times g\) 和 \(g\) 能快速求和。
单纯的 \(\varphi\) 的前缀和可以用 \(\varphi * 1\) 的杜教筛处理,但是这里的 \(f\) 多了一个 \(\operatorname{id}^2\),那么我们就卷一个 \(\operatorname{id}^2\) 上去,让它变成常数:
化一下卷积:
再化一下 \(S(n)\):
分块求解即可。
代码实现
#include <cmath>
#include <cstdio>
#include <map>
using namespace std;
const int N = 5e6, NP = 5e6, SZ = N;
long long n, P, inv2, inv6, s[N];
int phi[N], p[NP], cnt, pn;
bool bp[N];
map<long long, long long> s_map;
long long ksm(long long a, long long m) { // 求逆元用
long long res = 1;
while (m) {
if (m & 1) res = res * a % P;
a = a * a % P, m >>= 1;
}
return res;
}
void prime_work(int k) { // 线性筛phi,s
bp[0] = bp[1] = 1, phi[1] = 1;
for (int i = 2; i <= k; i++) {
if (!bp[i]) p[++cnt] = i, phi[i] = i - 1;
for (int j = 1; j <= cnt && i * p[j] <= k; j++) {
bp[i * p[j]] = 1;
if (i % p[j] == 0) {
phi[i * p[j]] = phi[i] * p[j];
break;
} else
phi[i * p[j]] = phi[i] * phi[p[j]];
}
}
for (int i = 1; i <= k; i++)
s[i] = (1ll * i * i % P * phi[i] % P + s[i - 1]) % P;
}
long long s3(long long k) { // 立方和
return k %= P, (k * (k + 1) / 2) % P * ((k * (k + 1) / 2) % P) % P;
}
long long s2(long long k) { // 平方和
return k %= P, k * (k + 1) % P * (k * 2 + 1) % P * inv6 % P;
}
long long calc(long long k) { // 计算S(k)
if (k <= pn) return s[k];
if (s_map[k]) return s_map[k]; // 对于超过pn的用map离散存储
long long res = s3(k), pre = 1, cur;
for (long long i = 2, j; i <= k; i = j + 1)
j = k / (k / i), cur = s2(j),
res = (res - calc(k / i) * (cur - pre) % P) % P, pre = cur;
return s_map[k] = (res + P) % P;
}
long long solve() {
long long res = 0, pre = 0, cur;
for (long long i = 1, j; i <= n; i = j + 1) {
j = n / (n / i);
cur = calc(j);
res = (res + (s3(n / i) * (cur - pre)) % P) % P;
pre = cur;
}
return (res + P) % P;
}
int main() {
scanf("%lld%lld", &P, &n);
inv2 = ksm(2, P - 2), inv6 = ksm(6, P - 2);
pn = (long long)pow(n, 0.666667); // n^(2/3)
prime_work(pn);
printf("%lld", solve());
return 0;
} // 不要为了省什么内存把数组开小,会卡80