杜教筛学习笔记

杜教筛

杜教筛的基本形式

对于积性函数$g(n)$我们希望求他的前缀和$S_g(n)$，如果有另一积性函数$f(n)$满足$f\ast g=h$，且$fh$的前缀和易求，那么我们可以通过$S_f(n)S_h(n)$快速的求出$S_g(n)$。

$$\begin{array}{rl}{S}_h(n)& =\sum _{i=1}^n\sum _{d|i}f(d)\cdot g(\frac{n}{d})\\ & =\sum _{d=1}^n\sum _{i=1}^{\lfloor \frac{n}{d}\rfloor }f(d)\cdot g(i)\\ & =\sum _{d=1}^{n}f(d)\cdot S_g\left(\lfloor \frac{n}{d}\rfloor \right)\\ S_g(n)& =S_h(n)-\sum _{d=2}^{n}f(d)\cdot S_g\left(\lfloor \frac{n}{d}\rfloor \right)\end{array}$$

对后面的$S_g\left(\lfloor \frac{n}{d}\rfloor \right)$进行数论分块就可以更快速的求出$g$的前缀和。

时间复杂度分析：

在求解过程中我们只需要关注$S_g$在$\{\lfloor \frac{n}{k}\rfloor |k\subseteq {\mathbb{N}}^{\ast }\}$处的取值，这部分一共只有$O\left(\sqrt{n}\right)$种取值，因此整个时间复杂度是：

$$\begin{array}{rl}& \sum _{i=1}^{\sqrt{n}}O(\sqrt{i})+\sum _{i=1}^{\sqrt{n}}O(\sqrt{\frac{n}{i}})\ge {\int }_0^{\sqrt{n}}(\sqrt{x}+\sqrt{\frac{n}{x}})dx=O\left(n^{\frac{3}{4}}\right)\\ & \sum _{i=1}^{\sqrt{n}}O(\sqrt{i})+\sum _{i=1}^{\sqrt{n}}O(\sqrt{\frac{n}{i}})\le \sqrt{n}\cdot 2\cdot \sqrt{\sqrt{i}\cdot \sqrt{\frac{n}{i}}}=O\left(n^{\frac{3}{4}}\right)\end{array}$$

通过上面的分析我们可以知道杜教筛的时间复杂度是$O\left(n^{\frac{3}{4}}\right)$

如果能预处理出前$n^{\frac{2}{3}}$项的前缀和，复杂度可以优化到$O\left(n^{\frac{2}{3}}\right)$

一些例子：

• $\mu$的前缀和，有$\mu \ast I=e$
  所以，$S_{\mu }(n)=1-\sum _{i=2}^n{S}_{\mu }\left(\lfloor \frac{n}{i}\rfloor \right)$

• $\phi$的前缀和，有$\phi \ast I=id$
  所以，$S_{\phi }(n)=\frac{n\cdot (n+1)}{2}-\sum _{i=2}^n{S}_{\phi }\left(\lfloor \frac{n}{i}\rfloor \right)$

贝尔级数

我们发现如果要用杜教筛求解积性函数前缀和问题，如何构造合适的狄利克雷卷积式子是很重要的一环，贝尔级数就是一个强有力的构造工具。

对于积性函数我们重点关注其在质数的幂处的取值，积性函数$f(n)$在质数$p$处的贝尔级数定义为$F_p(x)=\sum _{i=0}^{\mathrm{\infty }}f(p^i)x^i$

然后数论函数的狄利克雷卷积就可以变成一般多项式的卷积：

$$h=f\ast g\iff H_p(x)=\sum _{i=0}^{\mathrm{\infty }}h(p^i)x^i=\sum _{i=0}^{\mathrm{\infty }}\sum _{k=0}^{i}f(p^k)\cdot x^k\cdot g(p^{i-k})\cdot x^{i-k}=G_p(x)\ast F_p(x)$$

一些常见积性函数的贝尔级数：

• $e\to 1$
• $I\to \sum _{i=0}^{\mathrm{\infty }}{x}^i=\frac{1}{1-x}$
• $id\to \sum _{i=0}^{\mathrm{\infty }}{p}^i\cdot x^i=\frac{1}{1-p\cdot x}$
• $i{d}_k\to \sum _{i=0}^{\mathrm{\infty }}{p}^{k\cdot i}\cdot x^i=\frac{1}{1-p^k\cdot x}$
• $\mu \to \sum _{i=0}^{\mathrm{\infty }}\mu \left(p^i\right)\cdot x^i=1-x$
• $\phi \to 1+\sum _{i=1}^{\mathrm{\infty }}{p}^{i-1}\cdot (p-1)\cdot x^i=1+\frac{p-1}{p}\cdot \sum _{i=1}^{\mathrm{\infty }}{p}^i\cdot x^i=\frac{1-x}{1-p\cdot x}$

