杜教筛

发表于 分类于 ，

概论:

杜教筛可以做到在亚线性的时间复杂度下求出部分数论函数的前缀和.

事实上感觉和” 筛” 好像没啥关系

实现:

假设需要求解的前缀和为 \(S_f(n) = \sum_{i=1}^n f(i)\), 引入一个函数 \(g(n)\) 并考虑如下式子:

\[ \begin{aligned} & \sum_{i=1}^n (f*g)(i) \\ = & \sum_{i=1}^n \sum_{d \mid i} g(d) f(\frac{i}{d}) \\ = & \sum_{d=1}^n g(d) \sum_{i=1}^{\lfloor n/d \rfloor} f(i) \\ = & \sum_{i=1}^n g(i)S_f(\left\lfloor \frac{n}{i} \right\rfloor) \end{aligned} \]

将 \(g(1)S_f(n)\) 分离出来，便有:

\[ g(1)S_f(n) = \sum_{i=1}^n (f*g)(i) - \sum_{i=2}^{n}g(i)S_f(\left\lfloor \frac{n}{i} \right\rfloor) \]

(注意后面这个和式是从 2 开始的)

这个等式便是杜教筛的核心式。若我们能找到函数 \(g\), 使得:

  • \(\sum_{i=1}^n (f*g)(i)\) 易求
  • \(g(i)\) 的区间和易求

则可以使用数论分块加速求解 \(S_f(n)\)

一些例子:

  • \(\mu\)

已知 \(\mu * I = \epsilon\), 取 \(g(n) = I\), 有:

\[ S_{\mu}(n) = 1 - \sum_{i=2}^n S_{\mu}(\left\lfloor \frac{n}{i} \right\rfloor) \]

  • \(\varphi\)

已知 \(\varphi * I = id\), 取 \(g(n) = I\), 有:

\[ S_{\varphi}(n) = \frac{n(n+1)}{2} - \sum_{i=2}^n S_{\varphi}(\left\lfloor \frac{n}{i} \right\rfloor) \]

  • \(id \cdot \varphi\)

根据点积对狄利克雷卷积的分配律，有:

\[ \begin{aligned} (id \cdot \varphi)*id & = (id\cdot\varphi) * (id\cdot I)\\ & = id\cdot (\varphi*I) = id_2 \end{aligned} \]

取 \(g(n)=id\), 有:

\[ S_{id\cdot\varphi}(n) = \frac{n(n+1)(2n+1)}{6} - \sum_{i=2}^n iS_{id\cdot\varphi}(\left\lfloor \frac{n}{i} \right\rfloor) \]

复杂度分析:

若使用递归的方法实现杜教筛，暂时忽略前面的卷积和项， 设时间复杂度为 \(T(n)\), 则有:

\[ T(n) = \sum_{i=1}^{\lfloor \sqrt{n} \rfloor} T(i)+T(\left\lfloor \frac{n}{i} \right\rfloor) \]

同时对于 \(i \in \left[1,\lfloor \sqrt{n} \rfloor \right]\), 有 \(i < n/i\), 所以:

\[ T(n) = O\left(\sum_{i=1}^{\sqrt{n}} T\left(\frac{n}{i}\right)\right) \]

对于 \(T(n/i)\), 由于继续展开后均为高阶无穷小，所以只展开一层， 即取数论分块的时间复杂度 \(O(\sqrt{n/i})\):

\[\begin{aligned} T(n) & = O\left(\sum_{i=1}^{\sqrt{n}} \sqrt{\frac{n}{i}}\right) \\ & = O\left(\int_1^{\sqrt{n}} \sqrt{\frac{n}{x}}dx\right) \\ & = O(n^{3/4}) \end{aligned}\]

若线性筛出前 \(m\) 个数的前缀和， 则有:

\[\begin{aligned} T(n) & = O(m) + O\left(\sum_{i=1}^{n/m} \sqrt{\frac{n}{i}}\right) \\ & = O\left(m +\int_1^{n/m} \sqrt{\frac{n}{x}}dx\right) \\ & = O\left(m + \frac{n}{\sqrt{m}}\right) \end{aligned}\]

