关于杜教筛

杜教筛是对某种函数快速求前缀和（不是所有，是固定参数的）的做法，经过预处理优化，一般能够优化到 $O(n^{\frac{2}{3}})$ 的时间复杂度。它由杜瑜皓引进，所以称为杜教筛。

原理

假设函数 $f(n),g(n)$ 是积性函数， $h=f\ast g$，则函数 $h(n)$ 是积性函数。

$F,G,H$ 分别是 $f,g,h$ 的前缀和。

此时我们给定 $n$ ，需要快速求出 $F(n)$ 。

因为
$$\begin{array}{rl}h(n)& =\sum _{d|n}f(d)\cdot g(\frac{n}{d})\\ & =f(n)\cdot g(1)+\sum _{d|n,d<n}f(d)\cdot g(\frac{n}{d})\end{array}$$

所以

$$\sum _{i=1}^{n}h(i)=\sum _{i=1}^{n}f(i)\cdot g(1)+\sum _{i=1}^n\sum _{d|n,d<i}f(d)\cdot g(\frac{i}{d})$$

此时，我们设 $j=i/d$ ，则有

$$\begin{array}{rl}\sum _{i=1}^{n}h(i)& =H(n)\\ & =\sum _{i=1}^{n}f(i)\cdot g(1)+\sum _{j=2}^{n}g(j)\sum _{d=1}^{n/j}f(d)\\ & =F(n)\cdot g(1)+\sum _{j=2}^{n}g(j)\cdot F(\frac{n}{j})\end{array}$$

因为 $g(1)=1$ ，则可以写为

$$F(n)=H(n)-\sum _{j=2}^{n}g(j)\cdot F(\frac{n}{j})$$

如果 $H(n)$ 和 $G(j)$ 可以快速求（比如 $O(1)$ 求），我们就得到了一个关于 $F(n)$ 的递推式。

此时用数论分块可以降到 $O(n^{\frac{3}{4}})$ 的时间复杂度。

如果我们将前 $n^{\frac{2}{3}}$ 个 $F$ 数组预处理出来，则时间可以降为 $O(n^{\frac{2}{3}})$

例题

例题一

求 ${\sum }_{i=1}^n\varphi (i)$ 和 ${\sum }_{i=1}^n\mu (i)$ 。

其中 $n\le 2^{31}-1$

思路

板题，简单说说。

因为我们知道 $\varphi \ast I=Id$ ，并且 $I,Id$ 的前缀和都可以 $O(1)$ 求，此时就有

$$\sum _{i=1}^n\varphi (i)=\frac{n\cdot (n+1)}{2}-\sum _{j=2}^{n}I(j)\sum _{i=1}^{n/j}\varphi (i)$$

此时设 $f(n)={\sum }_{i=1}^n\varphi (i)$ ，则有

$$f(n)=\frac{n\cdot (n+1)}{2}-\sum _{j=2}^{n}I(j)\cdot f(n/j)$$

因为 $\mu \ast I=ϵ$ ，所以同上推导就可以了。

这就是杜教筛的模板。

代码

#include <bits/stdc++.h>
using namespace std;
typedef long long LL;
const int MAXN = 5000000;
int T, n, ptot, prime[MAXN + 10], phi[MAXN + 10], mu[MAXN + 10];
bool bz[MAXN + 10];
LL premu[MAXN + 10], prephi[MAXN + 10];
map<LL, LL> mpphi, mpmu;
void get(int n) {
    phi[1] = mu[1] = 1;
    for (int i = 2; i <= n; i++) {
        if (!bz[i])
            ptot++, prime[ptot] = i, phi[i] = i - 1, mu[i] = -1;
        for (int j = 1; j <= ptot && prime[j] * i <= n; j++) {
            bz[prime[j] * i] = true;
            if (i % prime[j] == 0) {
                phi[i * prime[j]] = phi[i] * prime[j];
                mu[i * prime[j]] = 0;
                break;
            } else {
                phi[i * prime[j]] = phi[i] * (prime[j] - 1);
                mu[i * prime[j]] = -mu[i];
            }
        }
    }
    for (int i = 1; i <= n; i++) premu[i] = premu[i - 1] + mu[i], prephi[i] = prephi[i - 1] + phi[i];
}
LL dfsphi(int n) {
    if (n <= MAXN)
        return prephi[n];
    if (mpphi.count(n))//如果之前求过，就把值拿来用
        return mpphi[n];
    LL sum = 0;
    for (int i = 2; i <= n; i++) {
        int last = i - 1;
        i = n / (n / i);
        sum = sum + dfsphi(n / i) * (i - last);
        if (i == n)//因为i加1后有可能会爆int
            break;
    }
    return mpphi[n] = n * (n + 1ll) / 2 - sum;
}
LL dfsmu(int n) {
    if (n <= MAXN)
        return premu[n];
    if (mpmu.count(n))//如果之前求过，就把值拿来用
        return mpmu[n];
    LL sum = 0;
    for (int i = 2; i <= n; i++) {
        int last = i - 1;
        i = n / (n / i);
        sum = sum + dfsmu(n / i) * (i - last);
        if (i == n)//因为i加1后有可能会爆int
            break;
    }
    return mpmu[n] = 1 - sum;
}
int main() {
    scanf("%d", &T);
    get(MAXN);//预处理处前n个前缀和
    while (T--) {
        scanf("%d", &n);
        printf("%lld %lld\n", dfsphi(n), dfsmu(n));
    }
    return 0;
}