一些例子：

• 求积性函数$f(p^k)=p^k+[k>0]\cdot (-1{)}^k$的前缀和。
  $$\begin{array}{rl}& \mathrm{先}\mathrm{写}\mathrm{出}f\mathrm{的}\mathrm{贝}\mathrm{尔}\mathrm{级}\mathrm{数}{F}_p(x)=1+\sum _{i=1}^{\mathrm{\infty }}(p^i+(-1{)}^i)\cdot x^i=\frac{1+p\cdot x^2}{(1-p\cdot x)\cdot (1+x)}\\ & F_p(x)\ast (1-p\cdot x)\cdot (1+x)=1+p\cdot x^2,\mathrm{其}\mathrm{中}(1-p\cdot x)\mathrm{对}\mathrm{应}\mu \cdot id(\mathrm{点}\mathrm{乘})(1+x)\mathrm{对}\mathrm{应}\mu \cdot \mu \\ & 1+p\cdot x^2\mathrm{对}\mathrm{应}\mathrm{的}\mathrm{函}\mathrm{数}\mathrm{满}\mathrm{足}h(n)=[\mathrm{\exists }d,n=d^2\wedge {\mu }^2(\sqrt{n})=1]\cdot \sqrt{n}\\ & \mathrm{令}g(n)=(\mu \cdot \mu )\ast (\mu \cdot id),\mathrm{则}{S}_f(n)=S_h(n)-\sum _{d=2}^{n}g(d)\cdot S_f\left(\lfloor \frac{n}{d}\rfloor \right)\\ & \sum _{i=1}^{n}h(i)=\sum _{i=1}^{\sqrt{n}}i\cdot {\mu }^2(i)\\ & \mathrm{现}\mathrm{在}\mathrm{需}\mathrm{要}\mathrm{求}g(n)\mathrm{的}\mathrm{前}\mathrm{缀}\mathrm{和}，\mathrm{注}\mathrm{意}\mathrm{到}(1-p\cdot x)\cdot (1+x)\cdot \frac{1}{1-px}=1+x\\ & S_g(n)=S_{\mu \cdot \mu }(n)-\sum _{d=2}^n\frac{d\cdot (d+1)}{2}\cdot S_g(\lfloor \frac{n}{d}\rfloor )\\ & \mathrm{然}\mathrm{后}\mathrm{需}\mathrm{要}\mathrm{求}\mu \cdot \mu \mathrm{的}\mathrm{前}\mathrm{缀}\mathrm{和}\\ & (\mu \cdot \mu )(n)=[n\mathrm{的}\mathrm{最}\mathrm{大}\mathrm{平}\mathrm{方}\mathrm{因}\mathrm{子}=1]=\sum _{d|n\mathrm{的}\mathrm{最}\mathrm{大}\mathrm{平}\mathrm{方}\mathrm{因}\mathrm{子}}\mu (d)=\sum _{d^2|n}\mu (d)\\ & S_{\mu \cdot \mu }(n)=\sum _{i=1}^n\sum _{d^2|i}\mu (d)=\sum _{d=1}^{\sqrt{n}}\mu (d)\cdot \lfloor \frac{n}{d^2}\rfloor \end{array}$$

Powerful Number

PN的定义是不含有次数为$1$的质因子的正整数。即考虑正整数$n$的标准分解形式$n=\prod _{i=1}^k{p}_i^{a_i},\mathrm{\forall }i{a}_i\ge 2$

这种形式的正整数一定可以写成$a^2\cdot b^3$，所以我们可以大概估计一下PN的数量，在${\int }_1^{\sqrt{n}}\sqrt[3]{\frac{n}{x^2}}dx=O\left(\sqrt{n}\right)$级别。这样的性质会给求一些积性函数的前缀和带来便利。

假设我们要求积性函数$f(n)$的前缀和，我们可以先构造一个新的积性函数$g(n)$满足$\mathrm{\forall }p\subseteq prime,f(p)=g(p)$。

此时必然存在一个积性函数$h(n)$满足$f=g\ast h$对于$\mathrm{\forall }p\subseteq prime,f(p)=g(1)\cdot h(p)+g(p)\cdot h(1)=h(p)+g(p)$

而根据$g$的构造我们可以得出$\mathrm{\forall }p\subseteq prime,h(p)=0$也就是说$h$只在PN处对答案有贡献，而PN的数量只有$O\left(\sqrt{n}\right)$!

$$S_f(n)=\sum _{i=1}^n\sum _{d|i}g(d)\cdot h(\frac{i}{d})=\sum _{d=1}^{n}h(d)\cdot S_g(\lfloor \frac{n}{d}\rfloor )=\sum _{d\subseteq PN}h(d)\cdot S_g(\lfloor \frac{n}{d}\rfloor )$$

