杜教筛和狄利克雷卷积

零、前置知识

1. 积性函数

积性函数的定义：若 $(a,b)=1$，则 $f(a\cdot b)=f(a)\cdot f(b)$。

常见的积性函数有：$\phi$ 函数，$\mu$
函数等。

积性函数有以下性质：

若 $f(x),g(x)$ 均为积性函数，则 $h(x)=f(x)\cdot g(x)$ 也为积性函数。

一、介绍

1.狄利克雷卷积

定义 $f\ast g(n)$ （读作：$f$ 卷 $g$，$n$ 为大小）为：$\sum _{i|n}f(i)\cdot g\left(\frac{n}{i}\right)$。

狄利克雷卷积有着交换律、结合律，请读者自证。

下面是一些常见函数：

$\phi ,\mu ,id,I,\epsilon$

其中 $id(n)=n$，$I(n)=1$，$\epsilon (n)$ 在 $n=1$ 时为 $1$，否则为 $0$。

现在来计算一下：$\phi \ast I(n)$

下面为计算过程：

φ ∗ I ( n ) = ∑ i ∣ n φ ( i ) = n = i d

===​φ∗I(n)i∣n∑​φ(i)nid​

还有一个常见的卷积：$\mu \ast I(n)=\epsilon$

下面为证明：

在 $n=1$ 时，显然成立

不妨设：$n={\alpha }_1^{p_1}\cdot {\alpha }_2^{p_2}\cdots \cdot {\alpha }_k^{p_k}$

$LHS=\sum _{i|n}\mu (i)$

由于 $\mu (i)$ 只有在不包含平方因子的时候才不是 $0$，所以：

L H S = ( n 0 ) − ( n 1 ) + ( n 2 ) − ⋯ = ∑ i = 0 k ( − 1 ) i ⋅ 1 k − i ⋅ ( n i ) = 0

===​LHS(0n​)−(1n​)+(2n​)−⋯i=0∑k​(−1)i⋅1k−i⋅(in​)0​

与 $\epsilon$ 的定义正好符合。

证毕。

相信你已熟知狄利克雷卷积，接下来我们进入杜教筛部分。

2.杜教筛

杜教筛是在高效的复杂度内求出一个积性函数的前缀和。

不妨设该函数为 $f(n)$，其前缀和为 $S(n)$，我们构造了另一个积性函数 $g(n)$。

接下来我们推一个式子：

∑ i = 1 n g ∗ f ( i ) = ∑ i = 1 n ∑ j ∣ i g ( i ) ⋅ f ( n i ) = ∑ j = 1 n ∑ i = 1 ⌊ n j ⌋ g [ j ] ⋅ f ( i ) = ∑ i = 1 n g ( i ) ⋅ S ( ⌊ n i ⌋ )

===​i=1∑n​g∗f(i)i=1∑n​j∣i∑​g(i)⋅f(in​)j=1∑n​i=1∑⌊jn​⌋​g[j]⋅f(i)i=1∑n​g(i)⋅S(⌊in​⌋)​

接下来考虑 $g(1)\cdot S(n)$，即得到：

g ( 1 ) ⋅ S ( n ) = ∑ i = 1 n g ( i ) ⋅ S ( ⌊ n i ⌋ ) − ∑ i = 2 n g ( i ) ⋅ S ( ⌊ n i ⌋ ) = ∑ i = 1 n g ∗ f ( i ) − ∑ i = 2 n g ( i ) ⋅ S ( ⌊ n i ⌋ )

==​g(1)⋅S(n)i=1∑n​g(i)⋅S(⌊in​⌋)−i=2∑n​g(i)⋅S(⌊in​⌋)i=1∑n​g∗f(i)−i=2∑n​g(i)⋅S(⌊in​⌋)​

所以，如果我们构造的 $g$ 满足如下条件：

• 可以快速得到 $g(n)$
• 可以快速得到 $g\ast f$ 的前缀和

然后注意到后面所减去的 $S(\lfloor \frac{n}{i}\rfloor )$ 是可以数论分块的，可以直接递归求得，这样复杂度为 $n^{\frac{3}{4}}$。

考虑到当我们要求的 $S(x)$ 小于阙值时，可以通过线性筛预处理出，大于等于阙值时，通过上述方式递归即可。

当阙值取 $n^{\frac{2}{3}}$ 时，复杂度最优，为 $n^{\frac{2}{3}}$。

二、习题

1.模板杜教筛

本题即求 $\phi$ 和 $\mu$ 的前缀和。

下面介绍求 $\phi$ 的前缀和，$\mu$ 同理，此处不分析。

我们取 $g(n)=I(n)$，上面已证：$\phi \ast I=id$

所以，可得：

