莫比乌斯反演
定义
数论函数：\(\mathbb{N}^* \to \mathbb{C}\) 。
积性函数：\(\forall \gcd(a,b) = 1, f(a)f(b) = f(ab)\) 的数论函数 \(f\) 。
由唯一分解定理可以写出：若 \(f\) 为积性函数，那么 \(f(n) = \prod\limits_{i = 1}^mf(p_i^{e_i})\)
完全积性函数：\(\forall a,b, f(a)f(b) = f(ab)\) 的数论函数 \(f\) 。
Dirichlet 卷积：两个数论函数 \(f, g\) 的 Dirichlet 卷积 \(h\) 满足：\((f * g)(n) = h(n) = \sum\limits_{d|n} f(d)g(\dfrac{n}{d})\)
满足交换律，结合律：\(f*g = g*f, f*(g*h) = (f*g)*h\) 。
满足对点值加法的分配律：\(f*(g+h) = f*g + f*h\) 。
若 \(f,g\) 为积性函数，那么 \(h = f*g\) 仍然是个积性函数。
存在乘法幺元函数 \(\epsilon(n) = [n = 1]\) （这里的中括号是 Iverson Bracket）满足 \(f*\epsilon = f\) 。
你可以简单的理解为，有一个数论函数集合 \(R\) ，在其上定义两个二元运算：点值加法和 Dirichlet 卷积，可以构成一个交换环 \((R;+,*)\) ，\(\epsilon\) 即为乘法幺元。
类似的，如果有 \(f*g = \epsilon\) ，那么 \(g\) 是 \(f\) 的 Dirichlet 逆元。
很多时候在知道 \(f,h\) 的情况下可以反推 \(g\) ，这就是我们等会要说的莫比乌斯反演。
常见的完全积性函数：
幺元函数 \(\epsilon(n) = [n = 1]\) 。
恒等函数 \(Id_k(n) = n^k\) ，一般来说就用 \(k = 0, k = 1\) 。
常见的积性函数：
约数函数 \(\sigma_k(n) = \sum\limits_{d | n} d^k = \prod\limits_i(\sum\limits_{j = 0}^{e_i}p_i^{jk})\) 。
\(\sigma_0\) 表示约数个数，\(\sigma_1\) 表示约数和，后面那个等号就是积性函数的性质的一个变形？欧拉函数 \(\varphi(n)\) ，这个在之前的博客 里提到过。
一个结论：\(\varphi * Id_0 = Id_1\) ，这个是显然的，写一下就是 \((\varphi * Id_0)(n) = \sum\limits_{d | n} \varphi(d) Id_0(\dfrac{n}{d}) = \sum\limits_{d | n} \varphi(d) = n = Id_1(n)\) 。
莫比乌斯函数 （重头戏）：\(\mu(n = \prod\limits_{i = 1}^k p_i^{e_i}) = \begin{cases}(-1)^k &\forall i, e_i = 1 \\ 0 &\text{otherwise.}\end{cases}\)
线性筛求积性函数
莫比乌斯函数：
考虑线性筛的过程。
如果进入第一层循环，发现 \(vis(i) = 0\) 即 \(\text{PRIME}(i)\) 时，\(\mu(i) = -1\) 。
然后在利用 \(pri(j)\) 标记 \(i \times pri(j)\) 时。
如果 \(pri(j) | i\) ，证明 \(i \times pri(j)\) 存在 \(e_i \ge 2\) ，\(\mu(i \times pri(j)) = 0\) 。
否则就是 \(\omega(i \times pri(j)) = \omega(i) + 1\) ，取反即可。
Code
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16 int n , m ; int mu [ si ]; int pri [ si ], vis [ si ]; void sieve ( int n ) { mu [ 1 ] = 1 , m = 0 ; memset ( vis , 0 , sizeof vis ); for ( int i = 2 ; i <= n ; ++ i ) { if ( ! vis [ i ]) pri [ ++ m ] = i , vis [ i ] = i , mu [ i ] = -1 ; for ( int j = 1 ; j <= m ; ++ 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 ; break ; } mu [ i * pri [ j ]] = - mu [ i ]; } } }
通用：
其实就是在标记的时候，考虑一下 \(i\) 和 \(pri(j)\) 的关系，并根据积性函数的性质递推即可。
比如约数个数 \(\sigma_0\) ：
如果 \(i \not\equiv 0 \pmod{pri(j)}\) ，那么按照积性函数的性质递推即可。
否则，展开式子，需要记录最小质因子的指数 \(e_i\) 。
Code
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19 int n , m ; int f [ si ], d [ si ]; int pri [ si ], vis [ si ]; void sieve ( int n ) { f [ 1 ] = 0 , d [ 1 ] = 1 , m = 0 ; memset ( vis , 0 , sizeof vis ); for ( int i = 2 ; i <= n ; ++ i ) { if ( ! vis [ i ]) pri [ ++ m ] = i , vis [ i ] = i , f [ i ] = 1 , d [ i ] = 2 ; for ( int j = 1 ; j <= m ; ++ j ) { if ( pri [ j ] > vis [ i ] || pri [ j ] * i > n ) break ; vis [ i * pri [ j ]] = pri [ j ]; if ( i % pri [ j ] == 0 ) {
f [ i * pri [ j ]] = f [ i ] + 1 ,
d [ i * pri [ j ]] = d [ i ] / ( f [ i ] + 1 ) * ( f [ i * pri [ j ]] + 1 ); break ;
} f [ i * pri [ j ]] = 1 , d [ i * pri [ j ]] = d [ i ] * d [ pri [ j ]]; } } }
另外有一个结论：\(\sum\limits_{i = 1}^n\sigma_0(i) = \sum\limits_{i = 1}^n\lfloor\dfrac n i\rfloor\) 。
这个在杜教筛或者是后面莫反的时候也许会有用。
欧拉函数：link
这个也是类似的，如果 \(i \equiv 0 \pmod{pri(j)}\) 的话，展开式子尝试找到关系即可。
莫比乌斯反演
重要性质
\(\mu * Id_0 = \epsilon\) ，即 \(Id_0\) 的 Dirichlet 逆元是 \(\mu\)
证明 ：
若 \(n = 1\) ，\((\mu * Id_0)(n) = 1\) ，这是显然的。
若 \(n \not= 1\) ，记 \(k = \omega(n) \iff \exists p\ s.t.\ \text{PRIME}(p),p^k | n, p^{k + 1}\not|\; n\) ，则有：
因为 \(\mu, Id_0\) 为积性函数，所以 \((\mu * Id_0)(n) = (\mu * Id_0)(p^k) \cdot (\mu * Id_0)(\dfrac{n}{p^k}) = (\sum\limits_{d | p^k} \mu(d)) \cdot (\sum\limits_{d | \frac{n}{p^k}} \mu(d)) = 0\) 。
（前者由于 \(d | p^k\) 的 \(d\) 都是 \(p^t\) 形式的，而 \(p\) 为素数，所以所有这样的 \(\mu(d)\) 都等于 \(0\) ，求和自然也是 \(0\) ）.
Q.E.D.
炫酷反演魔术
刚才我们提到了，\(Id_0\) 的 Dirichlet 逆元是 \(\mu\) 。
于是我们有这样的一个想法：
如果有一个函数 \(g = Id_0 * f\) ，即 \(g(n) = \sum\limits_{d | n} f(d)\) 。
如果知道 \(f\) 求 \(g\) 是容易的，但是如果知道 \(g\) 要求 \(f\) 呢？
换句话说我们要让 \(f\) 独立出来。
不妨考虑两边同时卷上一个 \(\mu\) ，构造 \(\epsilon\) 和 \(f\) 的 Dirichlet 卷积形式：
\[
g * \mu = (Id_0 * \mu) * f \iff f = g * \mu
\]
这个式子被叫做莫比乌斯反演公式，非常好用！
例题
我们来看一个简单的练手题：
Ex.1
求 \(\sum\limits_{i = 1}^n \sum\limits_{j = 1}^m [\gcd(i, j) = 1](m < n)\) 。
\(n, m \le 10^7\) 。
注意到这个 Iverson Bracket 神似 \(\epsilon(n)\) 的形式。
就是：\(\sum\limits_{i = 1}^n \sum\limits_{j = 1}^m \epsilon(\gcd(i, j))\) 。
因为 \(Id_0 * \mu = \epsilon\) ，不难想到使用 \(\mu\) 来替换：
\(\epsilon(\gcd(i, j)) = \sum\limits_{d | \gcd(i, j)} \mu(d)Id_0(\dfrac{\gcd(i, j)}{d}) = \sum\limits_{d | \gcd(i, j)} \mu(d)\) 。
原式：\(\sum\limits_{i = 1}^n \sum\limits_{j = 1}^m\sum\limits_{d | \gcd(i, j)} \mu(d)\)
这种一堆求和号的式子有一个很显然的技巧，尝试观察变量之间的关系，观察是否能交换 \(\sum\) 。
发现 \(d | i \land d | j\) ，并且 \(d \in [1,m]\) （\(\gcd(i, j)\) 怎么可能比 \(\min(i, j)\) 大？）（令 \(i = j\) 即可构造出所有 \([1, m]\) 中的 \(d\) ），所以原式：
\[
\sum\limits_{d = 1}^m\mu(d)\sum\limits_{d | i}\sum\limits_{d | j} 1
\]
后面直接写成整除的形式：
\[
\sum\limits_{d = 1}^m\mu(d)\lfloor\dfrac n d\rfloor\lfloor\dfrac m d\rfloor
\]
这部分就是个标准的二维整除分块，但是还需要求 \(\sum \mu\) ，这个需要筛一下。
如果 \(n, m\) 再大一点就要杜教筛了。
Tips: 另外一种理解方式是，考虑枚举 \(d\) ，然后加上 \(d = \gcd(i, j)\) 的条件，把 \(d\) 扔到最外层，然后把条件下放到内层循环里。
Ex.2
如果是 \([\gcd(i, j) = k]\) 呢？
多次询问，\(5e4\) 。
Source: Luogu3455 - Zap
这里设 \(n \le m\) ，我们尝试转化为 Ex.1 的形式：
\[
\begin{aligned}
&=\sum\limits_{i = 1}^n\sum\limits_{j = 1}^m[\gcd(i, j) = k]
\end{aligned}
\]
不妨从 \(i, j\) 中扣掉 \(k\) 这个因子，就是把 \(i, j\) 看作 \(p\cdot k, q \cdot k\) ，然后我们枚举 \(p, q\) 即可，条件变成了 \(p, q\) 互质，这样才能保证 \(d\) 是 \(\gcd\) 。
\[
\begin{aligned}
&=\sum\limits_{i = 1}^{\lfloor\frac n k\rfloor}\sum\limits_{j = 1}^{\lfloor\frac m k\rfloor}[\gcd(i, j) = 1]
\end{aligned}
\]
然后就没了，注意应用这个条件的时候不要只是傻傻的除过去，想清楚它的本质是什么。
Code
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 // author : black_trees
#include <cmath> #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 = 5e4 + 10 ; int n , m ; int mu [ si ], s [ si ]; int pri [ si ], vis [ si ]; void sieve ( int n ) { mu [ 1 ] = 1 , m = 0 ; memset ( vis , 0 , sizeof vis ); for ( int i = 2 ; i <= n ; ++ i ) { if ( ! vis [ i ]) pri [ ++ m ] = i , vis [ i ] = i , mu [ i ] = -1 ; for ( int j = 1 ; j <= m ; ++ 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 ; break ; } mu [ i * pri [ j ]] = - mu [ i ]; } } s [ 1 ] = mu [ 1 ]; for ( int i = 2 ; i <= n ; ++ i ) s [ i ] = s [ i - 1 ] + mu [ i ]; } int sum ( int l , int r ) { return s [ r ] - s [ l - 1 ]; }
int F ( int n , int m ) { if ( n > m ) swap ( n , m ); int ret = 0 ; int l = 1 , r = 0 ; for (; l <= n && n / l && m / l ; ) { r = min (( n / ( n / l )), ( m / ( m / l ))), r = min ( r , n ); ret += sum ( l , r ) * ( n / l ) * ( m / l ), l = r + 1 ; } return 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 -- ) { int a , b , k ; cin >> a >> b >> k ; cout << F ( a / k , b / k ) << endl ; } return 0 ; }
题外话，理解了这些我以前看都不敢看的式子之后打代码特别轻快，不知道为啥。
Ex.3
求 \(\sum\limits_{i = 1}^n\sum\limits_{j = 1}^m \gcd(i, j)\) 。
\(n \le m \le 10^7\) 。
思考一下这个 \(\gcd(i, j)\) 怎么拆。
如果是分类尝试枚举看起来太 shaber 了，我们考虑一个高级点的做法。
我们之前提到过恒等函数 \(Id_1(n) = n\) ，也就是说我们可以用他替换任意的变量。
而注意到 \(Id_1 = Id_0 * \varphi\) ，所以我们可以这么做：
原式：
\[
\begin{aligned}
&=\sum\limits_{i = 1}^n\sum\limits_{j = 1}^m Id_1(\gcd(i, j)) \\
&=\sum\limits_{i = 1}^n\sum\limits_{j = 1}^m \sum\limits_{d | \gcd(i, j)} Id_0(d)\varphi(d) \\
&= \sum\limits_{d = 1}^n \varphi(d)\sum\limits_{d | i}\sum\limits_{d | j} 1\\
&= \sum\limits_{d = 1}^n \varphi(d)\lfloor\dfrac n d\rfloor\lfloor\dfrac m d\rfloor
\end{aligned}
\]
又是个数论分块板子了，筛一下欧拉函数即可。
Ex.4
\(\sum\limits_{i = 1}^n\sum\limits_{j = 1}^m \operatorname{lcm}(i, j)\) ？
\(n \le m \le 10^7\) 。
Source: Luogu1829 - Crash的数字表格
好炫酷，上强度了。
\[
\begin{aligned}
&=\sum\limits_{i = 1}^n\sum\limits_{j = 1}^m \operatorname{lcm}(i, j) \\
&=\sum\limits_{i = 1}^n\sum\limits_{j = 1}^m \dfrac{ij}{\gcd(i, j)} \\
\end{aligned}
\]
直接拆 \(\gcd\) 太鲁莽了，这里是在分母，拆了也没用，先尝试把 \(\gcd(i, j)\) 拿出来。
考虑经典套路，枚举所有可能的 \(d\) ，替换 \(\gcd(i, j)\) ，然后用 Iverson Bracket 限制转化（就是 Ex.1 里 Tips 说的理解方式）。
\[
\begin{aligned}
&=\sum\limits_{d = 1}^n\sum\limits_{i = 1}^n\sum\limits_{j = 1}^m \dfrac{ij}{d}\cdot[\gcd(i, j) = d] \\
&=\sum\limits_{d = 1}^n\sum\limits_{i = 1}^{\lfloor\frac n d\rfloor}\sum\limits_{j = 1}^{\lfloor\frac m d\rfloor}\dfrac{ij}{d} \cdot d^2 \cdot [\gcd(i, j) = 1] \\
&=\sum\limits_{d = 1}^nd\sum\limits_{i = 1}^{\lfloor\frac n d\rfloor}\sum\limits_{j = 1}^{\lfloor\frac m d\rfloor}ij\cdot[\gcd(i, j) = 1]
\end{aligned}
\]
这里要乘上 \(d^2\) 是因为这里有 \(ij\) 乘积项的贡献，我们缩短 \(i, j\) 的枚举范围之后也会少算。
然后此时 \(d\) 其实就和 \(\gcd\) 分离了，它现在只是一个系数了。
我们考虑 \(\sum\limits_{i = 1}^{\lfloor\frac n d\rfloor}\sum\limits_{j = 1}^{\lfloor\frac m d\rfloor}ij\cdot[\gcd(i, j) = 1]\) 这一坨怎么算（为了少写一些，并且会清楚很多）。
反演！狠狠的反演！记这一块为 \(F(\lfloor\dfrac n d\rfloor,\lfloor\dfrac m d\rfloor)\) 。
\[
\begin{aligned}
F(n, m) &= \sum\limits_{i = 1}^n\sum\limits_{j = 1}^mij\cdot[\gcd(i,j) = 1]\\
&= \sum\limits_{i = 1}^n\sum\limits_{j = 1}^mij\cdot\sum\limits_{d | \gcd(i, j)}\mu(d) \\
&= \sum\limits_{d = 1}^n\mu(d)\cdot\sum\limits_{i = 1}^n\sum\limits_{j = 1}^mij,(d | i \land d | j) \\
&= \sum\limits_{d = 1}^n\mu(d)\cdot\sum\limits_{i = 1}^ni\sum\limits_{j = 1}^mj,(d | i \land d | j) \\
&= \sum\limits_{d = 1}^n\mu(d)\cdot\sum\limits_{i = 1}^{\lfloor\frac n d\rfloor} id\sum\limits_{j = 1}^{\lfloor\frac m d\rfloor}jd\\
&= \sum\limits_{d = 1}^n\mu(d)\cdot d^2 \sum\limits_{i = 1}^{\lfloor\frac n d\rfloor}i\sum\limits_{j = 1}^{\lfloor\frac m d\rfloor}j
\end{aligned}
\]
后面那一坨是简单求和，前面 \(\mu(d) d^2\) 可以筛出来，然后 \(F(n, m)\) 就可以整除分块了（因为 \(i,j\) 上界是整除，所以这一段在 \(d\) 属于某个块的时候全是相等的，就是整除分块的形式）。
带回原式：
\[
\sum\limits_{d = 1}^nd\cdot F(\lfloor\dfrac n d\rfloor,\lfloor\dfrac m d\rfloor)
\]
这部分也可以整除分块（下标里有整除，值自然也会分块）。
两次做完，复杂度 \(O(n + m)\) 。
Code
没绷住，sieve 的时候 cnt 用的是 m，全局，然后 main 里 sieve 完就寄了。
怪不得我改 sieve 的 limit 的时候结果会变。
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 // author : black_trees
#include <cmath> #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 = 1e7 + 10 ; const int mod = 20101009 ; int Qpow ( int a , int b ) { int ret = 1 % mod ; for (; b ; b >>= 1 ) { if ( b & 1 ) ret = ret * a % mod ; a = a * a % mod ; } return ret ; } int inv2 ; int n , m , t ; int mu [ si ], s [ si ]; int pri [ si ], vis [ si ]; void sieve ( int n ) { mu [ 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 ; for ( int j = 1 ; j <= t ; ++ j ) { if ( pri [ j ] > vis [ i ] || pri [ j ] * i > n ) break ; vis [ i * pri [ j ]] = pri [ j ]; if ( i % pri [ j ] == 0 ) { mu [ i * pri [ j ]] = 0 ; break ; } mu [ i * pri [ j ]] = - mu [ i ]; } } s [ 1 ] = mu [ 1 ], s [ 1 ] = ( s [ 1 ] + mod ) % mod ; for ( int i = 2 ; i <= n ; ++ i ) s [ i ] = ( s [ i - 1 ] % mod + 1l l * i % mod * i % mod * ( mu [ i ] + mod ) % mod ) % mod ;
} int sum ( int n ) { return ( n % mod ) * ( n + 1 ) % mod * inv2 % mod ; } int sum2 ( int n , int m ) { return sum ( n ) % mod * sum ( m ) % mod ; } int F ( int n , int m ) { int ret = 0 ; int l = 1 , r = 0 ; for (; l <= n && n / l && m / l ; ) { r = min (( n / ( n / l )), ( m / ( m / l ))), r = min ( r , n ); ret = ( ret + ( s [ r ] - s [ l - 1 ] + mod ) % mod * sum2 ( n / l , m / l ) % mod ) % mod ; l = r + 1 ; } return ret ; } int G ( int n , int m ) { int ret = 0 ; int l = 1 , r = 0 ; for (; l <= n && n / l && m / l ; ) { r = min (( n / ( n / l )), ( m / ( m / l ))), r = min ( r , n ); ret = ( ret + ( r - l + 1 ) % mod * ( r + l ) % mod * inv2 % mod * F ( n / l , m / l ) % mod ) % mod ; l = r + 1 ; } return ret ; } signed main () { cin . tie ( 0 ) -> sync_with_stdio ( false ); cin . exceptions ( cin . failbit | cin . badbit ); cin >> n >> m ; if ( n > m ) swap ( n , m ); sieve ( n ), inv2 = Qpow ( 2 , mod - 2 ); cout << G ( n , m ) % mod << endl ; return 0 ; }
Ex.5
求 \(\left(\sum\limits_{i = 1}^n\sum\limits_{j = 1}^n ij\gcd(i, j)\right)\bmod p\)
\(\text{PRIME(p)}, p \in [5 \times 10^8, 1.1\times 10^9], n \le 10^{10}\) 。
Source: Luogu3768 - 简单的数学题
\[
\begin{aligned}
&=\sum\limits_{i = 1}^n\sum\limits_{j = 1}^n i \cdot j \cdot \gcd(i, j) \\
&= \sum\limits_{i = 1}^n\sum\limits_{j = 1}^n i \cdot j \cdot Id_1(\gcd(i, j)) \\
&= \sum\limits_{i = 1}^n\sum\limits_{j = 1}^n i \cdot j \cdot \sum\limits_{d | \gcd(i, j)} \varphi(d) \\
&= \sum\limits_{d = 1}^n\varphi(d)\sum\limits_{d | i}^{i \in [1, n]}\sum\limits_{d | j}^{j \in [1, n]}i\cdot j \\
&= \sum\limits_{d = 1}^n\varphi(d)\sum\limits_{i = 1}^{\lfloor\frac n d\rfloor}\sum\limits_{j = 1}^{\lfloor\frac n d\rfloor} i\cdot j \cdot d^2 \\
&= \sum\limits_{d = 1}^n\varphi(d)\cdot d^2 \sum\limits_{i = 1}^{\lfloor\frac n d\rfloor}\sum\limits_{j = 1}^{\lfloor\frac n d\rfloor}i\cdot j
\end{aligned}
\]
记 \(F(n) = \sum\limits_{i = 1}^n\sum\limits_{j = 1}^ni\cdot j \iff F(n) = \left(\dfrac{n (n + 1)}{2}\right)^2\)
原式：
\[
\begin{aligned}
&= \sum\limits_{d = 1}^n \varphi(d) \cdot d^2 \cdot F(\lfloor\dfrac{n}{d}\rfloor)
\end{aligned}
\]
一眼丁真，鉴定为，纯纯的整除分块。
问题在于怎么在 \(n = 10^{10}\) 的时候求 \(\varphi(d)\cdot d^2\) 的前缀和。
这部分需要杜教筛 ，然后就做完了。
具体来说因为要求前缀和的 \(f = \varphi(d) \cdot i^2\) ，我们构造 \(g = Id_2\) 卷上去就行了。
这种后面跟一个变量的几次方，一般都是考虑 Dirichlet 卷积的形式构造 \(Id_k\) 卷上去。
Code
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 // 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 = 1e7 + 10 ; int n , m , t , mod , inv6 , inv4 ; int Qpow ( int a , int b ) { int ret = 1 % mod ; for (; b ; b >>= 1 ) { if ( b & 1 ) ret = ret * a % mod ; a = a * a % mod ; } return ret % mod ; } int phi [ si ], sphi [ si ]; int pri [ si ], vis [ si ]; void sieve ( int n ) { t = 0 , phi [ 1 ] = 1 ; memset ( vis , 0 , sizeof vis ); for ( int i = 2 ; i <= n ; ++ i ) { if ( ! vis [ i ]) pri [ ++ t ] = i , vis [ i ] = i , 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 ) { phi [ i * pri [ j ]] = phi [ i ] * pri [ j ]; break ; } phi [ i * pri [ j ]] = phi [ i ] * ( pri [ j ] - 1 ); } phi [ i ] = 1l l * i % mod * i % mod * phi [ i ] % mod ; ( phi [ i ] += phi [ i - 1 ] % mod ) %= mod ; } } std :: bitset < si > vi ; int id ( int x ) { return n / x ; } int sum2 ( int n ) { return n % mod * ( n + 1 ) % mod * ( 2 * n % mod + 1 ) % mod * inv6 % mod ; } int sum3 ( int n ) { return ( n % mod * ( n + 1 ) % mod ) % mod * ( n % mod * ( n + 1 ) % mod ) % mod * inv4 % mod ; } int Phi ( int n ) { if ( n <= m ) return phi [ n ] % mod ; if ( vi [ id ( n )]) return sphi [ id ( n )] % mod ; int ret = sum3 ( n ) % mod ; int l = 2 , r = 0 ; for (; l <= n ; ) { r = min ( n , ( n / ( n / l ))); ret = ( ret - (( sum2 ( r ) - sum2 ( l - 1 ) + mod ) % mod * Phi ( n / l ) % mod ) + mod ) % mod ; l = r + 1 ; } vi [ id ( n )] = true ; return sphi [ id ( n )] = ret % mod ; } int solve ( int n ) { int ret = 0 ; int l = 1 , r = 0 ; for (; l <= n ; ) { r = min ( n , ( n / ( n / l ))); ( ret += ( Phi ( r ) - Phi ( l - 1 ) + mod ) % mod * sum3 ( n / l ) % mod ) %= mod ;
l = r + 1 ; } return ret ; } signed main () { cin . tie ( 0 ) -> sync_with_stdio ( false ); cin . exceptions ( cin . failbit | cin . badbit ); cin >> mod >> n ; sieve ( m = si - 10 ); // 1e10^{2/3} ~ 1e7
inv4 = Qpow ( 4 , mod - 2 ); inv6 = Qpow ( 6 , mod - 2 ); cout << solve ( n ) << endl ; return 0 ; }
习题
Luogu2522 - Problem b
Zap 那题的升级版，现在是 \(i \in[a,b], j \in[c,d]\) 。
很简单啊，考虑一个广义的二维前缀和即可。
Code
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 // author : black_trees
#include <cmath> #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 = 5e4 + 10 ; int n , m ; int mu [ si ], s [ si ]; int pri [ si ], vis [ si ]; void sieve ( int n ) { mu [ 1 ] = 1 , m = 0 ; memset ( vis , 0 , sizeof vis ); for ( int i = 2 ; i <= n ; ++ i ) { if ( ! vis [ i ]) pri [ ++ m ] = i , vis [ i ] = i , mu [ i ] = -1 ; for ( int j = 1 ; j <= m ; ++ 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 ; break ; } mu [ i * pri [ j ]] = - mu [ i ]; } } s [ 1 ] = mu [ 1 ]; for ( int i = 2 ; i <= n ; ++ i ) s [ i ] = s [ i - 1 ] + mu [ i ]; } int sum ( int l , int r ) { return s [ r ] - s [ l - 1 ]; }
int F ( int n , int m , int k ) { if ( n > m ) swap ( n , m ); int ret = 0 ; int l = 1 , r = 0 ; for (; l <= n / k ; ) { r = min (( n / ( n / l )), ( m / ( m / l ))), r = min ( r , n / k ); ret += sum ( l , r ) * ( n / ( l * k )) * ( m / ( l * k )); l = r + 1 ; } return 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 -- ) { int a , b , c , d , k ; cin >> a >> b >> c >> d >> k ; cout << F ( b , d , k ) - F ( a - 1 , d , k ) - F ( b , c - 1 , k ) + F ( a - 1 , c - 1 , k ) << endl ; } return 0 ; }
HDU1695 - GCD
还是 Zap 的升级版，但是现在 \((1, 2),(2, 1)\) 算成一对。
并且 \(n, m\) 不一样。
注意到 \(n, m\) 不一样，所以我们不能简单的除以二去重，我们需要考虑别的方式。
对于这种重复计数，最简单的处理方法就是钦定大小，因为 \(= k\) 和 \(= 1\) 就是自变量除一下所以我直接考虑 \(= 1\) 了。
不妨设 \(m < n\) ，答案为：
\[
\begin{aligned}
F(n,m)&=\sum\limits_{i = 1}^n\sum\limits_{j = 1}^m[\gcd(i, j) = 1]\cdot [j \le i] \\
&=\sum\limits_{d = 1}^n\mu(d) \sum\limits_{i = 1}^{\lfloor\frac n d\rfloor}\sum\limits_{j = 1}^{\min(i, \lfloor\frac m d\rfloor)}1 \\
&=\sum\limits_{d = 1}^n\mu(d)(\sum\limits_{i = 1}^{\lfloor\frac m d\rfloor}\sum\limits_{j = 1}^{i}1+\sum\limits_{i = \lfloor\frac m d\rfloor + 1}^{\lfloor\frac n d\rfloor}\sum\limits_{j = 1}^{\lfloor\frac n d\rfloor}1)\\
&=\sum\limits_{d = 1}^n\mu(d)(\sum\limits_{i = 1}^{\lfloor\frac m d\rfloor}i + \lfloor\dfrac m d\rfloor(\lfloor\dfrac n d\rfloor - \lfloor\dfrac m d\rfloor))
\end{aligned}
\]
答案是 \(F(\lfloor\dfrac{n}{k}\rfloor,\lfloor\dfrac{m}{k}\rfloor)\) 。
Code
记得推式子的时候不要莫名其妙加 \(\sum\) 和变量！一步一步来不要跳！
如果整除分块的函数有两个元记得 \(r\) 的取值，cnmd.
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 // author : black_trees
#include <cmath> #include <cctype> #include <cstdio> #include <cstring> #include <iostream> #include <algorithm> #define endl '\n'
#define int __int128_t
using namespace std ; using i64 = long long ; template < typename __Tp > void read ( __Tp & x ) { int f = x = 0 ; char ch = getchar (); for (; ! isdigit ( ch ); ch = getchar ()) if ( ch == '-' ) f = 1 ; for (; isdigit ( ch ); ch = getchar ()) x = ( x << 1 ) + ( x << 3 ) + ( ch ^ 48 ); if ( f ) x = - x ; } void read ( char & ch ) { for ( ch = getchar (); isspace ( ch ); ch = getchar ()); } template < typename __Tp1 , typename ... __Tp2 > void read ( __Tp1 & x , __Tp2 & ... y ) { read ( x ), read ( y ...); } template < typename __Tp > void write ( __Tp x ) { if ( x < 0 ) putchar ( '-' ), x = - x ; if ( x > 9 ) write ( x / 10 ); putchar ( x % 10 + 48 ); } void write ( char ch ) { putchar ( ch ); } void write ( const char * s ) { for (; * s ; ++ s ) putchar ( * s ); } template < typename __Tp1 , typename ... __Tp2 > void write ( __Tp1 x , __Tp2 ... y ) { write ( x ), write ( y ...); } const int si = 1e5 + 10 ; int n , m , t ; int mu [ si ], s [ si ]; int pri [ si ], vis [ si ]; void sieve ( int n ) { mu [ 1 ] = 1 , t = 0 , s [ 1 ] = 1 ; memset ( vis , 0 , sizeof vis ); for ( int i = 2 ; i <= n ; ++ i ) { if ( ! vis [ i ]) pri [ ++ t ] = i , vis [ i ] = i , mu [ i ] = -1 ; for ( int j = 1 ; j <= t ; ++ j ) { if ( pri [ j ] > vis [ i ] || pri [ j ] * i > n ) break ; vis [ i * pri [ j ]] = pri [ j ]; if ( i % pri [ j ] == 0 ) { mu [ i * pri [ j ]] = 0 ; break ; } mu [ i * pri [ j ]] = - mu [ i ]; } s [ i ] = s [ i - 1 ] + mu [ i ]; } } int get ( int n , int m , int d ) { int u = n / d , v = m / d ; return ( v * ( v + 1 )) / 2 + ( v * ( u - v )); } int F ( int n , int m ) { if ( n < m ) swap ( n , m ); int ret = 0 ; int l = 1 , r = 0 ; for (; l <= m && m / l && n / l ; ) { r = min ({ m , m / ( m / l ), n / ( n / l )}); ret += ( s [ r ] - s [ l - 1 ]) * get ( n , m , l ), l = r + 1 ; } return ret ; } signed main () { sieve ( si - 10 ); int T ; read ( T ); for ( int nw = 1 ; nw <= T ; ++ nw ) { int a , b , c , d , k ; read ( a , b , c , d , k ); write ( "Case " , nw , ": " ); if ( ! k ) write ( "0" , endl ); else write ( F ( b / k , d / k ), endl ); } return 0 ; }
比较趣味，这题都多少年了，网上居然没有我这个显而易见的做法（）
可见大部分选手的独立思考能力有待（欸好像我也是）
Luogu2303 - Longge 的问题
求 \(\sum\limits_{i = 1}^n\gcd(i, n)\) ，\(1\le n < 2^{32}\) 。
很经典的问题啊！
显然就是 \(\sum\limits_{d | n}\varphi(d)\cdot\dfrac{n}{d}\)
\(\varphi\) 就不要线性筛了，直接 \(O(\sqrt n)\) 求就行。
这题似乎可以再次化简，但是我不会。
Code
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 // author : black_trees
#include <cmath> #include <cstdio> #include <cstring> #include <iostream> #include <algorithm> #define endl '\n'
#define int long long
using namespace std ; using i64 = long long ; int n , ans = 0 ; int phi ( int x ) {
int ans = 1 ; bool flag = false ; for ( int i = 2 ; i * i <= x ; i ++ )
if ( ! ( x % i )) { flag = true ; break ; } if ( ! flag && x != 1 ) return x - 1 ;
for ( int i = 2 ; i * i <= x ; i ++ ) { if ( ! ( x % i )) { ans *= ( i - 1 ), x /= i ;
while ( ! ( x % i )) x /= i , ans *= i ; } } if ( x != 1 ) ans *= x - 1 ;
return ans ;
} signed main () { cin . tie ( 0 ) -> sync_with_stdio ( false ); cin . exceptions ( cin . failbit | cin . badbit ); cin >> n ; for ( int i = 1 ; i * i <= n ; i ++ ) { if ( ! ( n % i )) { ans += phi ( i ) * ( n / i );
if ( i * i < n ) ans += phi ( n / i ) * i ;
} } cout << ans << endl ; return 0 ; }
Luogu3327 - 约数个数和
求 \(\sum\limits_{i = 1}^n\sum\limits_{j = 1}^m\sigma_0(ij)\) 。
\(n, m, T \le 50000\) 。
有一个比较蜜蜂的想法，据说叫 rng_58 - clj 等式：
\[
\sigma_0(ij) = \sum\limits_{x | i}\sum\limits_{y | j}[\gcd(x, y) = 1]
\]
然后我们先不着急带进去。
如果直接带进去反演，有 \(i, j\) 在前面管着，很烦躁，所以我们先交换一下 Sigma：
\[
\begin{aligned}
ans &= \sum\limits_{i = 1}^n\sum\limits_{j = 1}^m\sum\limits_{x = 1}^i\sum\limits_{y = 1}^j[\gcd(x, y) = 1] \\
&= \sum\limits_{x = 1}^n\sum\limits_{y = 1}^m\sum\limits_{x | i}^{i \in[1, n]}\sum\limits_{y | j}^{j\in[1,m]}[\gcd(x, y) = 1] \\
&= \sum\limits_{x = 1}^n\sum\limits_{y = 1}^m\lfloor\dfrac n x\rfloor\lfloor\dfrac m y\rfloor[\gcd(x, y) = 1] \\
&= \sum\limits_{d = 1}^n\mu(d)\sum\limits_{d | x}\sum\limits_{d | y}\lfloor\dfrac n x\rfloor\lfloor\dfrac m y\rfloor\\
&= \sum\limits_{d= 1}^n\mu(d)\cdot g(\lfloor\dfrac n d\rfloor)\cdot g(\lfloor\dfrac m d\rfloor)
\end{aligned}
\]
\(g\) 可以提前预处理，整除分块即可。
Code
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 // author : black_trees
#include <cmath> #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 = 5e4 + 10 ; int n , m , t ; int mu [ si ], s [ si ]; int pri [ si ], vis [ si ]; void sieve ( int n ) { t = 0 , mu [ 1 ] = 1 ; memset ( vis , 0 , sizeof vis ); for ( int i = 2 ; i <= n ; ++ i ) { if ( ! vis [ i ]) pri [ ++ t ] = i , vis [ i ] = i , mu [ 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 ; break ; } mu [ i * pri [ j ]] = - mu [ i ]; } mu [ i ] += mu [ i - 1 ]; } for ( int i = 1 ; i <= n ; ++ i ) for ( int j = i ; j <= n ; j += i ) s [ j ] ++ ; for ( int i = 1 ; i <= n ; ++ i ) s [ i ] += s [ i - 1 ]; } int F ( int n , int m ) { if ( n > m ) swap ( n , m ); int ret = 0 ; int l = 1 , r = 0 ; for (; l <= n && n / l && m / l ; ) { r = min ({ n , n / ( n / l ), m / ( m / l )}); ret += ( mu [ r ] - mu [ l - 1 ]) * s [ n / l ] * s [ m / l ]; l = r + 1 ; } return ret ; } signed main () { cin . tie ( 0 ) -> sync_with_stdio ( false ); cin . exceptions ( cin . failbit | cin . badbit ); sieve ( si - 10 );
int T ; cin >> T ; while ( T -- ) { cin >> n >> m ; cout << F ( n , m ) << endl ; } return 0 ; }
Luogu2257 - YY的GCD
求 \(\sum\limits_{i = 1}^{n}\sum\limits_{j = 1}^{m}[\gcd(i, j) \in \text{prime}]\) 。
\(T = 10^4, n, m \le 10^7\) 。
注意到对 \(\gcd\) 有限制，为了快速换出反演式子，我们钦定 \(\gcd\) ，并将其提出（枚举 \(\gcd\) ）。
\[
\begin{aligned}
&=\sum\limits_{d \in \text{prime}}\sum\limits_{i = 1}^{n}\sum\limits_{j = 1}^{m}[\gcd(i, j) = d] \\
&=\sum\limits_{d \in \text{prime}}^{}\sum\limits_{i = 1}^{\lfloor\frac{n}{d}\rfloor}\sum\limits_{j = 1}^{\lfloor\frac{m}{d}\rfloor}[\gcd(i, j) = 1]\\
&=\sum\limits_{d \in \text{prime}}^{}\sum\limits_{i = 1}^{\lfloor\frac{n}{d}\rfloor}\sum\limits_{j = 1}^{\lfloor\frac{m}{d}\rfloor}\sum\limits_{t | \gcd(i, j)}^{}\mu(t)\\
&=\sum\limits_{d \in \text{prime}}^{}\sum\limits_{t = 1}^{\lfloor\frac{n}{d}\rfloor}\mu(t)\sum\limits_{i = 1}^{\lfloor\frac{n}{dt}\rfloor}\sum\limits_{j = 1}^{\lfloor\frac{m}{dt}\rfloor}1\\
&=\sum\limits_{d \in \text{prime}}^{}\sum\limits_{t = 1}^{\lfloor\frac{n}{d}\rfloor}\mu(t)\lfloor\dfrac{n}{dt}\rfloor\lfloor\dfrac{m}{dt}\rfloor\\
\end{aligned}
\]
反演那一步不要记录成 \(t | d\) ，这样会增加变量之间的耦合度，时刻记住我们要做的是降低耦合度拆分变量，而我们知道 \(i, j\) 一定会被化掉，所以我们记录成 \(t | \gcd(i, j)\)
这个看起来没法很快做，但是可以记录 \(T = dt\) 。
然后式子变成：
\[
\begin{aligned}
&=\sum\limits_{T = 1}^n \lfloor\dfrac{n}{T}\rfloor\lfloor\dfrac{m}{T}\rfloor\cdot \sum\limits_{d \in \text{prime}, d | T}\sum\limits_{t = 1}^{\lfloor\frac{n}{d}\rfloor}\mu(t)
\end{aligned}
\]
后面这个就可以预处理了，对每个质数，给他的倍数的贡献加一下就可以了。
记得开 O2，不然 TLE/cy
Code
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 // author : black_trees
#include <cmath> #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 = 1e7 + 10 ; int n , m , t ; int mu [ si ], f [ si ]; int pri [ si ], vis [ si ]; void sieve ( int n ) { t = 0 , mu [ 1 ] = 1 ; memset ( vis , 0 , sizeof vis ); for ( int i = 2 ; i <= n ; ++ i ) { if ( ! vis [ i ]) pri [ ++ t ] = i , vis [ i ] = i , mu [ 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 ; break ; } mu [ i * pri [ j ]] = - mu [ i ]; } } for ( int i = 1 ; i <= t ; ++ i ) { for ( int j = 1 ; j * pri [ i ] <= n ; ++ j ) { f [ pri [ i ] * j ] += mu [ j ]; } } f [ 0 ] = 0 , mu [ 0 ] = 0 ; for ( int i = 1 ; i <= n ; ++ i )
f [ i ] += f [ i - 1 ], mu [ i ] += mu [ i - 1 ]; } int F ( int n , int m ) { if ( n > m ) swap ( n , m ); int ret = 0 ; int l = 1 , r = 0 ; for (; l <= n ; ) { r = min ({ n , n / ( n / l ), m / ( m / l )}); ret += ( f [ r ] - f [ l - 1 ]) * ( n / l ) * ( m / l ); l = r + 1 ; } return ret ; } signed main () { cin . tie ( 0 ) -> sync_with_stdio ( false ); cin . exceptions ( cin . failbit | cin . badbit ); sieve ( si - 10 ); int T ; cin >> T ; while ( T -- ) { cin >> n >> m ; cout << F ( n , m ) << endl ; } return 0 ; }
SPOJ - LCMSUM
求 \(\sum\limits_{i = 1}^n \operatorname{lcm}(i, n)\) 。
\(T \le 3 \times 10^5, 1 \le n \le 10^6\) 。
老办法。
\[
\begin{aligned}
answer &= \sum\limits_{i = 1}^{n}\dfrac{in}{\gcd(i, n)} \\
&= n\sum\limits_{i = 1}^{n}\dfrac{i}{\gcd(i, n)} \\
&= n\sum\limits_{d = 1}^n\sum\limits_{i = 1}^n\dfrac{i}{d}\cdot[\gcd(i, n) = d] \\
&= n\sum\limits_{d = 1}^{n}\sum\limits_{i = 1}^{\lfloor\frac{n}{d}\rfloor}i\cdot[\gcd(i, n / d) = 1]
\end{aligned}
\]
后面这部分提出来反演化简即可，前面的部分就是整除分块。
注意除掉 \(d\) 那个地方不要忘记这个操作本身的含义，所以 \(n \to n / d\) 。
懒得推了。
Luogu3704 - 数字表格
求 \(\prod\limits_{i = 1}^{n}\prod\limits_{j = 1}^{m}F(\gcd(i, j))\) 。
其中 \(F\) 为斐波那契数列，答案对 \(p\) 取模。
\(1\le T \le 10^3, 1\le n, m \le 10^6\) 。
Simple Problem.
首先 \(\gcd(i, j)\) 在下标里面，所以枚举 \(\gcd(i, j)\) 将其提出来的做法是显然的。
但是一个问题在于，现在这里求的是 \(\prod\) 不是 \(\sum\) ，如果我们用 Iverson Bracket 转化，那最后一定会出现 \(0\) ，整个式子就变成 \(0\) 没法做了。
所以我们思考一下怎么很好的表达这个限制。
其实本质是，我们希望不合法的元素在 \((\oplus, \mathbb Z_p)\) 意义下为单位元，如果求 \(\sum\) ，那么 \(\oplus = +\) ，\(e = 0\) 。
所以现在是 \((\times, \mathbb Z_p) \iff e = 1\) 。
那么我们不难想到 \(x^0 = 1 = e\) ，所以我们把 Iverson Bracket 放在指数。
\[
\begin{aligned}
answer &= \prod\limits_{d = 1}^{n}\prod\limits_{i = 1}^{n}\prod\limits_{j = 1}^{m} F(d)^{[\gcd(i, j) = d]}
\end{aligned}
\]
注意到其实后面两个 \(\prod\) 就是在数 \((1, 1) \sim (n, m)\) 中有多少 \((i, j)\) 满足 \(\gcd(i, j) = d\) ，我们记这个为 \(G(n, m, d)\) 。
原式：
\[
\begin{aligned}
&=\prod\limits_{d = 1}^{n}F(d)^{G(n, m, d)}
\end{aligned}
\]
\(G\) 就是我们的老朋友（Ex.2 中的答案），可以反演。
然后式子变为：
\[
\begin{aligned}
&=\prod\limits_{d = 1}^{n}F(d)^{\sum\limits_{t = 1}^{\lfloor\frac{n}{d}\rfloor}\mu(t)\lfloor\frac{n}{dt}\rfloor\lfloor\frac{m}{dt}\rfloor}
\end{aligned}
\]
看起来已经可以做了对吧？但是 \(T = 10^3\) ，我们需要单次询问更快的做法，这里单次询问是 \(O(n)\) 的。
然后定睛一看，这个 \(dt\) 的形式，哎呀，这不是我们 YY 的 GCD 吗，枚举 \(T = dt\) ！
\[
\begin{aligned}
&=\prod\limits_{T = 1}^{n}f(T)^{\lfloor\frac{n}{T}\rfloor\lfloor\frac{m}{T}\rfloor}
\end{aligned}
\]
其中 \(f(T) = \prod\limits_{d | T}^{}F(d)^{\mu(\frac{T}{d})}\) 。
\(f\) 通过倍数法可以算，然后整除分块即可。
代码暂略，之后写。
Luogu5518 - 幽灵乐团 / 莫比乌斯反演基础练习题
东风谷 早苗(Kochiya Sanae)非常喜欢幽灵乐团的演奏，她想对她们的演奏评分。
因为幽灵乐团有 \(3\) 个人，所以我们可以用 \(3\) 个正整数 \(A,B,C\) 来表示出乐团演奏的分数，她们的演奏分数可以表示为
\[\prod_{i=1}^{A}\prod_{j=1}^{B}\prod_{k=1}^{C}\left(\frac{\text{lcm}(i,j)}{\gcd(i,k)}\right)^{f(type)}\]
因为音乐在不同的部分会有不同的听觉感受，所以 \(type\) 会在 \(\{0,1,2\}\) 中发生变化，其中：
\[ \begin{aligned} f(0)&=1 \cr f(1)&=i \times j \times k \cr f(2)&=\gcd(i,j,k) \end{aligned} \]
\[ 1\leq A,B,C\leq 10^5 \ \ \ \ 10^7 \leq p \leq 1.05\times 10^9\ \ \ \ p\in \{ prime\} \ \ \ \ T =70\]
咕咕咕。
题外话
Source: BZOJ2498/HDU4093 Xavier is Learning to Count
咕咕咕，这个有点难，之后写。
但是它引入了一个结论：
一个满足 \(f(1) \not = 0\) 的数论函数，用它做 Dirichlet 卷积是可逆的。
即存在 \(g\) 使得 \(f * g = \epsilon\) 。
这个其实是显然的，因为 \(\epsilon(1) = 1\) ，如果 \(f(1) = 0\) 了怎么卷都卷不回去。
以上大部分内容参考 zyw 学姐的讲义，例题参考 lsj 学长的博客，习题来自 CYJian 的题单，在此感谢。
最后更新:
July 10, 2023