杜教筛

大概是能求很多很多积性函数的前缀和的亚线性筛法。

假设我们要求的是 \(\sum\limits_{i = 1}^nf(i) = F(i)\)。

我们找一个积性函数 \(g\) 和 \(f\) 卷起来，看看会发生什么：

\[ \begin{aligned} (f * g)(n) &= \sum\limits_{d | n}f(d)g(\dfrac{n}{d}) \\ &= \sum\limits_{d | n}f(\dfrac{n}{d})g(d)\iff \\ \sum\limits_{i = 1}^n (f * g)(i) &= \sum\limits_{i = 1}^n \sum\limits_{d | i} f(\dfrac{i}{d})g(d) \\ &= \sum\limits_{d = 1}^n\sum\limits_{i = 1\land d | i}^{n} f(\dfrac i d)g(d) \\ &= \sum\limits_{d = 1}^ng(d)\sum\limits_{i = 1}^{\lfloor\frac n d\rfloor} f(i) \\ &= \sum\limits_{d = 1}^ng(d)F(\lfloor\dfrac n d\rfloor) \end{aligned} \]

如果我们能知道 \(f * g\) 的前缀和 \(g\) 的前缀和，似乎能快速求 \(F\)，但是直接做是不行的，我们考虑分割子问题：

\[ \begin{aligned} g(1)F(n) &= \sum\limits_{d = 1}^ng(d)F(\lfloor\dfrac n d\rfloor) -\sum\limits_{d = 2}^ng(d)F(\lfloor\dfrac n d\rfloor) \\ &= \sum\limits_{i = 1}^n (f * g)(i) - \sum\limits_{d = 2}g(d)F(\lfloor\dfrac n d\rfloor) \end{aligned} \]

也就是说如果我们知道 \((f * g)\) 和 \(g\) 的前缀和，我们可以通过整除分块来计算 \(F\)，式子当中的 \(F(\lfloor\dfrac n d\rfloor)\) 则可以递归处理。

伪代码大概是这样的：