S ( n ) = ∑ i = 1 n g ∗ f ( i ) − ∑ i = 2 n g ( i ) ⋅ S ( ⌊ n i ⌋ ) = ∑ i = 1 n i − ∑ i = 2 n S ( ⌊ n i ⌋ ) = n ⋅ ( n + 1 ) 2 − ∑ i = 2 n S ( ⌊ n i ⌋ )

===​S(n)i=1∑n​g∗f(i)−i=2∑n​g(i)⋅S(⌊in​⌋)i=1∑n​i−i=2∑n​S(⌊in​⌋)2n⋅(n+1)​−i=2∑n​S(⌊in​⌋)​

AC code：

#include 
using namespace std;
#define int long long
const int N = 3000300;
int T, n, B;
int a[N], cnt;
bool b[N];
int mu[N], phi[N], summu[N], sumphi[N];
unordered_map<int, int> smu, sphi;
void Get_Prime(int n) {
    phi[1] = 1;
    mu[1] = 1;
    for (int i = 2; i <= n; ++i) {
        if (!b[i]) {
            a[++cnt] = i;
            phi[i] = i - 1;
            mu[i] = -1;
        }
        for (int j = 1; j <= cnt && a[j] * i <= n; ++j) {
            b[a[j] * i] = 1;
            phi[a[j] * i] = phi[i] * (a[j] - 1);
            mu[a[j] * i] = -mu[i];
            if (i % a[j] == 0) {
                phi[a[j] * i] = phi[i] * a[j];
                mu[a[j] * i] = 0;
                break;
            }
        }
    }
    for (int i = 1; i <= n; ++i)
        summu[i] = summu[i - 1] + mu[i], sumphi[i] = sumphi[i - 1] + phi[i];
    return;
}
int GetSphi(int n) {
    if (n <= B)
        return sumphi[n];
    if (sphi.count(n) > 0)
        return sphi[n];
    int ans = n * (n + 1) / 2;
    for (int l = 2, r; l <= n;) {
        int val = n / l;
        r = n / val;
        ans -= (r - l + 1) * GetSphi(val);
        l = r + 1;
    }
    sphi[n] = ans;
    return ans;
}
int GetSmu(int n) {
    if (n <= B)
        return summu[n];
    if (smu.count(n) > 0)
        return smu[n];
    int ans = 1;
    for (int l = 2, r; l <= n;) {
        int val = n / l;
        r = n / val;
        ans -= (r - l + 1) * GetSmu(val);
        l = r + 1;
    }
    smu[n] = ans;
    return ans;
}
signed main() {
    scanf("%lld", &T);
    Get_Prime(N - 1);
    for (int __ = 1; __ <= T; ++__) {
        scanf("%lld", &n);
        B = pow(n, 2.0 / 3);
        printf("%lld %lld\n", GetSphi(n), GetSmu(n));
    }
    return 0;
}
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

2.四元组计数

题意

给定一个 $n$，求出有多少个整数 $a,b,c,d\in [1,n]$，满足：$a\cdot b=c\cdot d$。

数据范围：$1\le n\le 10^{11}$。

思路

注意到原式等价于：$\frac{a}{c}=\frac{b}{d}$，考虑将这一类分数约分后归为一类：$\frac{p}{q}(p,q)=1$。那么方案数：

∑ p = 1 n ∑ q = 1 n [ ( p , q ) = 1 ] ⋅ n max ⁡ ( p , q ) 2 = 2 ⋅ ∑ p = 1 n ∑ q = 1 p [ ( p , q ) = 1 ] ⋅ ⌊ n p ⌋ 2 − n 2 = 2 ⋅ ∑ p = 1 n ⌊ n p ⌋ 2 ⋅ φ ( p ) − n 2

==​p=1∑n​q=1∑n​[(p,q)=1]⋅max(p,q)n​22⋅p=1∑n​q=1∑p​[(p,q)=1]⋅⌊pn​⌋2−n22⋅p=1∑n​⌊pn​⌋2⋅φ(p)−n2​

我们记 $k={\sum }_{p=1}^n{\lfloor \frac{n}{p}\rfloor }^2\cdot \phi (p)$，则答案为 $2\cdot k-n^2$

k = ∑ i = 1 n φ ( i ) ⋅ ⌊ n i ⌋ 2 = ∑ i = 1 n φ ( i ) ∑ j = 1 ⌊ n i ⌋ ( 2 j − 1 ) = 2 ⋅ ∑ i ⋅ j ≤ n φ ( i ) ⋅ j − ∑ i ⋅ j ≤ n φ ( i ) = 2 ⋅ ∑ T = 1 n ∑ d ∣ T φ ( d ) ⋅ T d − ∑ T = 1 n ∑ d ∣ T φ ( d ) = 2 ⋅ ∑ i = 1 n φ ∗ i d ( i ) − ∑ i = 1 n i

