[模板] 杜教筛

引言

在竞赛中我们经常遇到这样一类求函数的前缀和问题，即

$$
\boldsymbol{S}(n) = \sum_{i=1}^n \boldsymbol{f}(x)
$$

其中 $\boldsymbol{f}(x)$ 一般是积性函数，通常来说普通的线性筛法已经足够，但是仍可以以更快的速度求解

杜教筛

在学习杜教筛之前需要了解 狄利克雷卷积与莫比乌斯反演

我们考虑构造一个数论函数 $\boldsymbol{g}(x)$

$$
\begin{aligned}
\sum_{i=1}^n (\boldsymbol{f} * \boldsymbol{g})(i) & = \sum_{i=1}^n \sum_{xy = i}\boldsymbol{f}(x) \boldsymbol{g}(y) \\
&= \sum_{y=1}^n \boldsymbol{g}(y) \sum_{xy \leq n} \boldsymbol{f}(x) \\
&= \sum_{y=1}^n \boldsymbol{g}(y) \sum_{x=1}^{\lfloor \frac{n}{x}\rfloor} \boldsymbol{f}(x) \\
&= \sum_{y=1}^n \boldsymbol{g}(y) \boldsymbol{S}(\lfloor \frac{n}{x} \rfloor)
\end{aligned}
$$

移项可得
$$
\boldsymbol{g}(1) \boldsymbol{S}(n) = \sum_{i=1}^n (\boldsymbol{f} * \boldsymbol{g})(i) – \sum_{y=2}^n \boldsymbol{g}(y) \boldsymbol{S}(\lfloor \frac{n}{x} \rfloor)
$$

我们发现式子是一个递归的形式，可以将 $S(n)$ 递归求得

我们只需要使得 $\boldsymbol{f} * \boldsymbol{g}$ 的前缀和能够快速求出，即可使用杜教筛

同时我们可以将 $\boldsymbol{S}(N)$ 在较小的可计算的范围 $N \leq n^{\frac{2}{3}}$ 内进行预处理，剩下 $\frac{1}{3}$ 使用 $\mathrm{map}$ 或者数组保存，复杂度为 $O(n^{\frac{2}{3}})$，使用 $\mathrm{map}$ 则会多一个 $\log$，可以使用 $\mathrm{unordered_map}$，但是常数依旧巨大

我们可以使用 sum2[x] 来表示在 $\frac{N}{x}$ 时的答案

求欧拉函数的前缀和

因为
$$
\boldsymbol{\varphi} * \boldsymbol{1} = \boldsymbol{id}
$$

那么我们令 $\boldsymbol{g}(n) = 1$，即可快速求出 $(\boldsymbol{f} * \boldsymbol{g})(n) = n$

后半部分使用整除分块即可求解

inline long long getPhi(long long n){
    if(n < SIZ){
        return sumA[n];
    }

    int x = N / n;

    if(visA2[x])
        return sumA2[x];

    visA2[x] = true;
    long long &ans = sumA2[x];

    ans = n * (n + 1) / 2;

    int r;
    for(register int l=2; l<=n; l = r + 1){
        r = n / (n / l);
        ans -= (r - l + 1) * getPhi(n / l);
    }

    return ans;
}

求莫比乌斯函数的前缀和

由

$$
\boldsymbol{\mu} * \boldsymbol{1} = \boldsymbol {\epsilon}
$$

我们令 $\boldsymbol{g}(n) = 1$，也可以快速求出 $(\boldsymbol{f} * \boldsymbol{g})(n)$ 的值

inline int getMu(int n){
    if(n < SIZ)
        return sumB[n];

    int x = N / n;

    if(visB2[x])
        return sumB2[x];

    visB2[x] = true;
    long long &ans = sumB2[x];

    ans = 1;

    int r;
    for(register int l=2; l<=n; l = r + 1){
        r = n / (n / l);
        ans -= (r - l + 1) * getMu(n / l);
    }

    return ans; 
}