接下来只需要求解$h$在PN处的取值，$S_g$可以用杜教筛计算。

考虑积性函数的取值只需要考虑其在质数的幂次处的取值：

$$f(p^k)=\sum _{i=0}^{k}g(p^i)\cdot h(p^{k-i})\Rightarrow h(p^k)=f(p^k)-\sum _{i=1}^{k}g(p^i)\cdot h(p^{k-i})$$

由于是PN这里只需要考虑小于$\sqrt{n}$的质数，这部分暴力计算即可。

例子：定义积性函数$f(p^k)=p^k\cdot (p^k-1)p\subseteq prime$求$S_f(n)$。

构造$g(n)=n\cdot \phi (n)$有$g(p)=p\cdot (p-1)=f(p)p\subseteq prime$

写出$g$的贝尔级数$1+\sum _{i=1}^{\mathrm{\infty }}{p}^i\cdot p^{i-1}\cdot (p-1)\cdot x^i=\frac{1-p\cdot x}{1-p^2\cdot x}$

注意到$g\ast id=i{d}_2$,所以$S_g(n)=\frac{n\cdot (n+1)\cdot (2\cdot n+1)}{6}-\sum _{d=2}^n\frac{n\cdot (n+1)}{2}\cdot S_g(\lfloor \frac{n}{d}\rfloor )$

对于$h$的PN部分暴力计算。

code:

#include <cstdio>
#include <cstring>
#include <iostream>
#include <vector>
#define int long long

using namespace std;

const int mo = 1e9 + 7;
const int inv2 = (mo + 1) >> 1;
const int inv3 = (mo + 1) / 3;
const int inv6 = 1ll * inv2 * inv3 % mo;

int prime[200005], cnt = 0, g[2000005], n, sum[100005], ans = 0;
bool a[2000005];
vector<int> h[10005];

int mi(int o, int p) {
    int yu = 1;
    while (p) {
        if (p & 1) yu = 1ll * yu * o % mo;
        o = 1ll * o * o % mo;
        p >>= 1;
    }
    return yu;
}

int S(int N) { return 1ll * N % mo * (N % mo + 1) % mo * inv2 % mo; }

int cal(int N) {
    if (N <= 2000000) return g[N];
    if (sum[n / N] != -1) return sum[n / N];
    int yu = 1ll * N % mo * (N % mo + 1) % mo * ((N + N) % mo + 1) % mo;
    yu = 1ll * yu * inv6 % mo;
    for (int l = 2, r; l <= N; l = r + 1) {
        r = N / (N / l);
        yu -= 1ll * (S(r) - S(l - 1)) * cal(N / l) % mo;
        yu = (yu % mo + mo) % mo;
    }
    return sum[n / N] = yu;
}

void dfs(int o, int p, int q) {
    ans = (ans + 1ll * q * cal(n / o) % mo) % mo;
    for (int tt = p; tt <= cnt; tt++) {
        if (n / o < prime[tt] * prime[tt]) break;
        int yu = o * prime[tt] * prime[tt];
        for (int t = 2; yu <= n; yu *= prime[tt], t++) {
            dfs(yu, tt + 1, 1ll * q * h[tt][t] % mo);
        }
    }
}

signed main() {
    memset(sum, -1, sizeof(sum));
    for (int i = 2; i <= 2000000; i++) {
        if (!a[i]) {
            prime[++cnt] = i;
            g[i] = 1ll * i * (i - 1) % mo;
        }
        for (int j = 1; j <= cnt && i * prime[j] <= 2000000; j++) {
            a[i * prime[j]] = 1;
            if (i % prime[j] == 0) {
                g[i * prime[j]] = 1ll * g[i] * prime[j] % mo * prime[j] % mo;
                break;
            }
            g[i * prime[j]] = 1ll * g[i] * prime[j] % mo * (prime[j] - 1) % mo;
        }
    }
    g[1] = 1;
    for (int i = 2; i <= 2000000; i++) {
        g[i] = (g[i - 1] + g[i]) % mo;
    }
    scanf("%lld", &n);
    for (int i = 1; i <= cnt; i++) {
        if (prime[i] > 100000) break;
        h[i].push_back(1);
        for (int j = prime[i], k = 1; j <= n; j *= prime[i], k++) {
            int yu = 1ll * j % mo * (j % mo - 1) % mo;
            for (int s = 0; s < k; s++) {
                yu -= 1ll * h[i][s] * mi(prime[i], 2 * k - 2 * s - 1) % mo *
                      (prime[i] - 1) % mo;
                yu = (yu % mo + mo) % mo;
            }
            h[i].push_back(yu);
        }
    }
    dfs(1, 1, 1);
    cout << ans << endl;
    return 0;
}