k=====​i=1∑n​φ(i)⋅⌊in​⌋2i=1∑n​φ(i)j=1∑⌊in​⌋​(2j−1)2⋅i⋅j≤n∑​φ(i)⋅j−i⋅j≤n∑​φ(i)2⋅T=1∑n​d∣T∑​φ(d)⋅dT​−T=1∑n​d∣T∑​φ(d)2⋅i=1∑n​φ∗id(i)−i=1∑n​i​

现在已经很明朗了。

我们再记 $ans={\sum }_{i=1}^n\phi \ast id(i)$，在 $k=2\cdot ans-\frac{n\ast (n+1)}{2}$。

那么 $ans$ 也就是积性函数 $f(n)=\phi \ast id(n)$ 的前缀和。

我们考虑给他卷上一个 $I$，注意到狄利克雷卷积具有交换律，所以其等价于 $\phi \ast I\ast id(n)$，即 $id\ast id(n)={\sum }_{d|n}d\cdot \frac{n}{d}=n$，所以它的值就是 $n\cdot p(n)$，$p(n)$ 为 $n$ 的约数个数。考虑如何快速求出 $id\ast id(n)$，继续推式子。

∑ i = 1 n i d ∗ i d ( i ) = ∑ i = 1 n ∑ j = 1 ⌊ n i ⌋ i ⋅ j = ∑ i = 1 n i ⌊ n i ⌋ ⋅ ( ⌊ n i ⌋ + 1 ) 2

==​i=1∑n​id∗id(i)i=1∑n​j=1∑⌊in​⌋​i⋅ji=1∑n​i2⌊in​⌋⋅(⌊in​⌋+1)​​

这个可以通过整除分块，根号复杂度内求得，剩下的就是直接杜教筛即可。

式子：

$$su{m}_n=\sum _{i=1}^{n}id\ast id(i)-\sum _{i=2}^{n}su{m}_n/i$$

前半部分和后半部分皆用整除分块，复杂度还是 $\theta (n^{\frac{2}{3}})$。

Warning：此题数据范围有点大，需要 __int128。

AC code：

#include 
using namespace std;
#define int __int128
const int N = 4000010, M = 400010;
int T, B, n;
int a[M], cnt;
int v[N];
int pwr[N];
int tim[N];
int pid[N];
map<int, int> mp;
void Get_Prime(int n) {
    for (int i = 1; i <= n; ++i)
        v[i] = i;
    pid[1] = 1;
    for (int i = 2; i <= n; ++i) {
        if (v[i] == i) {
            a[++cnt] = i;
            pwr[i] = 1;
            tim[i] = i;
            pid[i] = 2 * i - 1;
        }
        for (int j = 1; j <= cnt && a[j] * i <= n; ++j) {
            v[a[j] * i] = a[j];
            pwr[a[j] * i] = 1;
            tim[a[j] * i] = a[j];
            pid[a[j] * i] = pid[a[j]] * pid[i];
            if (i % a[j] == 0) {
                tim[a[j] * i] = tim[i] * a[j];
                pwr[a[j] * i] = pwr[i] + 1;
                int val = tim[a[j] * i], pval = val + pwr[a[j] * i] * (v[a[j] * i] - 1) * val / v[a[j] * i];
                pid[a[j] * i] = pid[a[j] * i / val] * pval;
                break;
            }
        }
    }
    for (int i = 1; i <= n; ++i)
        pid[i] += pid[i - 1];
    return;
}
int ask(int n) {
    if (n <= B)
        return pid[n];
    if (mp.count(n) > 0)
        return mp[n];
    int ans = 0;
    for (int l = 1, r; l <= n;) {
        int val = n / l;
        r = n / val;
        ans += (r - l + 1) * (l + r) / 2 * val * (val + 1) / 2;
        l = r + 1;
    }
    for (int l = 2, r; l <= n;) {
        int val = n / l;
        r = n / val;
        ans -= (r - l + 1) * ask(val);
        l = r + 1;
    }
    mp[n] = ans;
    return ans;
}
__int128 readin() {
    __int128 x = 0, f = 1;
    char ch = getchar();
    while (ch < '0' || ch > '9') {
        if (ch == '-')
            f = -1;
        ch = getchar();
    }
    while (ch >= '0' && ch <= '9') {
        x = x * 10 + (ch ^ 48);
        ch = getchar();
    }
    return x * f;
}
void print(__int128 x) {
    if (x < 0)
        puts("-"), x = -x;
    if (!x) {
        return;
    }
    print(x / 10);
    putchar(x % 10 + '0');
    return;
}
signed main() {
    Get_Prime(N - 1);
    B = N - 1;
    T = readin();
    for (int __ = 1; __ <= T; ++__) {
        n = readin();
        int ans = ask(n);
        ans = ans * 2 - n * (n + 1) / 2;
        ans = ans * 2 - n * n;
        print(ans);
        puts("");
    }
    return 0;
}
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