杜教筛学习笔记

杜教筛学习笔记

闲话

感觉以前根本没学懂杜教筛，于是重学了一遍，写个笔记记录一下。

前置知识

依赖于迪利克雷卷积、莫比乌斯反演、整除分块相关知识。

记号约定及基本性质

约定：

• $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;
}