杜教筛学习笔记
杜教筛学习笔记
闲话
感觉以前根本没学懂杜教筛，于是重学了一遍，写个笔记记录一下。
前置知识
依赖于迪利克雷卷积、莫比乌斯反演、整除分块相关知识。
记号约定及基本性质
约定：
- \(f*g\) 表示 \(f\) 与 \(g\) 的迪利克雷卷积，即 \((f*g)(n)=\sum\limits_{ij=n}f(i)g(j)\)。
- \(f\cdot g\) 表示 \(f\) 与 \(g\) 的点积，即 \((f\cdot g)(n)=f(n)g(n)\)。
- \(f^k\) 表示 \(f\) 点积的幂次，即 \(f^k=\underbrace{f\cdot f\cdot\cdots\cdot f}_{k\textrm{ times}}\)。
- \(S_f\) 表示 \(f\) 的前缀和，即 \(S_f(n)=\sum\limits_{i=1}^nf(i)\)。
- \(\epsilon(x)=[x=1]\)。
- \(\operatorname{I}(x)=1\)。
- \(\operatorname{id}(x)=x\)。
- \(\varphi\) 为欧拉函数，即 \(\varphi(n)=\sum\limits_{i=1}^n[i\perp n]\)。
- \(\mu\) 为莫比乌斯函数，是 \(\operatorname{I}\) 的反函数，设 \(n=\prod_{i=1}^kp_i^{\alpha_i}\)，即 \(\mu(n)=\begin{cases}1,&n=1\\0,&\exists i,\alpha_i\ge 2\\(-1)^k,&\forall i,\alpha_i=1\end{cases}\)。
- \(\sigma_k\) 为因数幂和函数，即 \(\sigma_k(n)=\sum\limits_{d\mid n}d^k\)。特别地，\(\sigma_0\) 为因数个数函数，\(\sigma_1\) 为因数和函数。
有性质：
- \(\varphi*\operatorname{I}=\operatorname{id}\)，即 \(\varphi=\mu*\operatorname{id}\)。
- \(\mu*\operatorname{I}=\epsilon\)。
- \((f\cdot g)*(f\cdot h)=f\cdot(g*h)\)。
杜教筛算法流程
杜教筛用于求一类数论函数的前缀和，并不要求积性。
假设要求 \(S_f(n)\)，如果构造出数论函数 \(g\) 满足 \(S_g,S_h\) 可以快速求出（记 \(h=f*g\)），即可快速计算出 \(S_f(n)\)。
进行以下推导：
\[\begin{aligned}
S_h(n)&=\sum_{i=1}^nh(i)\\
&=\sum_{i=1}^n\sum_{d\mid i}f\left(\frac{i}{d}\right)g(d)\\
&=\sum_{d=1}^n\sum_{i=1}^{\frac{n}{d}}f(i)g(d)\\
&=\sum_{d=1}^ng(d)S_f\left(\left\lfloor\frac{n}{d}\right\rfloor\right)\\
&=g(1)S_f(n)+\sum_{d=2}^ng(d)S_f\left(\left\lfloor\frac{n}{d}\right\rfloor\right)\\
S_f(n)&=\frac{1}{g(1)}\left(S_h(n)-\sum_{d=2}^ng(d)S_f\left(\left\lfloor\frac{n}{d}\right\rfloor\right)\right)
\end{aligned}
\]
利用可以快速求出的 \(S_h\) 和 \(S_g\)，可以整除分块求解。
复杂度证明
显然只会递归求 \(S_f(\lfloor\frac{n}{d}\rfloor)\) 的值，共 \(O(\sqrt{n})\) 个，可以记忆化之。
复杂度为：
\[\begin{aligned}
T(n)&=O\left(\sum_{i=1}^{\sqrt{n}}\sqrt{i}+\sum_{i=1}^{\sqrt{n}}\sqrt{\left\lfloor\frac{n}{i}\right\rfloor}\right)\\
&=O\left(\sum_{i=1}^{\sqrt{n}}\sqrt{\left\lfloor\frac{n}{i}\right\rfloor}\right)\\
&=O\left(\int_1^{\sqrt{n}}\sqrt{\frac{n}{x}}\textrm{d}x\right)\\
&=O\left(n^{\frac{3}{4}}\right)
\end{aligned}
\]
如果我们通过线性筛等方式预处理出 \(1\sim n^c\) 的答案（其中 \(c > \frac{1}{2}\)），那么只需要递归计算 \(\lfloor\frac{n}{1}\rfloor,\lfloor\frac{n}{2}\rfloor,\cdots,\lfloor\frac{n}{n^{1-c}}\rfloor\)，复杂度为：
\[\begin{aligned}
T(n)&=O\left(n^c+\sum_{i=1}^{n^{1-c}}\sqrt{\left\lfloor\frac{n}{i}\right\rfloor}\right)\\
&=O\left(n^c+\int_1^{n^{1-c}}\sqrt{\frac{n}{x}}\textrm{d}x\right)\\
&=O\left(n^c+n^{1-\frac{1}{2}c}\right)\\
\end{aligned}
\]
令 \(c=\frac{2}{3}\) 时达到最优复杂度 \(T(n)=O(n^{\frac{2}{3}})\)。实际应用时，考虑到多测等原因，一般不预处理到 \(O(n^{\frac{2}{3}})\)，而是预处理到可接受的范围（例如 \(10^7\)）。
记忆化时使用 map 的话多一个 \(\log\)，如果利用整除分块的性质在数组里记忆化，就不会带 \(\log\)。
杜教筛套杜教筛
有时求 \(S_g,S_h\) 时也需要进行杜教筛，但是不会影响复杂度，因为都只会递归求 \(\lfloor\frac{n}{d}\rfloor\) 处的值，总复杂度是相加关系。
例题
P4213 【模板】杜教筛（Sum）
单点求 \(S_\varphi\) 和 \(S_\mu\)。
对于 \(S_\varphi\)，使用 \(\varphi*\operatorname{I}=\operatorname{id}\)。显然有 \(S_{\operatorname{I}}(n)=n\)，\(S_{\operatorname{id}}(n)=\frac{n(n+1)}{2}\)。
对于 \(S_\mu\)，使用 \(\mu*\operatorname{I}=\epsilon\)。显然有 \(S_{\operatorname{I}}(n)=n\)，\(S_\epsilon(n)=1\)。
代码
//By: OIer rui_er
#include <bits/stdc++.h>
#define rep(x,y,z) for(ll x=(y);x<=(z);x++)
#define per(x,y,z) for(ll x=(y);x>=(z);x--)
#define debug(format...) fprintf(stderr, format)
#define fileIO(s) do{freopen(s".in","r",stdin);freopen(s".out","w",stdout);}while(false)
using namespace std;
typedef long long ll;
mt19937 rnd(std::chrono::duration_cast<std::chrono::nanoseconds>(std::chrono::system_clock::now().time_since_epoch()).count());
ll randint(ll L, ll R) {
uniform_int_distribution<ll> dist(L, R);
return dist(rnd);
}
template<typename T> void chkmin(T& x, T y) {if(x > y) x = y;}
template<typename T> void chkmax(T& x, T y) {if(x < y) x = y;}
const ll N = 1e7+5;
ll T, n, p[N / 10], pcnt, sphi[N], smu[N];
bool tab[N];
map<ll, ll> Sphi, Smu;
void sieve(ll lim) {
sphi[1] = 1;
smu[1] = 1;
rep(i, 2, lim) {
if(!tab[i]) {
p[++pcnt] = i;
sphi[i] = i - 1;
smu[i] = -1;
}
for(ll j = 1; j <= pcnt && i * p[j] <= lim; j++) {
tab[i * p[j]] = true;
if(i % p[j] != 0) {
sphi[i * p[j]] = sphi[i] * sphi[p[j]];
smu[i * p[j]] = -smu[i];
}
else {
sphi[i * p[j]] = sphi[i] * p[j];
smu[i * p[j]] = 0;
break;
}
}
}
rep(i, 2, lim) {
sphi[i] += sphi[i-1];
smu[i] += smu[i-1];
}
}
ll djs_phi(ll x) {
if(x < N) return sphi[x];
if(Sphi.count(x)) return Sphi[x];
ll ans = x * (x + 1) / 2;
for(ll L = 2, R; L <= x; L = R + 1) {
R = x / (x / L);
ans -= (R - L + 1) * djs_phi(x / L);
}
return Sphi[x] = ans;
}
ll djs_mu(ll x) {
if(x < N) return smu[x];
if(Smu.count(x)) return Smu[x];
ll ans = 1;
for(ll L = 2, R; L <= x; L = R + 1) {
R = x / (x / L);
ans -= (R - L + 1) * djs_mu(x / L);
}
return Smu[x] = ans;
}
int main() {
sieve(N-1);
for(scanf("%lld", &T); T; T--) {
scanf("%lld", &n);
printf("%lld %lld\n", djs_phi(n), djs_mu(n));
}
return 0;
}
P3768 简单的数学题
求 \(\left(\sum\limits_{i=1}^n\sum\limits_{j=1}^nij\gcd(i,j)\right)\bmod p\)。
易知：
\[\begin{aligned}
\sum_{i=1}^n\sum_{j=1}^nij\gcd(i,j)&=\sum_{i=1}^n\sum_{j=1}^nij\sum_{d\mid i,d\mid j}\varphi(d)\\
&=\sum_{d=1}^n\varphi(d)\sum_{i=1}^{\frac{n}{d}}id\sum_{j=1}^{\frac{n}{d}}jd\\
&=\sum_{d=1}^nd^2\varphi(d)\left(\frac{\left\lfloor\frac{n}{d}\right\rfloor\left(\left\lfloor\frac{n}{d}\right\rfloor+1\right)}{2}\right)^2
\end{aligned}
\]
设 \(f(x)=x^2\varphi(x)\)，只需快速求 \(S_f\) 即可整除分块求解。
注意到 \(f=\operatorname{id}^2\cdot\varphi\)，由前文性质构造出 \(g=\operatorname{id}^2\cdot\operatorname{I}=\operatorname{id}^2\)，则 \(h=f*g=\operatorname{id}^3\)。
显然有 \(S_g(n)=\frac{n(n+1)(2n+1)}{6}\)，\(S_h(n)=\left(\frac{n(n+1)}{2}\right)^2\)，杜教筛即可。
代码
//By: OIer rui_er
#include <bits/stdc++.h>
#define rep(x,y,z) for(ll x=(y);x<=(z);x++)
#define per(x,y,z) for(ll x=(y);x>=(z);x--)
#define debug(format...) fprintf(stderr, format)
#define fileIO(s) do{freopen(s".in","r",stdin);freopen(s".out","w",stdout);}while(false)
using namespace std;
typedef long long ll;
template<typename T> void chkmin(T& x, T y) {if(x > y) x = y;}
template<typename T> void chkmax(T& x, T y) {if(x < y) x = y;}
const ll N = 1e7+5;
ll mod, n, inv2, inv6, p[N], pcnt, sf[N];
bool tab[N];
map<ll, ll> Sf;
// f = id ` id ` phi
// g = id ` id ` I
// f * g = id ` id ` id
inline void sieve(ll lim) {
sf[1] = 1;
rep(i, 2, lim) {
if(!tab[i]) {
p[++pcnt] = i;
sf[i] = i - 1;
}
for(ll j = 1; j <= pcnt && i * p[j] <= lim; j++) {
tab[i * p[j]] = true;
if(i % p[j] != 0) {
sf[i * p[j]] = sf[i] * sf[p[j]] % mod;
}
else {
sf[i * p[j]] = sf[i] * p[j] % mod;
break;
}
}
}
rep(i, 1, lim) {
sf[i] = sf[i] * i % mod * i % mod;
sf[i] = (sf[i-1] + sf[i]) % mod;
}
}
inline ll S2(ll x) {x %= mod; return x * (x + 1) % mod * inv2 % mod;}
inline ll Sg(ll x) {x %= mod; return x * (x + 1) % mod * (2 * x + 1) % mod * inv6 % mod;}
inline ll Sh(ll x) {x %= mod; return S2(x) * S2(x) % mod;}
ll djs(ll x) {
if(x < N) return sf[x];
if(Sf.count(x)) return Sf[x];
ll ans = Sh(x);
for(ll L = 2, R; L <= x; L = R + 1) {
R = x / (x / L);
ans = (ans - (Sg(R) - Sg(L-1) + mod) % mod * djs(x / L) % mod + mod) % mod;
}
return Sf[x] = ans;
}
inline ll qpow(ll x, ll y) {
ll ans = 1;
for(; y; y >>= 1, x = x * x % mod) if(y & 1) ans = ans * x % mod;
return ans;
}
inline ll inv(ll x) {return qpow(x, mod-2);}
int main() {
scanf("%lld%lld", &mod, &n);
inv2 = inv(2); inv6 = inv(6);
sieve(N-1);
ll ans = 0;
for(ll L = 1, R; L <= n; L = R + 1) {
R = n / (n / L);
ans = (ans + Sh(n / L) * (djs(R) - djs(L-1) + mod) % mod) % mod;
}
printf("%lld\n", ans);
return 0;
}
浙公网安备 33010602011771号