\[ \begin{array}{ll} & \textbf{Du's Sieve method for} \\ & \textbf{Calculating prefix sum of a multiplicative function} \\ & \\ 1 & \textbf{Input.} \text{An integer } n\\ 2 & \textbf{Output.} \text{The sum of} \sum\limits_{i = 1}^n f(i)\\ 3 & \textbf{Function.}\ \text{Sieve}(n) \\ 4 & \qquad answer \gets \sum\limits_{i = 1}^n(f * g)(i) \\ 5 & \qquad l \gets 2, r\gets 0 \\ 6 & \qquad \textbf{while}\text{ l} \le n: \\ 7 & \qquad \qquad r \gets \min(n, \left\lfloor\dfrac{n}{\lfloor\dfrac{n}{l}\rfloor}\right\rfloor) \\ 8 & \qquad \qquad answer \gets answer - (\sum\limits_{i = l}^r g(i) \times \text{Sieve}(\lfloor\dfrac{n}{l}\rfloor)) \\ 9 & \qquad \textbf{return}\ answer / g(1) \end{array} \]

然后有几个常数上的小优化：

  1. 小于 \(B\) 的直接线性筛，这部分 \(O(B)\)。
  2. 记忆化，剩下没线性筛的部分直接开一个 map 或者 Hash Table 就行，但是这太慢了，我们不能接受，注意到这里求和会调用的部分全是 \(\lfloor\dfrac n i\rfloor\) 的形式，于是我们开一个数组记这些位置即可，也就是说，把 \(\lfloor\dfrac n i\rfloor\) 映射成 \(\left\lfloor\dfrac{n}{\lfloor\dfrac{n}{i}\rfloor}\right\rfloor\)，这部分的数组开个 \(O(n^{\frac 1 3})\) 差不多了，保险一点也可以开 \(O(\sqrt n)\)。

\(B\) 取 \(\sqrt n\) 不一定是最优的。

这个东西的复杂度是：\(\sum\limits_{i = 1}^{\lfloor\frac n B\rfloor}\sqrt{\dfrac n i} = O(\dfrac{n}{\sqrt B})\) 的。

由均值不等式：\(B = n^{\frac{2}{3}}\) 时，复杂度最优，为 \(O(n^{\frac 2 3})\)

一些小应用：

  1. 筛 \(\mu\)，注意到 \(\mu * Id_0 = \epsilon\)，这是 Trivial 的。
  2. 筛 \(\varphi\)，注意到 \(\varphi * Id_0 = Id_1\)，这是 Trivial 的。
上面两个玩意儿的 Code

Source: 杜教筛模板

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
// author : black_trees

#include <cmath>
#include <bitset>
#include <cstdio>
#include <cstring>
#include <iostream>
#include <algorithm>

#define endl '\n'
#define int long long

using namespace std;
using i64 = long long;

const int si = 5e6 + 10;

int n, m, t;
int smu[si], mu[si];
int sphi[si], phi[si];
int pri[si], vis[si];
void sieve(int n) {
    mu[1] = 1, phi[1] = 1, t = 0;
    memset(vis, 0, sizeof vis);
    for(int i = 2; i <= n; ++i) {
        if(!vis[i]) pri[++t] = i, vis[i] = i, mu[i] = -1, phi[i] = i - 1;
        for(int j = 1; j <= t; ++j) {
            if(pri[j] > vis[i] || i * pri[j] > n) break;
            vis[i * pri[j]] = pri[j];
            if(i % pri[j] == 0) { mu[i * pri[j]] = 0, phi[i * pri[j]] = phi[i] * pri[j];  break; }
            mu[i * pri[j]] = -mu[i], phi[i * pri[j]] = phi[i] * (pri[j] - 1);
        }
        mu[i] += mu[i - 1], phi[i] += phi[i - 1];
    }
}
int id(int x) { return n / x; }
std::bitset<si> vm, vp;
int Mobius(int n) {
    if(n <= m) return mu[n];
    if(vm[id(n)]) return smu[id(n)];
    int ret = 1; // epsilon
    int l = 2, r = 0;
    for(; l <= n; ) {
        r = min((n / (n / l)), n);
        ret -= (r - l + 1) * Mobius(n / l), l = r + 1;
    }
    vm[id(n)] = true;
    return smu[id(n)] = ret; // Id_0(1) = 1
} 
int Phi(int n) {
    if(n <= m) return phi[n];
    if(vp[id(n)]) return sphi[id(n)];
    int ret = n * (n + 1) / 2; // Id_1
    int l = 2, r = 0;
    for(; l <= n; ) {
        r = min((n / (n / l)), n);
        ret -= (r - l + 1) * Phi(n / l), l = r + 1;
    }
    vp[id(n)] = true;
    return sphi[id(n)] = ret;
}

signed main() {

    cin.tie(0) -> sync_with_stdio(false);
    cin.exceptions(cin.failbit | cin.badbit);

    int T; cin >> T;
    sieve(si - 10);
    while(T--) {
        vm.reset(), vp.reset();
        cin >> n, m = si - 10;
        cout << Phi(n) << " " << Mobius(n) << endl;
    }

    return 0;
}

强一点的应用：

\(\sum\limits_{i = 1}^n\varphi(i)\cdot i\)。

不难注意到 \(f(i) = \varphi(i) \cdot i\) 依旧是一个积性函数。

考虑构造一个 \(g\) 使得 \(f * g\) 的前缀和比较好算（最好是 \(O(1)\)）。

不妨令 \(g = Id_1\)。

\[ \begin{aligned} (f * Id_1)(n)&= \sum\limits_{d | n}f(d)Id_1(\dfrac n d) \\ &= \sum\limits_{d | n}\varphi(d) \cdot d \cdot \dfrac n d \\ &= n \end{aligned} \]

所以 \(\sum\limits_{i = 1}^n (f * g)(i) = n^2\)。

\(g\) 的前缀和是 Trivial 的，然后就可以做了。

这里选择 \(g = Id_1\) 的原因大概是，我们想尝试消去 \(\varphi(d)\) 或者 \(d\) 来快速求和，\(\varphi(d)\) 的话如果卷上 \(Id_0\) 是没用的，展不开，所以考虑消去 \(d\)，注意到 Dirichlet 卷积形式里有 \(\dfrac n d\)，所以构造一个 \(Id_1\) 卷上去就能消掉了。

最后更新: June 27, 2023