例题二

求 ${\sum }_{i=1}^n\varphi (i^2)$ 模 $10^9+7$ 的结果。

其中 $n\le 10^9$ 。

思路

我们要知道一个性质，即 $\varphi (i^2)=i\cdot \varphi (i)$ 。

所以原式就变成了 ${\sum }_{i=1}^{n}i\cdot \varphi (i)$ 。

因为 $n^2=n{\sum }_{i|n}phi(i)={\sum }_{i|n}\varphi (i)\cdot i\cdot \frac{n}{i}$ 。

所以 $n\cdot \varphi (n)=n^2-{\sum }_{i|n,i<n}\varphi (i)\cdot i\cdot \frac{n}{i}$ 。

此时我们设 $f(i)=\varphi (i^2)=i\cdot \varphi (i)$ ，则有 $f(n)=n^2-{\sum }_{i|n,i<n}f(i)\cdot \frac{n}{i}$ 。

我们设 $F(n)={\sum }_{i=1}^{n}f(i)$ ，则

$$\begin{array}{rl}F(n)& =\sum _{i=1}^{n}f(i)\\ & =\sum _{i=1}^n{i}^2-\sum _{i=1}^n\sum _{j|i,j<i}f(j)\cdot \frac{i}{j}\\ & =\sum _{i=1}^n{i}^2-\sum _{d=2}^{n}d\sum _{i=1}^{n/d}f(i)\\ & =\sum _{i=1}^n{i}^2-\sum _{d=2}^{n}d\cdot F(\frac{n}{d})\end{array}$$

此时就用杜教筛来做就行了。

代码

#include <bits/stdc++.h>
using namespace std;
typedef long long LL;
const int MAXN = 5000000;
const LL MOD = 1e9 + 7;
const LL rev = 166666668;
int T, n, ptot, prime[MAXN + 10];
bool bz[MAXN + 10];
LL prephi[MAXN + 10], phi[MAXN + 10];
map<LL, LL> mpphi;
void get(int n) {
    phi[1] = 1;
    for (int i = 2; i <= n; i++) {
        if (!bz[i])
            ptot++, prime[ptot] = i, phi[i] = 1ll * (i - 1) * i % MOD;
        for (int j = 1; j <= ptot && prime[j] * i <= n; j++) {
            bz[prime[j] * i] = true;
            if (i % prime[j] == 0) {
                phi[i * prime[j]] = phi[i] * prime[j] % MOD * prime[j] % MOD;
                break;
            } else {
                phi[i * prime[j]] = phi[i] * (prime[j] - 1) % MOD * prime[j] % MOD;
            }
        }
    }
    for (int i = 1; i <= n; i++) prephi[i] = (prephi[i - 1] + phi[i]) % MOD;
}
LL dfsphi(int n) {
    if (n <= MAXN)
        return prephi[n];
    if (mpphi.count(n))
        return mpphi[n];
    LL sum = 0;
    for (int i = 2; i <= n; i++) {
        int last = i - 1;
        i = n / (n / i);
        sum = (sum + (1ll * (i - last) * (i + last + 1) / 2) % MOD * dfsphi(n / i) % MOD) % MOD;
    }
    return mpphi[n] = (1ll * n * (n + 1ll) % MOD * (2ll * n + 1ll) % MOD * rev % MOD - sum + MOD) % MOD;
}
int main() {
    get(MAXN);//这里求的是phi(i)的平方的前缀和
    scanf("%d", &n);
    printf("%lld\n", dfsphi(n));
    return 0;
}