当 \(m = n/\sqrt{m}\) 时取得最小值， 此时 \(m=n^{2/3}\), 时间复杂度为 \(O(n^{2/3})\)

事实上从上面的分析也可以看出，省略的项比较多，所以常数稍大.

同时需要将计算 \(\sum (f\*g)\) 的时间复杂度考虑进去，若 \(\sum (f\*g)\) 的时间复杂度小于等于 \(O(n^{2/3})\) 则总体时间复杂度不变.

注意到计算过程中 \(S_f(\lfloor n/i \rfloor)\) 也将被计算，同时使用数论分块时计算 \(g\) 的区间和只需要计算全体 \(S_g(\lfloor n/i \rfloor)\) 即可 (可以参见数论分块来证明), 因此若事先使用杜教筛将 \(S_g(\lfloor n/i \rfloor)\) 算出并储存，总体时间复杂度也不变.

综上，使用杜教筛在 \(O(n^{2/3})\) 时间内计算 \(S_f\) 的条件如下:

  • \(f\) 能够使用线性筛计算得到.
  • 能够找到函数 \(g\), 使得:
    • 可使用杜教筛在 \(O(n^{2/3})\) 时间内计算 \(S_g\)
    • 计算 \(S_{f*g}\) 的时间复杂度小于等于 \(O(n^{2/3})\)

至于如何存储 \(S(\lfloor n/i \rfloor)\), 可以使用 unordered_map, 也可以使用根号分治的方法将 \(\lfloor n/i \rfloor\) 分为两类存储.

示例代码:

luogu4213 【模板】杜教筛（Sum）

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
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
#include <cstdio>
#include <cctype>
#include <vector>
#include <cstring>
#include <unordered_map>
using namespace std;
typedef long long LL;
const int Mn(1000500);

template<typename T>
void qrd(T& x) {
    x = 0;
    char c = getchar();
    while(!isdigit(c)) {
        c = getchar();
    }
    while(isdigit(c)) {
        x = x*10 + c-'0';
        c = getchar();
    }
}

/*线性筛部分*/
LL smu[Mn],sph[Mn];
bool isp[Mn];
vector<int> prm;
void init() {
    memset(isp,-1,sizeof isp);
    sph[1] = smu[1] = 1;
    isp[1] = false;
    for(int i(2);i<Mn;++i) {
        if(isp[i]) {
            prm.push_back(i);
            smu[i] = -1;
            sph[i] = i-1;
        }
        for(int a: prm) {
            if(1ll*a*i>=Mn) {
                break;
            }
            isp[i*a] = false;
            if(i%a==0) {
                smu[i*a] = 0;
                sph[i*a] = sph[i] * a;
                break;
            } else {
                smu[i*a] = -smu[i];
                sph[i*a] = sph[i] * (a-1);
            }
        }
        smu[i] += smu[i-1]; 
        sph[i] += sph[i-1];
    }
}

/*杜教筛 mu*/
unordered_map<int,LL> ansmu,ansph;
LL getMu(LL x) {
    if(x<Mn) {  //线性筛预处理优化
        return smu[x];
    }
    if(ansmu[x]) {
        return ansmu[x];
    }
    LL ret = 1;
    LL l(2),r;  //注意起始位置
    for(;l<=x;l = r+1) {
        r = x/(x/l);
        ret -= getMu(x/l)*(r-l+1);  //公式
    }
    return ansmu[x] = ret;  //存入unordered_map
}
/*杜教筛 phi*/
LL getPh(LL x) {
    if(x<Mn) {
        return sph[x];
    }
    if(ansph[x]) {
        return ansph[x];
    }
    LL ret = (1ll*x*(x+1)/2ll);
    LL l(2),r;
    for(;l<=x;l = r+1) {
        r = x/(x/l);
        ret -= getPh(x/l)*(r-l+1);  //公式too
    }
    return ansph[x] = ret;
}

int main() {
    init();
    int T;
    qrd(T);
    while(T--) {
        LL n; qrd(n);
        printf("%lld %lld\n",getPh(n),getMu(n));
    }
    return 0;
}