杜教筛学习笔记

· · 个人记录

杜教筛可以在亚线性时间求出一些数论函数的前缀和。

\

基本思路

比如我们要求出 \sum_{i=1}^n f(i)，设 S(n) = \sum_{i=1}^n f(i)。

我们考虑一个函数 g 与 f 的狄利克雷卷积前缀和：

\begin{aligned} \sum_{i=1}^n (f * g) (i) &= \sum_{i=1}^n \sum_{d | i} g(d) f\left( \frac{i}{d} \right)\\ &= \sum_{d=1}^n \sum_{i=1}^{\lfloor \frac{n}{d} \rfloor} g(d) f(i)\\ &= \sum_{d=1}^n g(d) \sum_{i=1}^{\lfloor \frac{n}{d} \rfloor} f(i)\\ &= \sum_{i=1}^n g(i) S\left( \lfloor \frac{n}{i} \rfloor \right) \end{aligned}

所以有：

\begin{aligned} g(1) S\left( n \right) &= \sum_{i=1}^n g(i) S\left( \lfloor \frac{n}{i} \rfloor \right) - \sum_{i=2}^n g(i) S\left( \lfloor \frac{n}{i} \rfloor \right)\\ &= \sum_{i=1}^n (f * g) (i) - \sum_{i=2}^n g(i) S\left( \lfloor \frac{n}{i} \rfloor \right)\\ \end{aligned}

如果我们能快速求出 f 和 g 的狄利克雷卷积，还能快速求出 g(n)，那么我们就可以递归计算这个式子。数论分块即可。

\

模板：P4213 【模板】杜教筛

给定正整数 n，求 \sum_{i=1}^n \mu(i) 和 \sum_{i=1}^n \varphi(i)。

对于 \mu，我们知道 \mu * 1 = \epsilon，所以有：

\begin{aligned} 1 \times S\left( n \right) = \sum_{i=1}^n \epsilon (i) - \sum_{i=2}^n 1 \times S\left( \lfloor \frac{n}{i} \rfloor \right)\\ S\left( n \right) = 1 - \sum_{i=2}^n S\left( \lfloor \frac{n}{i} \rfloor \right) \end{aligned}

对于 \varphi，我们知道 \varphi * 1 = id，所以有：

\begin{aligned} 1 \times S\left( n \right) = \sum_{i=1}^n i - \sum_{i=2}^n 1 \times S\left( \lfloor \frac{n}{i} \rfloor \right)\\ S\left( n \right) = \frac{n(n+1)}{2} - \sum_{i=2}^n S\left( \lfloor \frac{n}{i} \rfloor \right) \end{aligned}

但是也可以用莫比乌斯反演。还记得仪仗队这道题吗？我们有一个式子：

\begin{aligned} \sum_{i=1}^n \sum_{j=1}^n [\gcd(i, j) = 1] &= 2 \sum_{i=1}^n \varphi(i) - 1\\ &= 2 S(n) - 1 \end{aligned}

所以：

S(n) = \frac{\sum_{i=1}^n \sum_{j=1}^n [\gcd(i, j) = 1] + 1}{2}

所以只要求出 \sum_{i=1}^n \sum_{j=1}^n [\gcd(i, j) = 1] 可以了。

用莫比乌斯反演：

\begin{aligned} \sum_{i=1}^n \sum_{j=1}^n [\gcd(i, j) = 1] &= \sum_{i=1}^n \sum_{j=1}^n \sum_{d | i, d | j} \mu(d)\\ &= \sum_{d=1}^n \mu(d) {\lfloor \frac{n}{d} \rfloor}^2 \end{aligned}

所以也可以用 \mu 的前缀和来计算，同样是数论分块。

CODE：

#include <iostream>
#include <cstdio>
#include <map>
using namespace std;
typedef long long ll;

const int maxn = 2e6 + 5;
const int maxk = 2e6;

int tot;
int mu[maxn], prime[maxn>>1];
ll smu[maxn];
map<int, ll> mp;
bool not_prime[maxn];

void prework() {
    mu[1] = 1;
    for(int i = 2; i <= maxk; i++) {
        if(!not_prime[i]) prime[++tot] = i, mu[i] = -1;
        for(int j = 1; j <= tot && i*prime[j] <= maxk; j++) {
            not_prime[i*prime[j]] = true;
            if(i % prime[j] == 0) {
                mu[i*prime[j]] = 0;
                break;
            }
            mu[i*prime[j]] = -mu[i];
        }
    }
    for(int i = 1; i <= maxk; i++)
        smu[i] = smu[i-1] + mu[i];
}

ll summu(int n) {
    if(n <= maxk) return smu[n];
    if(mp[n]) return mp[n];
    ll res = 1;
    for(ll l = 2, r; l <= n; l = r+1) {
        r = n/(n/l);
        res -= summu(n/l) * (r-l+1);
    }
    return mp[n] = res;
}

ll sumphi(int n) {
    ll res = 0;
    for(ll l = 1, r; l <= n; l = r+1) {
        r = n/(n/l);
        res += (summu(r) - summu(l-1)) * (n/l) * (n/l);
    }
    return (res + 1) >> 1ll;
}

int main() {
    prework();

    int T, n;
    scanf("%d", &T);
    while(T--) {
        scanf("%d", &n);
        printf("%lld %lld\n", sumphi(n), summu(n));
    }
    return 0;
}