莫比乌斯函数

莫比乌斯函数

主要是 莫比乌斯反演，前面已经 用过几次了，但是没有 专门讲

建议阅读 数论函数基础，数论分块，欧拉函数 章节，了解基本概念与先要知识

此处获取本节调试数据 / 代码包

全文 绝大多数 内容是对 [0] 中讲述的 粗略抄写 和 胡乱加工

1. 概念

关于 莫比乌斯函数，前面 数论函数基础 章节中有所提到，这里只给出 定义式

\[\begin {aligned} \mu (n) = \begin {cases} 1 & n = 1 \\ 0 & \exists d > 1, d ^ 2 \mid n \\ (-1) ^ {\omega (n)} & \operatorname {otherwise} \end {cases} \end {aligned} \]

判定是否等于 \(1\)，是否有 平方因子

但其实它最初被发现的时候长这样（在 数论函数基础 中作为一个 结论 提出了）

\[\sum _ {d \mid n} \mu (d) = [n = 1] \]

一种转化过程在 数论函数基础 一节中已给出，考虑展开后 二项式定理 即可，这里不做赘述

下面会描述另外一种转化方式

上述式子经常以 \(\sum _ {d \mid \gcd (i, j)} \mu (d) = [i \perp j]\) 的形式出现，容易发现其本质相同

在 \(1857\) 年，\(\textsf {Richard Dedekind}\) 和 \(\textsf {Joseph Liouville}\) 注意到了如下重要的 反演原理

\[\begin {aligned} g (m) = \sum _ {d \mid m} f (d) \Longleftrightarrow f (m) = \sum _ {d \mid m} \mu (d) g \left ( \dfrac m d \right ) \end {aligned} \]

这就是我们常说的 莫比乌斯反演，考虑证明

先来考虑 \(g (m) = \sum _ {d \mid m} f (d) \Longrightarrow f (m) = \sum _ {d \mid m} \mu (d) g \left ( \dfrac m d \right )\)

\[\begin {aligned} \sum _ {d \mid m} { \mu (d) g \left ( \dfrac m d \right ) } &= \sum _ {d \mid m} { \mu \left ( \dfrac m d \right ) g (d) } \\ &= \sum _ {d \mid m} { \mu \left ( \dfrac m d \right ) \sum _ {k \mid d} f (k) } \\ &= \sum _ {k \mid d} { \sum _ { d \mid \frac m k } { \mu \left ( \dfrac m {kd} \right ) f (k) } } \\ &= \sum _ {k \mid d} { \sum _ { d \mid \frac m k } { \mu (d) f (k) } } \\ &= \sum _ {k \mid d} { \left [ \dfrac m k = 1 \right ] f (k) } \\ &= f (k) \end {aligned} \]

其中 \(Line ~ 3 \to Line ~ 4\) 比较有趣，注意到 \(d \mid \dfrac m k\)，故 \(d\) 与 \(\dfrac m {kd}\) 一定 对称且一一对应

然后考虑 \(g (m) = \sum _ {d \mid m} f (d) \Longleftarrow f (m) = \sum _ {d \mid m} \mu (d) g \left ( \dfrac m d \right )\)

\[\begin {aligned} \sum _ {d \mid m} f (d) &= \sum _ {d \mid m} { \sum _ {k \mid d} { \mu (k) g \left ( \dfrac d k \right ) } } \\ &= \sum _ {d \mid m} { \sum _ {k \mid d} { \mu \left ( \dfrac d k \right ) g (k) } } \\ &= \sum _ {k \mid m} { \sum _ { d \mid \frac m k } { \mu (d) g (k) } } \\ &= \sum _ {k \mid m} { \left [ \dfrac m k = 1 \right ] g (k) } \\ &= g (m) \end {aligned} \]

对于 \(Line ~ 2 \to Line ~ 3\)，容易发现，下面的 \(d\) 等效于上面的 \(dk\)，于是 \(\dfrac d k \to \dfrac {dk} {k} \to d\)

于是得证

我们可以用 狄利克雷卷积 很好的表示 莫比乌斯反演 的式子

\[\begin {aligned} g = f * 1 \Longleftrightarrow f = g * \mu \end {aligned} \]

于是回到 发现时的式子，我们发现其可以表示为

\[\mu * 1 = \epsilon \]

注意到 积性函数的逆元也是积性函数，故由此我们也可以说明 莫比乌斯函数是积性函数

根据 积性函数的性质，我们知道我们只需得出 \(\mu (p ^ k)\) 就能算出 \(\mu (n)\)

当 \(n = p ^ k, k \ge 1\) 时，根据 发现时的式子，我们有 \(\sum _ {i = 0} ^ k \mu (p ^ i) = [n = 1] = 0\)

容易得到 \(\mu (p) = -1, \mu (p ^ k) =0 ~ (k > 1)\) 的结论，于是 定义式 就出来了

\[\begin {aligned} \mu (n) = \begin {cases} 1 & n = 1 \\ 0 & \exists d > 1, d ^ 2 \mid n \\ (-1) ^ {\omega (n)} & \operatorname {otherwise} \end {cases} \end {aligned} \]

我们继续对 欧拉函数 一节中的 欧拉反演 加以探究（证明时笔者就采取了 莫比乌斯反演 的思路）

注意到欧拉反演可以写成 \(\sum _ {d \mid n} \varphi (d) = \operatorname {id} (n)\) 的形式

其中 \(\operatorname {id}\) 是 恒等函数 / 幂函数，在 数论函数基础 一节中有所介绍

于是使用 莫比乌斯反演 即可得到 \(\varphi (n) = \sum_{d \mid n} \mu (d) \operatorname {id} \left ( \dfrac n d \right)\) 即 \(\mu * \operatorname {id} = \varphi\)

我们也可以从 容斥系数 的角度来理解 莫比乌斯函数，具体请参考 [0] 中内容，此处不表

2. 线性筛莫比乌斯函数

我们知道了 莫比乌斯函数 是 积性函数，再根据 定义式，线性筛是简单的

inline void Sieve () {
	mu[1] = 1;
	for (int i = 2; i <= N; ++ i) {
		if (!vis[i]) pri[++ Cnt] = i, mu[i] = - 1;
		for (int k = 1; k <= Cnt && pri[k] * i <= N; ++ k) {
			vis[pri[k] * i] = 1;
			if (i % pri[k] == 0) break ;
			mu[i * pri[k]] = - mu[i];
		}
	}
}

3. Tricks

在应用 莫比乌斯反演 的时候，我们常常可以见到类似下面式子的形式

\[\begin {aligned} \sum _ {d \mid \gcd (i, j)} \mu (d) = [\gcd (i, j) = 1] \end {aligned} \]

这其实就是 上面提到过的 发现时的式子的等价形式

而有一个经典的套路，在需要枚举 \(i, j\) 的时候，我们可以枚举 \(\gcd\) 然后计算合法 \((i, j)\) 对数

显然，当且仅当 \(d \mid i\wedge d \mid j\) 时 \((i, j)\) 合法，这样我们就可以把 \(i, j\) 部分独立开了

注意，是 \(d\) 整除 \(i\) 与 \(d\) 整除 \(j\)，不是 \(d\) 或 \(i\) 与 \(d\) 或 \(j\)

Trick 1

我们可以这样解决下面的部分（这是很常见的）

\[\begin {aligned} \sum _ {i = 1} ^ n { \sum _ {j = 1} ^ m { [\gcd (i, j) = 1] } } &= \sum _ {i = 1} ^ n { \sum _ {j = 1} ^ m { \sum _ {d \mid \gcd (i, j)} \mu (d) } } \\ &= \sum _ {d = 1} ^ {\min (n, m)} { \mu (d) \sum _ {i = 1} ^ n { \sum _ {j = 1} ^ m { [d \mid i \wedge d \mid j] } } } \\ &= \sum _ {d = 1} ^ {\min (n, m)} { \mu (d) \left \lfloor \dfrac n d \right \rfloor \left \lfloor \dfrac m d \right \rfloor } \end {aligned} \]

这个式子就可以简单的用 整除分块 在 \(O (\sqrt n + \sqrt m)\) 的复杂度内解决，十分高效

对于上面的变换，我们有一种 感性的理解方式

相当于我们对 \(\gcd (i, j) = d\) 的 因数 进行了 容斥

即加上了 \(\gcd = 1\) 的数的贡献，减去 \(\gcd = p _ i ~ (p \in \mathbb P)\) 的贡献

加上 \(\gcd = p_i p_j ~ (p_i, p_j \in \mathbb P)\) 的贡献... 容易发现，这样的 容斥系数 就是 \(\mu\)

Trick 2

还有一个常用的 等式 如下

\[\begin {aligned} d (ij) = \sum _ {x \mid i} { \sum _ {y \mid j} { [x \perp y] } } \end {aligned} \]

一个证明思路是 分离质因子，考虑 \(i, j\) 的 唯一分解 \(i = \prod p_i ^ {k_i}, j = \prod p_i ^ {t_i}\)

于是显然有 \(ij\) 的唯一分解 \(ij = \prod p_i ^ {k_i + t_i}\)，对于 \(ij\) 的因数 \(d = \prod p_i ^ {c_i}\)

若有 \(c_i \le k_i\)，我们令 \(x\) 在 \(p_i\) 一项次数上取 \(c_i\)， \(y\) 对应次数为 \(0\)，符合互质条件

若有 \(c_i > k_i\)，我们令 \(x\) 在 \(p_i\) 一项次数上取 \(0\)，\(y\) 对应次数为 \(k_i + t_i - c_i\)，同样符合条件

显然当 \(c_i \in [0, k_i + t_i]\) 时，这样的构造形成一个 双射，且显然，每个质因子的贡献独立

故可以证得这个结论

Trick 2 - ex

根据上面的证明，我们可以容易将这个结论拓展到 \(n\) 元的情况

\[\begin {aligned} d (a_1a_2...a_n) = \sum _ {b_1 \mid a_1} { \sum _ {b_2 \mid a_2} { ... \sum _ {b_n \mid a_n} { [\gcd (b_i, b_j) = 1] ~ (i, j \in [1, n]) } } } \end {aligned} \]

同样对质因子 分开考虑，然后对 \(c_i\) 分段

我们设 \(a_j = \prod p_i ^ {k_{i, j}}\)，显然有 \(\sum _ j k_{i, j} = c_i\)，考虑设 \(s_{i, r} = \sum _ {j = 1} ^ r a_{i, j}\)（就是个 前缀和）

那么对于任意因数 \(d = \prod p_i ^ {e_i}\)

若 \(r_i \in [s_{i, j}, s_{i, j + 1}]\)，则我们令 \(a_{j + 1}\) 在 \(p_i\) 一项次数上取 \(e_i - s_{i, j}\)，其余 \(a\) 在 \(p_i\) 项次数为 \(0\)

容易证明满足 互质条件，同时 质因数间独立，这样的构造形成一个 双射，那么结论成立

为什么这样的构造是一个 双射（因数 与 取法 一一对应）

感性理解，一个质因数 \(p\)，显然只能被一个 \(a\) 包含（均不考虑 \(p ^ 0\) 的情况）

在 \(a_1 a_2 ... a_n\) 的因数中，这个 \(p\) 有 \(c\) 种可能状态（\(p, p ^ 2, ..., p_c\)）

恰好要使得一个 \(a\) 包含 \(p\)，我们仅有 \(c\) 种可能方式

即 \(a_i\) 包含 \(1, 2, ..., k_i\) 次，其余 \(a\) 不包含（\(\sum k = c\)），故我们可以 一一对应

Trick 3

\[\begin {aligned} \sum _ {i = 1} ^ n { i [\gcd (i, n) = 1] } = \dfrac {n (\varphi (n) + \epsilon (n))} 2 \end {aligned} \]

对左式 莫比乌斯反演

\[\begin {aligned} \sum _ {i = 1} ^ n { i \sum _ {d \mid i \wedge d \mid n} \mu (d) } \end {aligned} \]

转换枚举顺序

\[\begin {aligned} \sum _ {d \mid n} { \sum _ {i = 1} ^ { \frac n d } { \mu (d) id } } \end {aligned} \]

去掉一个 \(\sum\)

\[\begin {aligned} \sum _ {d \mid n} { \mu (d) d \dfrac { \frac n d \left ( \frac n d + 1 \right ) } 2 } \end {aligned} \]

再移动一下

\[\begin {aligned} \dfrac n 2 \sum _ {d \mid n} { \mu (d) \left ( \dfrac n d + 1 \right ) } \end {aligned} \]

后半部分就是一个 \(\mu * (\operatorname {id} + 1) = \varphi + \epsilon\)，于是原式就得证了

Trick 3 - ex

其实就是下面 Luogu P3700 [CQOI2017] 小 Q 的表格 这个题的 重要结论

但是这个式子确实 十分简洁，也比较常见，就在这里记录一下

\[\begin {aligned} \sum _ {i = 1} ^ n { \sum _ {j = 1} ^ n { [\gcd (i, j) = 1] ij } } = \sum _ {i = 1} ^ n { i ^ 2 \varphi (i) } \end {aligned} \]

显然 \(\gcd (a, b) = \gcd (b, a)\)，于是我们只需要维护一半的内容

\[\begin {aligned} 2 \sum _ {i = 1} ^ n { \sum _ {j = 1} ^ i { [\gcd (i, j) = 1] ij } } - \sum _ {i = 1} ^ n [i = 1] \end {aligned} \]

然后使用 Trick 3，带入到 \(\sum _ {j = 1} ^ i [\gcd (i, j) = 1]ij\) 里面

\[\begin {aligned} \left ( 2 \sum _ {i = 1} ^ n { i \dfrac {i (\varphi (i) + \epsilon (i))} 2 } \right ) - 1 \end {aligned} \]

简单化简一下，注意到 \(\epsilon (i) = [i = 1]\)，于是得到证明

4. 例题

巨大多多多多题，可能有各种技巧 / 算法的复合

由于例题过多，部分存在于 [0] 但 与其他问题相似 / 被包含 的问题就被 跳过了

可能后面有时间会补

Luogu P2522 [HAOI2011] Problem b

直接形式化题意是这样的

\[\begin {aligned} \sum _ {x = a} ^ b { \sum _ {y = c} ^ d { [\gcd (x, y) = k] } } \end {aligned} \]

多次询问上述式子的值，每次给定 \(a, b, c, d, k\)

乍一看感觉这东西 啥也不是，于是我们进行简单转化，显然，这个式子可以差分

于是我们只需考虑形如 \(\sum _ {x = 1} ^ n \sum _ {y = 1} ^ m [\gcd (x, y) = k]\) 的式子 怎么求即可

集中注意力，我们发现它和 Trick 1 十分相似，只有 \(\gcd\) 值的区别，容易想到上界同时除 \(k\)

\[\begin {aligned} & \sum _ {x = 1} ^ n { \sum _ {y = 1} ^ m { [\gcd (x, y) = k] } } \\ =& \sum _ {x = 1} ^ { \left \lfloor \frac n k \right \rfloor } { \sum _ {y = 1} ^ { \left \lfloor \frac m k \right \rfloor } { [\gcd (x, y) = 1] } } \end {aligned} \]

于是应用 Trick 1，就可以转化成下面的最终形式

\[\begin {aligned} \sum _ {d = 1} ^ { \min \left ( \left \lfloor \frac n k \right \rfloor , \left \lfloor \frac m k \right \rfloor \right ) } \mu (d) \left \lfloor \dfrac n {kd} \right \rfloor \left \lfloor \dfrac m {kd} \right \rfloor \end {aligned} \]

然后可以 整除分块 做掉，时间复杂度 \(O (n + T \sqrt n)\)

#include <iostream>
#include <cstdint>

const int MAXN = 50005;

int mu[MAXN], pri[MAXN], Cnt = 0;
bool vis[MAXN];

inline void Sieve (const int n = 50000) {
	mu[1] = 1;
	for (int i = 2; i <= n; ++ i) {
		if (!vis[i]) pri[++ Cnt] = i, mu[i] = - 1;
		for (int k = 1; k <= Cnt && pri[k] * i <= n; ++ k) {
			vis[pri[k] * i] = 1;
			if (i % pri[k] == 0) break ;
			mu[pri[k] * i] = - mu[i];
		}
	}
	for (int i = 1; i <= n; ++ i) mu[i] += mu[i - 1];
}

int a, b, c, d, k;

inline int64_t Calc (int n, int m) {
	int64_t d, v, ans = 0; n /= k, m /= k;
	for (int l = 1, r; l <= std::min (n, m); l = r + 1) {
		r = std::min (m / (m / l), n / (n / l)), d = m / l, v = n / l;
		ans += d * v * (mu[r] - mu[l - 1]);
	}
	return ans;
}

int N;

int main () {
	
	std::cin >> N, Sieve ();
	
	for (int i = 1; i <= N; ++ i) {
		std::cin >> a >> b >> c >> d >> k;
		std::cout << Calc (b, d) + Calc (a - 1, c - 1) - Calc (a - 1, d) - Calc (b, c - 1) << '\n';
	}
	
	return 0;
}

EXP - 1 - Luogu P3455 [POI2007] ZAP-Queries（\(a = c = 0\)）

Luogu P2257 YY的GCD

形式化题意是这样的

\[\begin {aligned} \sum _ {x = 1} ^ N { \sum _ {y = 1} ^ M { [\gcd (x, y) \in \mathbb P] } } \end {aligned} \]

多次询问，每次给定 \(N, M\)

这个 \([\gcd (x, y) \in \mathbb P]\) 这部分不好 整体刻画，考虑 枚举质数

\[\begin {aligned} \sum _ {p \in \mathbb P} { \sum _ {x = 1} ^ N { \sum _ {y = 1} ^ M { [\gcd (x, y) = p] } } } \end {aligned} \]

然后 套路上界除以 \(p\)，得到一个比较好的形式

\[\begin {aligned} \sum _ {p \in \mathbb P} { \sum _ {x = 1} ^ { \left \lfloor \frac N p \right \rfloor } { \sum _ {y = 1} ^ { \left \lfloor \frac M p \right \rfloor } { [\gcd (x, y) = 1] } } } \end {aligned} \]

使用 Trick 1，转化

\[\begin {aligned} \sum _ {p \in \mathbb P} { \sum _ {d = 1} ^ { \min \left ( \left \lfloor \frac N p \right \rfloor , \left \lfloor \frac M p \right \rfloor \right ) } { \mu (d) \left \lfloor \dfrac N {pd} \right \rfloor \left \lfloor \dfrac M {pd} \right \rfloor } } \end {aligned} \]

现在时间复杂度是 \(O (TN \ln N)\) 的，仍然无法通过此题

显然，这时候瓶颈就在于 枚举质数 的部分了，考虑怎么 消除其影响，我们设 \(k = pd\)

于是上式可以转化成

\[\begin {aligned} \sum _ {p \in \mathbb P} { \sum _ {d = 1} ^ { \min \left ( \left \lfloor \frac N p \right \rfloor , \left \lfloor \frac M p \right \rfloor \right ) } { \mu (d) \left \lfloor \dfrac N k \right \rfloor \left \lfloor \dfrac M k \right \rfloor } } \end {aligned} \]

更进一步的，交换枚举顺序后，我们有

\[\begin {aligned} \sum _ {k = 1} ^ { \min (n, m) } { \left \lfloor \dfrac N k \right \rfloor \left \lfloor \dfrac M k \right \rfloor \sum _ {p \mid k, p \in \mathbb P} { \mu \left ( \left \lfloor \dfrac k p \right \rfloor \right ) } } \end {aligned} \]

于是前半部分直接 整除分块 就可以在 \(O (T \sqrt N)\) 的时间内完成

后半部分看上去是 困难的，但事实上我们在 线性筛 中可以 预处理 其值，计算 前缀和

（看成一个函数 \(f (x) = \sum _ {p \mid x, p \in \mathbb P} \mu \left ( \left \lfloor \dfrac x p \right \rfloor \right)\)）

注意 线性筛 中 i % pri[k] == 0 的特殊处理

考虑此时 i * pri[k] 中一定有 pri[k] ^ 2 这个 平方因数

故 \(\forall p \neq pri_k, \mu \left ( \left \lfloor \dfrac x p \right \rfloor \right) = 0\)，只需考虑 \(p = pri_k\) 的情况

即 \(\mu (i)\) 的值

然后就可以 获取胜利 了，时间复杂度 \(O (N + T \sqrt N)\)

#include <iostream>
#include <cstdint>

const int MAXN = 10000005;

int pri[MAXN], mu[MAXN], ff[MAXN], Cnt = 0;
bool vis[MAXN];

inline void Sieve (const int n = MAXN - 5) {
	mu[1] = 1;
	for (int i = 2; i <= n; ++ i) {
		if (!vis[i]) pri[++ Cnt] = i, mu[i] = - 1, ff[i] = 1;
		for (int k = 1; k <= Cnt && pri[k] * i <= n; ++ k) {
			vis[pri[k] * i] = 1;
			if (i % pri[k] == 0) {
				ff[pri[k] * i] = mu[i];
				break ;
			}
			mu[pri[k] * i] = - mu[i];
			ff[pri[k] * i] = - ff[i] + mu[i];
		}
	}
	for (int i = 1; i <= n; ++ i) ff[i] += ff[i - 1];
}

inline int64_t Calc (const int n, const int m) {
	int64_t ans = 0, d, v;
	for (int l = 1, r; l <= std::min (n, m); l = r + 1) {
		r = std::min (n / (n / l), m / (m / l)), d = n / l, v = m / l;
		ans += d * v * (ff[r] - ff[l - 1]);
	}
	return ans;
}

int T;
int a, b;

int main () {
	
	std::cin >> T, Sieve ();
	
	for (int i = 1; i <= T; ++ i) {
		std::cin >> a >> b;
		std::cout << Calc (a, b) << '\n';
	}
	
	return 0;
}

EXP - 1 - Luogu P2568 GCD

• ##### SP5971 LCMSUM - LCM Sum

Luogu P4318 完全平方数

考虑设 \(f (n)\) 表示 \([1, n]\) 中有多少 小 \(X\) 不讨厌的数，显然其具有 单调性

可以考虑 求出来之后二分答案，考虑怎么快速求解

于是考虑每个 质数 \(p\)，容易发现我们需要 除去 \([1, n]\) 中所有 整除 \(p_i ^ 2\) 的数

而在此之后，显然整除 \((p_ip_j) ^ 2\) 的数 多算了一次，需要 加回来

以此类推，容易发现这样 容斥系数 就是 \(\mu\)，于是有下面的计算式

\[\begin {aligned} f (n) = \sum _ {i = 1} ^ {\sqrt n} { \mu (i) \left \lfloor \dfrac n {i ^ 2} \right \rfloor } \end {aligned} \]

直接做就可以了，不需要整除分块，套上二分，时间复杂度 \(O (\sqrt n \log n)\)

我们的 \(\mu\) 预处理到 \(\sqrt n\) 就行了（保证 \(K\) 能取到 \(10 ^ 9\) 的 \(n\)），时间正确

#include <bits/stdc++.h>

const int MAXN = 65540;

int pri[MAXN], mu[MAXN], cnt;
bool vis[MAXN];

inline void Sieve (const int n = 1 << 16) {
	mu[1] = 1;
	for (int i = 2; i <= n; ++ i) {
		if (!vis[i]) pri[++ cnt] = i, mu[i] = - 1;
		for (int k = 1; k <= cnt && pri[k] * i <= n; ++ k) {
			vis[i * pri[k]] = 1;
			if (i % pri[k] == 0) break ;
			mu[i * pri[k]] = - mu[i];
		}
	}
}

inline int  F (const int n) {
	int ans = 0;
	for (int i = 1; i <= sqrt (n); ++ i) ans += mu[i] * (n / i / i);
	return ans;
}

inline int  G (const int n) {
	int l = 0, r = INT_MAX, m = 0, ans = 0;
	while (r > l) {
		m = (0ll + l + r) >> 1;
		if (F (m) >= n) r = ans = m;
		else l = m + 1;
	}
	return ans;
}

int N, k;

int main () {
	
	std::cin >> N, Sieve ();
	
	for (int i = 1; i <= N; ++ i) {
		std::cin >> k;
		std::cout << G (k) << '\n';
	}
	
	return 0;
}

• Luogu P1829 [国家集训队] Crash的数字表格 / JZPTAB

• Luogu P3704 [SDOI2017] 数字表格

• Luogu P5221 Product

• Luogu P3911 最小公倍数之和

• [AGC038C] LCMs

• CF1285F Classical?

Luogu P3327 [SDOI2015] 约数个数和

题面好像就挺 形式化 的

\[\begin {aligned} \sum _ {i = 1} ^ n { \sum _ {j = 1} ^ m { d (ij) } } \end {aligned} \]

多次询问上式的取值，每次给出 \(n, m\)

\(d (ij)\) 我们并没有什么其它很好的转化，于是直接使用 Trick 2 将其转化成 下式

\[\begin {aligned} \sum _ {i = 1} ^ n { \sum _ {j = 1} ^ m { \sum _ {x \mid i} { \sum _ {y \mid j} { [x \perp y] } } } } \end {aligned} \]

显然 \([x \perp y]\) 等价于 \([\gcd (x, y) = 1]\)，于是套用 发现时的式子的变式，转化到

\[\begin {aligned} \sum _ {i = 1} ^ n { \sum _ {j = 1} ^ m { \sum _ {x \mid i} { \sum _ {y \mid j} { \sum _ {d \mid \gcd (x, y)} \mu (d) } } } } \end {aligned} \]

改变枚举顺序，容易得到等价形式

\[\begin {aligned} \sum _ {d = 1} ^ {\min (n, m)} { \mu (d) \sum _ {i = 1} ^ n { \sum _ {j = 1} ^ m { \sum _ {x \mid i} { \sum _ {y \mid j} { [d \mid \gcd (x, y)] } } } } } \end {aligned} \]

注意到后面的 \(\gcd\) 只与 \(x, y\) 有关

显然枚举 \(i, j\) 并没有优势，于是枚举 \(x, y\)，同时考虑每个 \(x, y\) 会被多少 \(i, j\) 枚举到

\[\begin {aligned} \sum _ {d = 1} ^ {\min (n, m)} { \mu (d) \sum _ {x = 1} ^ n { \sum _ {y = 1} ^ m { [d \mid \gcd (x, y)] \left \lfloor \dfrac n x \right \rfloor \left \lfloor \dfrac m y \right \rfloor } } } \end {aligned} \]

通过枚举 \(d\) 的倍数来省去 \([d \mid \gcd (x, y)]\) 这个 令人火大的条件

\[\begin {aligned} \sum _ {d = 1} ^ {\min (n, m)} { \mu (d) \sum _ {x = 1} ^ { \left \lfloor \frac n d \right \rfloor } { \sum _ {x = 1} ^ { \left \lfloor \frac m d \right \rfloor } { \left \lfloor \dfrac n {dx} \right \rfloor \left \lfloor \dfrac m {dy} \right \rfloor } } } \end {aligned} \]

容易发现，后面部分两个和式可以变成 同构的部分

\[\sum _ {d = 1} ^ {\min (n, m)} { \mu (d) \left ( \sum _ {x = 1} ^ { \left \lfloor \frac n d \right \rfloor } { \left \lfloor \dfrac n {dx} \right \rfloor } \right ) \left ( \sum _ {y = 1} ^ { \left \lfloor \frac m d \right \rfloor } { \left \lfloor \dfrac m {dy} \right \rfloor } \right ) } \]

显然，对于一段使得 \(\left \lfloor \dfrac n d \right \rfloor\) 相等的 \(d\)，\(g (\left \lfloor \frac n d \right \rfloor) = \left ( \sum \limits _ {x = 1} ^ {\left \lfloor \frac n d \right \rfloor} \left \lfloor \dfrac n {dx} \right \rfloor \right)\) 的值相等，后半部分同理

于是我们可以枚举 \(\left \lfloor \dfrac n d \right \rfloor\) 的值，预处理出对应部分的取值，由于上式也符合 整除分块 的结构

故对于固定 \(\left \lfloor \dfrac n d \right \rfloor\)，预处理时间是 \(O (\sqrt {\left \lfloor \dfrac n d \right \rfloor})\) 的，预处理总复杂度 \(O (n \sqrt n)\)

之后每组询问同样 整除分块（对 \(d\) 分段） 做就行了，这部分复杂度 \(O (T \sqrt n)\)

总时间复杂度 \(O (n \sqrt n + T \sqrt n) = O (n \sqrt n)\)

#include <iostream>
#include <cstdint>

const int MAXN = 50005;

int64_t f[MAXN];
int pri[MAXN], mu[MAXN], cnt = 0;
bool vis[MAXN];

inline void Sieve (const int n = 50000) {
	mu[1] = 1;
	for (int i = 2; i <= n; ++ i) {
		if (!vis[i]) pri[++ cnt] = i, mu[i] = - 1;
		for (int k = 1; k <= cnt && pri[k] * i <= n; ++ k) {
			vis[pri[k] * i] = 1;
			if (i % pri[k] == 0) break ;
			mu[i * pri[k]] = - mu[i];
		}
	}
	for (int i = 1; i <= n; ++ i) mu[i] += mu[i - 1];
	for (int i = 1; i <= n; ++ i) {
		int64_t ans = 0;
		for (int l = 1, r; l <= i; l = r + 1) 
			r = i / (i / l), ans += (r - l + 1) * (i / l);
		f[i] = ans;
	}
}

inline int64_t Calc (const int n, const int m) {
	int64_t ans = 0;
	for (int l = 1, r; l <= std::min (n, m); l = r + 1)
		r = std::min (n / (n / l), m / (m / l)), ans += f[n / l] * f[m / l] * (mu[r] - mu[l - 1]);
	return ans;
}

int T, a, b;

int main () {
	
	std::cin >> T, Sieve ();
	
	for (int i = 1; i <= T; ++ i) {
		std::cin >> a >> b;
		std::cout << Calc (a, b) << '\n';
	}
	
	return 0;
}

Luogu P4619 [SDOI2018] 旧试题

看上去很像 上个题 的加强版

\[\begin {aligned} \left ( \sum _ {i = 1} ^ A { \sum _ {j = 1} ^ B { \sum _ {k = 1} ^ C { d (ijk) } } } \right ) \bmod (10 ^ 9 + 7) \end {aligned} \]

运用 Trick 2 - ex，我们可以将之转化为下式

\[\begin {aligned} \sum _ {i = 1} ^ A { \sum _ {j = 1} ^ B { \sum _ {k = 1} ^ C { \sum _ {x \mid i} { \sum _ {y \mid j} { \sum _ {z \mid k} { [\gcd (i, j) = 1][\gcd (j, k) = 1][\gcd (i, k) = 1] } } } } } } \end {aligned} \]

注意到 \([\gcd (i, j) = 1][\gcd (j, k) = 1][\gcd (i, k) = 1]\) 与 \([\gcd (i, j, k) = 1]\) 不等价

与上面的题相同的套路，我们容易的将式子转化成下面的样子

\[\begin {aligned} \sum _ {x = 1} ^ A { \sum _ {y = 1} ^ B { \sum _ {z = 1} ^ C { \left \lfloor \dfrac A x \right \rfloor \left \lfloor \dfrac B y \right \rfloor \left \lfloor \dfrac C z \right \rfloor \sum _ { d_1 \mid x \wedge d_1 \mid y } \mu (d_1) \sum _ { d_2 \mid x \wedge d_2 \mid z } \mu (d_2) \sum _ { d_3 \mid y \wedge d_3 \mid z } \mu (d_3) } } } \end {aligned} \]

然后变换枚举顺序， 注意 整除条件的取，有

\[\sum _ {d_1 = 1} ^ {\min (A, B)} { \mu (d_1) \sum _ {d_2 = 1} ^ {\min (A, C)} { \mu (d_2) \sum _ {d_3 = 1} ^ {\min (B, C)} { \mu (d_3) \sum _ { d_1 \mid x \wedge d_2 \mid x } { \left \lfloor \dfrac A x \right \rfloor } \sum _ { d_1 \mid y \wedge d_3 \mid y } { \left \lfloor \dfrac B y \right \rfloor } \sum _ { d_2 \mid z \wedge d_3 \mid z } { \left \lfloor \dfrac C z \right \rfloor } } } } \]

同样的套路，注意到 \(d_1 \mid x \wedge d_2 \mid x \Rightarrow \operatorname {lcm} (d_1, d_2) \mid x\)，于是最终变成下式

\[\sum _ {d_1 = 1} ^ {\min (A, B)} { \mu (d_1) \sum _ {d_2 = 1} ^ {\min (A, C)} { \mu (d_2) \sum _ {d_3 = 1} ^ {\min (B, C)} { \mu (d_3) \sum _ {x = 1} ^ { \left \lfloor \frac A {\operatorname {lcm} (d_1, d_2)} \right \rfloor } { \left \lfloor \frac A {x \operatorname {lcm} (d_1, d_2)} \right \rfloor } \sum _ {y = 1} ^ { \left \lfloor \frac B {\operatorname {lcm} (d_1, d_3)} \right \rfloor } { \left \lfloor \frac B {y \operatorname {lcm} (d_1, d_3)} \right \rfloor } \sum _ {z = 1} ^ { \left \lfloor \frac C {\operatorname {lcm} (d_2, d_3)} \right \rfloor } { \left \lfloor \frac C {z \operatorname {lcm} (d_2, d_3)} \right \rfloor } } } } \]

显然，后面的部分又可以看作 \(g (x)\) 进行预处理，只考虑 前面部分怎么做

暴力枚举 \(O (n ^ 3)\)，反演就 白做了，仔细思考，发现 有贡献的三元组有限

具体而言，可能有贡献 当且仅当 \(\mu (x) \neq 0\) 且 \(\operatorname {lcm} \le \max (A, B, C)\)

为什么 \(\operatorname {lcm}\) 只需要限制小于 \(\max (A, B, C)\) 就行，不需要分成三组限制

分成三组 不便于统计，而所有不满足实际限制的，其 \(g (x)\) 部分一定为 \(0\)，故 贡献是对的

于是对于每个满足上述条件的 二元组，我们将其 连边

尝试一下，发现这样的边数至多 \(821535\)（一说 \(760741\)）条，好像很可做

暴力连边是 \(O (n ^ 2)\) 的，我们也不是很喜欢，考虑优化

仍然枚举 \(\gcd\)，然后在一个 \(\gcd\) 下枚举两数（同时 \(\operatorname {lcm} \le \max (A, B, C)\) 限制上界）

这样的话连边就是 \(O (m)\) 的了

显然构建的图中的每个 三元环 就代表了一组 符合条件的三元组

于是使用 神秘三元环计数 的套路，以 度数 为第一关键字，编号 为第二关键字 从大到小 枚举

可以把枚举复杂度降到 \(O (m \sqrt m)\)，注意到这样的三元环计数并不考虑 三个点的顺序

所以最后我们还需要对找出的每组 三元环三个点 交换顺序计算贡献，然后就做完了

时间复杂度是 \(O (T (n \ln n + m \sqrt m))\) 的

#include <iostream>
#include <cstdint>
#include <vector>
#include <unordered_map>
#include <algorithm>

const int MAXN = 100005;
const int MOD  = 1000000007;

int64_t f[MAXN];
int pri[MAXN], mu[MAXN], cnt = 0;
bool vis[MAXN];

inline void Sieve (const int n = 100000) {
	mu[1] = 1;
	for (int i = 2; i <= n; ++ i) {
		if (!vis[i]) pri[++ cnt] = i, mu[i] = - 1;
		for (int k = 1; k <= cnt && pri[k] * i <= n; ++ k) {
			vis[pri[k] * i] = 1;
			if (i % pri[k] == 0) break ;
			mu[i * pri[k]] = - mu[i];
		}
	}
	for (int i = 1; i <= n; ++ i) {
		int64_t ans = 0;
		for (int l = 1, r; l <= i; l = r + 1) 
			r = i / (i / l), ans += (r - l + 1) * (i / l);
		f[i] = ans;
	}
}

struct Node {
	int u, v, w;
} E[MAXN << 4];

int32_t Deg[MAXN], Vis[MAXN];
int32_t Max = 0, tot = 0;

std::vector <Node> G[MAXN];
std::unordered_map <int, bool> Mp[MAXN];


inline int64_t RetAns (const int a, const int b, const int c, const int x, const int y, const int z) {
	return (1ll * f[a / x] * f[b / y] % MOD * f[c / z] % MOD) % MOD;
}

int T, A, B, C;

inline int64_t GetAns (const int d1, const int d2, const int d3, const int lcm1, const int lcm2, const int lcm3) {
	int64_t ans = 0;
	if (d1 == d2 && d2 == d3) ans += RetAns (A, B, C, lcm1, lcm2, lcm3);
	else if (d1 == d2 || d2 == d3 || d1 == d3) ans += RetAns (A, B, C, lcm1, lcm2, lcm3), ans += RetAns (C, A, B, lcm1, lcm2, lcm3), ans += RetAns (B, C, A, lcm1, lcm2, lcm3);
	else ans += RetAns (A, B, C, lcm1, lcm2, lcm3), ans += RetAns (A, C, B, lcm1, lcm2, lcm3), ans += RetAns (B, C, A, lcm1, lcm2, lcm3), ans += RetAns (B, A, C, lcm1, lcm2, lcm3), ans += RetAns (C, A, B, lcm1, lcm2, lcm3), ans += RetAns (C, B, A, lcm1, lcm2, lcm3);
	ans = ans % MOD, ans = (mu[d1] * mu[d2] * mu[d3] * ans % MOD + MOD) % MOD;
	return ans;
}

inline int64_t calc () {
	int64_t ans = 0;
	
	for (int i = 1; i <= tot; ++ i) {
		__builtin_prefetch (& E[i]);
		if (Deg[E[i].u] == Deg[E[i].v] ? E[i].u > E[i].v : Deg[E[i].u] > Deg[E[i].v])
			std::swap (E[i].u, E[i].v);
		G[E[i].u].push_back (E[i]);
	}
	
	for (int i = 1; i <= Max; ++ i) {
		
		for (auto x : G[i]) Vis[x.v] = x.w;
		for (auto x : G[i])
			for (auto k : G[x.v])
				if (Vis[k.v])
					ans += GetAns (i, x.v, k.v, Vis[k.v], x.w, k.w), ans %= MOD;
		for (auto x : G[i]) Vis[x.v] = 0;
	}
	
	return ans;
}

inline int64_t Calc (const int n, const int m, const int k) {
	int64_t ans = 0;
	Max = std::max ({n, m, k}), tot = 0;
	
	for (int i = Max; i >= 1; -- i) {
		if (mu[i]) {
			for (int j = i; j <= Max; j += i) {
				if (mu[j]) {
					for (int k = j; k <= 1ll * Max * i / j; k += i) {
						if (mu[k]) {
							if (!Mp[j][k]) ++ tot, E[tot] = {j, k, j / i * k}, ++ Deg[j], ++ Deg[k], Mp[j][k] = 1;
						}
					}
				}
			}
		}
	}
	
	ans += calc ();
	
	for (int i = 1; i <= Max; ++ i) G[i].clear (), Mp[i].clear (), Deg[i] = 0;
	return (ans + MOD) % MOD;
}

int main () {
	
	std::cin >> T, Sieve ();
	
	for (int i = 1; i <= T; ++ i) {
		std::cin >> A >> B >> C;
		std::cout << Calc (A, B, C) << '\n';
	}
	
	return 0;
}

EXP - 1 - CF235E Number Challenge

EXP - 2 - CF236B Easy Number Challenge

Luogu P5518 [MtOI2019] 幽灵乐团 / 莫比乌斯反演基础练习题

👻🎶👻👻👻👻🎶🎶🎶🎶🎶👻🎶🎶🎶🎶👻🎶🎶🎶🎶🎶🎶🎶🎶🎶👻🎶🎶🎶🎶🎶🎶🎶🎶🎶👻👻👻👻👻👻👻👻🎶🐨

这个题在赛时怎么是评的 紫 啊，逆大天

看这玩意儿，三重循环，还是 \(\prod\)，自己做的心已经死了

做过 Luogu P5221 Product 的会发现 眼前这个题又是一个 加强版，但可惜我 没有做过

分了三个 \(Type\)，显然 \(Type ~ 0\) 最为简单，考虑从这里下手

\(Type ~ 0\)

\[\begin {aligned} \prod _ {i = 1} ^ A { \prod _ {j = 1} ^ B { \prod _ {k = 1} ^ C { \left ( \dfrac {\operatorname {lcm} (i, j)} {\gcd (i, k)} \right ) } } } \end {aligned} \]

\(O (n ^ {0.75} \log n)\)

我们 \(\operatorname {lcm}, \gcd\) 一上一下并不好搞，考虑将 \(\operatorname {lcm}\) 表示成 \(\dfrac x {\gcd}\) 的形式，让 分子干净一点

\[\begin {aligned} \prod _ {i = 1} ^ A { \prod _ {j = 1} ^ B { \prod _ {k = 1} ^ C { \dfrac {ij} {\gcd (i, j) \gcd (i, k)} } } } \end {aligned} \]

然后显然 分子这部分可以单独搞掉，就可以把分子化成 \(1\)

\[\begin {aligned} \left ( \prod _ {i = 1} ^ A i \right ) ^ {BC} \left ( \prod _ {j = 1} ^ B j \right ) ^ {AC} \dfrac 1 { \left ( \prod _ {i = 1} ^ A { \prod _ {j = 1} ^ B { \gcd (i, j) } } \right ) ^ C \left ( \prod _ {i = 1} ^ A { \prod _ {k = 1} ^ C { \gcd (i, k) } } \right ) ^ B } \end {aligned} \]

显然，\(\left ( \prod _ {i = 1} ^ A i \right ) ^ {BC}\)和 \(\left ( \prod _ {j = 1} ^ B j \right ) ^ {AC}\) 用 快速幂 加 预处理阶乘 就可以在 \(\log\) 时间简单做掉

现在的难点在于处理 下面的式子（即 分母 上两个 同构 的部分）

\[\begin {aligned} \prod _ {i = 1} ^ n { \prod _ {j = 1} ^ m { \gcd (i, j) } } \end {aligned} \]

只要能在 合理的时间复杂度内 做掉这部分，那么 \(Type ~ 0\) 就结束了

这时候，做过 Luogu P3704 [SDOI2017] 数字表格 的就会发现这部分和那个题 本质相同

但是我还是没有做过，乐，于是只能 开始推神秘狮子

和做 \(\sum\) 时的想法一样，仍然考虑枚举 \(\gcd\) 的值 \(d\)，于是得到

\[\begin {aligned} \prod _ {d = 1} ^ {\min (n, m)} { d ^ { \sum _ {i = 1} ^ { \left \lfloor \frac A d \right \rfloor } \sum _ {j = 1} ^ { \left \lfloor \frac B d \right \rfloor } [\gcd (i, j) =1] } } \end {aligned} \]

然后到了经典 Trick 1 时间，套上就行

\[\begin {aligned} \prod _ {d = 1} ^ {\min (n, m)} { d ^ { \sum _ {k = 1} ^ { \min \left ( \left \lfloor \frac A d \right \rfloor , \left \lfloor \frac B d \right \rfloor \right ) } \mu (k) \left \lfloor \frac A {dk} \right \rfloor \left \lfloor \frac B {dk} \right \rfloor } } \end {aligned} \]

\(\sum _ {k = 1} ^ {\min \left ( \left \lfloor \frac A k \right \rfloor, \left \lfloor \frac B k \right \rfloor \right)} \mu (k) \left \lfloor \frac A {dk} \right \rfloor \left \lfloor \frac B {dk} \right \rfloor\) 的计算显然是 \(O (\sqrt n)\) 的

接着，显然对于一段 \(d\), \(\left \lfloor \dfrac A d \right \rfloor\)，\(\left \lfloor \dfrac B d \right \rfloor\) 的值也是 相等 的，于是这里也可以 整除分块

然后我们就可以 \(O (n ^ {0.75} \log n)\) 解决这个部分

于是此时 \(Type ~ 0\) 已经可以做了，时间复杂度 \(O (Tn ^ {0.75} \log n)\)

这部分的代码

namespace Type_0_1 {
	
	int64_t f[MAXN], rf[MAXN];
	
	inline void Init (const int n = 100000) {
		f[0] = rf[0] = 1;
		for (int i = 1; i <= n; ++ i) f[i] = f[i - 1] * i % MOD;
		for (int i = 1; i <= n; ++ i) rf[i] = qpow (f[i]);
	}
	
	inline int64_t calc (const int n, const int m) {
		int64_t ans = 0, d = 0, v = 0;
		for (int l = 1, r; l <= std::min (n, m); l = r + 1) {
			r = std::min (n / (n / l), m / (m / l)), d = n / l, v = m / l;
			ans += 1ll * (s[r] - s[l - 1]) * d * v, ans = (ans % (MOD - 1) + (MOD - 1)) % (MOD - 1);
		}
		return ans % (MOD - 1);
	}
	
	inline int64_t Calc (const int n, const int m) {
		int64_t ans = 1, d = 0, v = 0;
		for (int l = 1, r; l <= std::min (n, m); l = r + 1) {
			r = std::min (n / (n / l), m / (m / l)), d = n / l, v = m / l;
			ans = ans * qpow (f[r] * rf[l - 1] % MOD, calc (d, v)) % MOD;
		}
		return qpow (ans); // NOTE 1 / ans
	}
	
	inline int64_t Solve (const int64_t a, const int64_t b, const int64_t c) {
		return qpow (f[a], b * c) * qpow (f[b], a * c) % MOD * qpow (Calc (a, b), c) % MOD * qpow (Calc (a, c), b) % MOD;
	}
	
}

\(f (n) = n !\)

这个东西的时间复杂度证明十分有趣，稍不注意就可能会当成 \(O (n)\) 甚至 \(O (n \log n)\) 的了

注意到我们 内层整除分块 calc1 由于每次用了 快速幂，所以时间是 \(O (\sqrt {n'} \log n)\) 的

而 外层整除分块 calc2 调用了 \(O (\sqrt n)\) 次 calc1，调用完之后快速幂

所以外层这玩意儿时间是 \(O (\sqrt n (\sqrt {n'} \log n + \log n))\) 的，重点就在于这个 \(n'\) 了

这个 \(n'\) 实际上等于式子里的 \(\left \lfloor \dfrac A d \right \rfloor, \left \lfloor \dfrac B d \right \rfloor\)，也就是程序里的 d = n / l， v = m / l

那么这东西取值一共 \(O (\sqrt n)\) 个，其中值 小于 \(\sqrt n\) 的有 \(\sqrt n\) 个

对于这部分，\(\max n' = \sqrt n\)，复杂度 \(O (\sqrt n (\sqrt {\sqrt n} \log n + \log n)) = O (n ^ {0.75} \log n)\)

而对于值大于 \(\sqrt n\) 的 \(\sqrt n\) 个东西，\(n' = \dfrac n i ~ (i \in [1, \sqrt n])\)

于是时间复杂度是 \(O (\sum _ {i = 1} ^ {\sqrt n} {\sqrt {\dfrac n i}}\log n) = O (\sqrt n \sum _ {i = 1} ^ {\sqrt n} \dfrac 1 {\sqrt i} \log n)\) 的

而 \(\sum _ {i = 1} ^ m {\dfrac 1 {\sqrt i}} = \dfrac {\pi ^ 4} {90} + 2 \sqrt n + \dfrac 1 {2 \sqrt n} - \dfrac 1 {24n ^ {1.5}} \sim O (\sqrt m)\)，故 \(\sum _ {i = 1} ^ {\sqrt n} \dfrac 1 {\sqrt {i}} \sim O (n ^ {0.25})\)

故这一部分时间复杂度也是 \(O (n ^ {0.75} \log n)\) 的，总时间复杂度就是 \(O (n ^ {0.75} \log n)\) 的了

\(O (\sqrt n \log n)\)

看到这个 Type_0_1，有人知道不对劲了，有什么没有完成的部分呢？

这东西 还不够快，特别是当有人告诉你 后面还要用 的时候，于是我们得 再进一步

我们把 \(k\) 拿下来，可以得到 下面的式子

\[\begin {aligned} \prod _ {d = 1} ^ {\min (n, m)} { \prod _ {k = 1} ^ {\min (n, m)} { d ^ { \mu (k) \left \lfloor \frac n {kd} \right \rfloor \left \lfloor \frac m {kd} \right \rfloor } } } \end {aligned} \]

为啥拿下来之后 上界就变了？

额，其实 有意义的上界 是没变的，只是这样写方便后续运算，而且

当 \(k > \min (\left \lfloor \dfrac n d \right \rfloor, \left \lfloor \dfrac m d \right \rfloor)\) 时，\(\left \lfloor \dfrac n {kd} \right \rfloor\)，\(\left \lfloor \dfrac m {kd} \right \rfloor\) 必有一个为 \(0\)，故 不做贡献，值是不变的

然后经典的令 \(T = kd\)，得到

\[\begin {aligned} \prod _ {d = 1} ^ {\min (n, m)} { \prod _ {k = 1} ^ {\min (n, m)} { d ^ { \mu \left ( \frac T d \right ) \left \lfloor \frac n T \right \rfloor \left \lfloor \frac m T \right \rfloor } } } \end {aligned} \]

我们发现 \(k\) 没存在感了，于是可以 转变枚举顺序 了，我们枚举 \(T\)，然后枚举 \(d\) 作为因数

\[\begin {aligned} \prod _ {T = 1} ^ {\min (n, m)} { \left ( \prod _ {d \mid T} { d ^ { \mu \left ( \frac T d \right ) } } \right ) ^ { \left \lfloor \frac n T \right \rfloor \left \lfloor \frac m T \right \rfloor } } \end {aligned} \]

这个就很牛了，我们枚举 \(d\)，可以容易的在 \(O (n \log n)\) 的 时间内预处理 \(\prod _ {d \mid T} d ^ {\mu \left ( \frac T d \right)}\) 这一部分

然后同样的 整除分块，我们就能在 \(O (\sqrt n \log n)\) 的时间内解决这个 \(Type ~ 0\) 的问题

注意 指数上的取模是对 MOD - 1 而非 MOD，然后 不要把样例的模数看错

哦对，n / l * m / l != (n / l) * (m / l)，最后是这部分代码

namespace Type_0_2 {
	int64_t f[MAXN], rf[MAXN], g[MAXN], rg[MAXN];
	
	inline void Init (const int n = 100000) {
		g[0] = 1, f[0] = rf[0] = 1;
		for (int i = 1; i <= n; ++ i) g[i] = g[i - 1] * i % MOD;
		for (int i = 1; i <= n; ++ i) rg[i] = qpow (g[i]);
		
		for (int i = 1; i <= n; ++ i) f[i] = 1;
		for (int i = 1; i <= n; ++ i) 
			for (int j = 1, k = i; k <= n; ++ j, k = i * j) {
				if (mu[j] == + 1) f[k] = f[k] * i % MOD;
				if (mu[j] == - 1) f[k] = f[k] * rg[i] % MOD * g[i - 1] % MOD;
			}
		for (int i = 2; i <= n; ++ i) f[i] = f[i] * f[i - 1] % MOD;
		for (int i = 1; i <= n; ++ i) rf[i] = qpow (f[i]);
	}
	
	inline int64_t Calc (const int n, const int m) {
		int64_t ans = 1, d = 0, v = 0;
		for (int l = 1, r; l <= std::min (n, m); l = r + 1) {
			r = std::min (n / (n / l), m / (m / l)), d = n / l, v = m / l;
			ans = ans * qpow (f[r] * rf[l - 1] % MOD, d * v % (MOD - 1)) % MOD;
		}
		return qpow (ans); // NOTE 1 / ans
	}
	
	inline int64_t Solve (const int64_t a, const int64_t b, const int64_t c) {
		return qpow (g[a], b * c) * qpow (g[b], a * c) % MOD * qpow (Calc (a, b), c) % MOD * qpow (Calc (a, c), b) % MOD;
	}
}

\(f (n) = \prod _ {T = 1} ^ n \prod _ {d \mid T} d ^ {\mu \left ( \frac T d \right)}\)，\(g (n) = n !\)

\(Type ~ 1\)

\[\begin {aligned} \prod _ {i = 1} ^ A { \prod _ {j = 1} ^ B { \prod _ {k = 1} ^ C { \left ( \dfrac {\operatorname {lcm} (i, j)} {\gcd (i, k)} \right ) ^ {ijk} } } } \end {aligned} \]

\(O (n ^ {0.75} \log n)\)

我们还是先把 \(\operatorname {lcm}\) 干掉，分子上干净一点

\[\begin {aligned} \prod _ {i = 1} ^ A { \prod _ {j = 1} ^ B { \prod _ {k = 1} ^ C { \left ( \dfrac {ij} {\gcd (i, j) \gcd (i, k)} \right ) ^ {ijk} } } } \end {aligned} \]

同样的思路，我们把 \(i\)，\(j\) 拿出来 单独处理，分子只留下 \(1\)（考虑清楚 \(i, j\) 的贡献次数）

\[\begin {aligned} \left ( \prod _ {i = 1} ^ A { i ^ { i \sum _ {j = 1} ^ B { \sum _ {k = 1} ^ C { jk } } } } \right ) \left ( \prod _ {j = 1} ^ B { j ^ { j \sum _ {i = 1} ^ A { \sum _ {k = 1} ^ C { ik } } } } \right ) \dfrac 1 { \left ( \prod _ {i = 1} ^ A { \prod _ {j = 1} ^ B { \gcd (i, j) ^ { ij \sum _ {k = 1} ^ C k } } } \right ) \left ( \prod _ {i = 1} ^ A { \prod _ {k = 1} ^ C { \gcd (i, k) ^ { ik \sum _ {j = 1} ^ B j } } } \right ) } \end {aligned} \]

前两项直接做是 \(O (n \log n)\) 的，可以过但没必要，简单变形一下就可以呈 下面的形式

\[\begin {aligned} \left ( \prod _ {i = 1} ^ A { i ^ i } \right ) ^ { \left ( \sum _ {j = 1} ^ B j \right ) \left ( \sum _ {k = 1} ^ C k \right ) } \end {aligned} \]

于是预处理 \(\prod _ {i = 1} ^ A {i ^ i}\) 就可以了，这样单次询问就是 \(O (\log n)\) 的了

然后就是 后面一大堆，显然，分母两个括号内 仍然同构，也就是下面的形式

\[\begin {aligned} \prod _ {i = 1} ^ n { \prod _ {j = 1} ^ m { \gcd (i, j) ^ {ij} } } \end {aligned} \]

据说这东西 后面没有用，于是这里只追求一个单次 $ \le O (n)$ 的形式就行了，不用 \(O (\sqrt n \log n)\)

继续枚举 \(\gcd\)，变成这样的形式，和上面 \(Type ~ 0\) 如出一辙

\[\begin {aligned} \prod _ {d = 1} ^ {\min (n, m)} { d ^ { \sum _ {i = 1} ^ { \left \lfloor \frac n d \right \rfloor } \sum _ {j = 1} ^ { \left \lfloor \frac m d \right \rfloor } [\gcd (i, j) = 1] ijd ^ 2 } } \end {aligned} \]

多了个 \(ijd ^ 2\)，还是可以套 Trick 1，只考虑 指数部分，就变成了下面的样子

\[\begin {aligned} d ^ 2 \sum _ {k = 1} ^ { \min \left ( \left \lfloor \frac n d \right \rfloor , \left \lfloor \frac m d \right \rfloor \right ) } { \mu (k) k ^ 2 \dfrac { \left \lfloor \frac n {kd} \right \rfloor \left ( \left \lfloor \frac n {kd} \right \rfloor + 1 \right ) } 2 \dfrac { \left \lfloor \frac m {kd} \right \rfloor \left ( \left \lfloor \frac m {kd} \right \rfloor + 1 \right ) } 2 } \end {aligned} \]

我们再预处理一下 \(\prod _ {d = 1} ^ {\min (n, m)} d ^ {d ^ 2}\)，然后对 \(\sum _ {k = 1}\) 这一部分 整除分块 即可

这样的话就和 \(Type ~ 0 - 1\) 一样，是 两次整除分块，询问复杂度 \(O (n ^ {0.75} \log n)\)

namespace Type_1_1 {
	int64_t f[MAXN], rf[MAXN], g[MAXN], h[MAXN];
	
	inline void Init (const int n = 100000) {
		f[0] = rf[0] = 1, h[0] = 1;
		for (int i = 1; i <= n; ++ i) h[i] = qpow (i, i);
		for (int i = 1; i <= n; ++ i) h[i] = h[i] * h[i - 1] % MOD;
		for (int i = 1; i <= n; ++ i) f[i] = qpow (i, 1ll * i * i);
		for (int i = 1; i <= n; ++ i) f[i] = f[i] * f[i - 1] % MOD;
		for (int i = 1; i <= n; ++ i) g[i] = mu[i] * i * i % (MOD - 1);
		for (int i = 1; i <= n; ++ i) g[i] = g[i] + g[i - 1], g[i] = g[i] % (MOD - 1);
		
		for (int i = 1; i <= n; ++ i) rf[i] = qpow (f[i]);
	}
	
	inline int64_t sum (const int64_t a) {
		return (a * (a + 1) / 2) % (MOD - 1);
	}
	
	inline int64_t calc (const int n, const int m) {
		int64_t ans = 0, d = 0, v = 0;
		for (int l = 1, r; l <= std::min (n, m); l = r + 1) {
			r = std::min (n / (n / l), m / (m / l)), d = n / l, v = m / l;
			ans += (g[r] - g[l - 1] + MOD - 1) % (MOD - 1) * sum (d) % (MOD - 1) * sum (v) % (MOD - 1), ans = ans % (MOD - 1);
		}
		return (ans + MOD - 1) % (MOD - 1);
	} 
	
	inline int64_t Calc (const int n, const int m) {
		int64_t ans = 1, d = 0, v = 0;
		for (int l = 1, r; l <= std::min (n, m); l = r + 1) {
			r = std::min (n / (n / l), m / (m / l)), d = n / l, v = m / l;
			ans = ans * qpow (f[r] * rf[l - 1] % MOD, calc (d, v)) % MOD;
		}
		return qpow (ans); // NOTE 1 / ans
	}
	
	inline int64_t Solve (const int64_t a, const int64_t b, const int64_t c) {
		return qpow (h[a], sum (b) * sum (c)) * qpow (h[b], sum (a) * sum (c)) % MOD * qpow (Calc (a, b), sum (c)) % MOD * qpow (Calc (a, c), sum (b)) % MOD;
	}	
}

注意这部分要预处理 三个东西，\(f (i) = \prod i ^ i\)，\(g (d) = \prod d ^ {d ^ 2}\) 和 \(h (k) = \sum = \mu (k) k ^ 2\)

\(O (\sqrt n \log n)\)

有 \(1\) 就有 \(2\)，我们同样可以优化这个式子，还是先把 \(\sum _ k\) 搞下来，变成 \(\prod _ k\)

\[\begin {aligned} \prod _ {d = 1} ^ {\min (n, m)} { \prod _ {k = 1} ^ {\min (n, m)} { d ^ { d ^ 2 \mu (k) k ^ 2 \frac { \left \lfloor \frac n {kd} \right \rfloor \left ( \left \lfloor \frac n {kd} \right \rfloor + 1 \right ) } 2 \frac { \left \lfloor \frac m {kd} \right \rfloor \left ( \left \lfloor \frac m {kd} \right \rfloor + 1 \right ) } 2 } } } \end {aligned} \]

对上界变了有疑问看上面 \(Type ~ 0 - 2\)

同样令 \(T = kd\)，那么显然 \(k ^ 2\) 与 \(d ^ 2\) 可以合并了，于是有

\[\begin {aligned} \prod _ {d = 1} ^ {\min (n, m)} { \prod _ {k = 1} ^ {\min (n, m)} { d ^ { \mu \left ( \frac T d \right ) T ^ 2 \frac { \left \lfloor \frac n T \right \rfloor \left ( \left \lfloor \frac n T \right \rfloor + 1 \right ) } 2 \frac { \left \lfloor \frac m T \right \rfloor \left ( \left \lfloor \frac m T \right \rfloor + 1 \right ) } 2 } } } \end {aligned} \]

于是 \(k\) 就又废了，转换枚举 \(T\)，\(d\) 为 \(T\) 因数，有

\[\begin {aligned} \prod _ {T = 1} ^ {\min (n, m)} { \left ( \left ( \prod _ {d \mid T} d ^ { \mu \left ( \frac T d \right ) } \right ) ^ {T ^ 2} \right ) ^ { \frac { \left \lfloor \frac n T \right \rfloor \left ( \left \lfloor \frac n T \right \rfloor + 1 \right ) } 2 \frac { \left \lfloor \frac m T \right \rfloor \left ( \left \lfloor \frac m T \right \rfloor + 1 \right ) } 2 } } \end {aligned} \]

显然 \(\prod _ {d \mid T} d ^ {\mu \left ( \frac T d \right)}\) 是和 \(Type ~ 0 - 1\) 一样的部分，容易预处理，之后乘上 \(T ^ 2\) 也是简单的

于是对 \(T\) 分段，同样的 整除分块 就可以做了，时间复杂度 \(O (\sqrt n \log n)\)

namespace Type_1_2 {
	int64_t f[MAXN], rf[MAXN], g[MAXN], rg[MAXN], h[MAXN];
	
	inline void Init (const int n = 100000) {
		h[0] = 1, g[0] = 1, f[0] = rf[0] = 1;
		for (int i = 1; i <= n; ++ i) g[i] = g[i - 1] * i % MOD;
		for (int i = 1; i <= n; ++ i) rg[i] = qpow (g[i]);
		for (int i = 1; i <= n; ++ i) h[i] = qpow (i, i);
		for (int i = 1; i <= n; ++ i) h[i] = h[i] * h[i - 1] % MOD;
		
		for (int i = 1; i <= n; ++ i) f[i] = 1;
		for (int i = 1; i <= n; ++ i) 
			for (int j = 1, k = i; k <= n; ++ j, k = i * j) {
				if (mu[j] == + 1) f[k] = f[k] * i % MOD;
				if (mu[j] == - 1) f[k] = f[k] * rg[i] % MOD * g[i - 1] % MOD;
			}
		
		for (int i = 1; i <= n; ++ i) f[i] = qpow (f[i], 1ll * i * i);
		for (int i = 2; i <= n; ++ i) f[i] = f[i] * f[i - 1] % MOD;
		for (int i = 1; i <= n; ++ i) rf[i] = qpow (f[i]);
	}
	
	inline int64_t sum (const int64_t a) {
		return (a * (a + 1) / 2) % (MOD - 1);
	}
	
	inline int64_t Calc (const int n, const int m) {
		int64_t ans = 1, d = 0, v = 0;
		for (int l = 1, r; l <= std::min (n, m); l = r + 1) {
			r = std::min (n / (n / l), m / (m / l)), d = n / l, v = m / l;
			ans = ans * qpow (f[r] * rf[l - 1] % MOD, sum (d) * sum (v) % (MOD - 1)) % MOD;
		}
		return qpow (ans); // NOTE 1 / ans
	}
	
	inline int64_t Solve (const int64_t a, const int64_t b, const int64_t c) {
		return qpow (h[a], sum (b) * sum (c)) * qpow (h[b], sum (a) * sum (c)) % MOD * qpow (Calc (a, b), sum (c)) % MOD * qpow (Calc (a, c), sum (b)) % MOD;
	}	
}

\(f (n) = \prod _ {T = 1} ^ n \prod _ {d \mid T} d ^ {\mu \left ( \frac T d \right)}\)，\(g (n) = n !\)，\(h (i) = \prod i ^ i\)

\(Type 2\)

\(O (n ^ {0.75} \log n)\)

\[\begin {aligned} \prod _ {i = 1} ^ A { \prod _ {j = 1} ^ B { \prod _ {k = 1} ^ C { \left ( \dfrac {\operatorname {lcm} (i, j)} {\gcd (i, k)} \right ) ^ {\gcd (i, j, k)} } } } \end {aligned} \]

前面把 \(\operatorname {lcm}\) 变成 \(\dfrac x {\gcd}\)，将分子变为 \(1\) 都是 老套路 了，变完之后就是下面的式子

\[\begin {aligned} \left ( \prod _ {i = 1} ^ A { i ^ { \sum _ {j = 1} ^ B { \sum _ {k = 1} ^ C { \gcd (i, j, k) } } } } \right ) \left ( \prod _ {j = 1} ^ B { j ^ { \sum _ {i = 1} ^ A { \sum _ {k = 1} ^ C { \gcd (i, j, k) } } } } \right ) \dfrac 1 { \left ( \prod _ {i = 1} ^ A { \prod _ {j = 1} ^ B { \gcd (i, j) ^ { \sum _ {k = 1} ^ C { \gcd (i, j, k) } } } } \right ) \left ( \prod _ {i = 1} ^ A { \prod _ {k = 1} ^ C { \gcd (i, k) ^ { \sum _ {j = 1} ^ C { \gcd (i, j, k) } } } } \right ) } \end {aligned} \]

这下前两项由于 指数上有 \(\gcd\)，所以也需要做一些变换了，这里以 第一项为例

我们先枚举 \(\gcd\) 的取值 \(d\)，并把它拿下来

\[\begin {aligned} \prod _ {d = 1} ^ {\min (A, B, C)} { \prod _ {i = 1} ^ A { i ^ { \sum _ {j = 1} ^ B { \sum _ {k = 1} ^ C { [\gcd (i, j, k) = d] d } } } } } \end {aligned} \]

对于 \(i, j, k\)，我们均只枚举 \(d\) 的倍数，这样可以把 \([\gcd (i, j, k) = d]\) 变成 \([\gcd (i, j, k) = 1]\)

\[\begin {aligned} \prod _ {d = 1} ^ {\min (A, B, C)} { \prod _ {i = 1} ^ { \left \lfloor \frac A d \right \rfloor } { (id) ^ { \sum _ {j = 1} ^ { \left \lfloor \frac B d \right \rfloor } { \sum _ {k = 1} ^ { \left \lfloor \frac C d \right \rfloor } { [gcd (i, j, k) = 1] d } } } } } \end {aligned} \]

\([\gcd (i, j, k)] = 1\) 出现了！于是可以 莫比乌斯反演

\[\begin {aligned} \prod _ {d = 1} ^ {\min (A, B, C)} { \prod _ {i = 1} ^ { \left \lfloor \frac A d \right \rfloor } { (id) ^ { \sum _ {j = 1} ^ { \left \lfloor \frac B d \right \rfloor } { \sum _ {k = 1} ^ { \left \lfloor \frac C d \right \rfloor } { \sum _ {t \mid \gcd (i, j, k)} \mu (t) d } } } } } \end {aligned} \]

再把 \(t\) 拿下来，放到前面

\[\begin {aligned} \prod _ {d = 1} ^ {\min (A, B, C)} { \prod _ {t = 1} ^ {\min (A, B, C)} { \prod _ {i = 1} ^ { \left \lfloor \frac A d \right \rfloor } { (id) ^ { \sum _ {j = 1} ^ { \left \lfloor \frac B d \right \rfloor } { \sum _ {k = 1} ^ { \left \lfloor \frac C d \right \rfloor } { [t \mid i \wedge t \mid j \wedge t \mid k] \mu (t) d } } } } } } \end {aligned} \]

同样的，我们只要 \(t\) 的倍数，于是

\[\begin {aligned} \prod _ {d = 1} ^ {\min (A, B, C)} { \prod _ {t = 1} ^ {\min (A, B, C)} { \prod _ {i = 1} ^ { \left \lfloor \frac A {td} \right \rfloor } { (itd) ^ { \mu (t) d \left \lfloor \frac B {td} \right \rfloor \left \lfloor \frac C {td} \right \rfloor } } } } \end {aligned} \]

然后看到 \(td\)，十分令人恶心，于是换成 \(T\)

\[\begin {aligned} \prod _ {d = 1} ^ {\min (A, B, C)} { \prod _ {t = 1} ^ {\min (A, B, C)} { \prod _ {i = 1} ^ { \left \lfloor \frac A T \right \rfloor } { (iT) ^ { \mu \left ( \frac T d \right ) d \left \lfloor \frac B T \right \rfloor \left \lfloor \frac C T \right \rfloor } } } } \end {aligned} \]

这样 \(t\) 又没用了，于是

\[\prod _ {T = 1} ^ {\min (A, B, C)} { \prod _ {d \mid T} { \prod _ {i = 1} ^ { \left \lfloor \frac A T \right \rfloor } { (iT) ^ { \mu \left ( \frac T d \right ) d \left \lfloor \frac B T \right \rfloor \left \lfloor \frac C T \right \rfloor } } } } \]

这里需要一些注意力，考虑指数上有 \(\sum _ {d \mid T} d \mu \left ( \frac T d \right)\) 的结构，于是用 \(\operatorname {id} * \mu = \varphi\)

可能上式指数上并没有这样的东西，但是把底数的 \(\prod _ {d \mid T}\) 放上去就行了

于是就成了

\[\begin {aligned} \prod _ {T = 1} ^ {\min (A, B, C)} { \prod _ {i = 1} ^ { \left \lfloor \frac A T \right \rfloor } { (iT) ^ { \varphi (T) \left \lfloor \frac B T \right \rfloor \left \lfloor \frac C T \right \rfloor } } } \end {aligned} \]

但现在这个 \(iT\) 就让人 不知所措，直接分块的话似乎 前缀积比较复杂，于是我们再拆一下

\[\begin {aligned} \left ( \prod _ {T = 1} ^ {\min (A, B, C)} { T ^ { \varphi (T) \left \lfloor \frac A T \right \rfloor \left \lfloor \frac B T \right \rfloor \left \lfloor \frac C T \right \rfloor } } \right ) \left ( \prod _ {T = 1} ^ {\min (A, B, C)} { \prod _ {i = 1} ^ { \left \lfloor \frac A T \right \rfloor } { i ^ { \varphi (T) \left \lfloor \frac B T \right \rfloor \left \lfloor \frac C T \right \rfloor } } } \right ) \end {aligned} \]

这样就很牛了，前面后面都可以 简单整除分块 完成，需要的只有 \(i !\) 和 \(\varphi\) 的前缀和了

这时候咱 先别急着想第一部分咋办，我们先来处理 分母上的东西，两个这样 同构的式子

\[\begin {aligned} \prod _ {i = 1} ^ {A} { \prod _ {j = 1} ^ B { \gcd (i, j) ^ { \sum _ {k = 1} ^ C { \gcd (i, j, k) } } } } \end {aligned} \]

我们继续枚举 \(\gcd\)，同样把 \(d\) 拿出来

\[\prod _ {d = 1} ^ {\min (A, B, C)} { \prod _ {i = 1} ^ A { \prod _ {j = 1} ^ B { \gcd (i, j) ^ { \sum _ {k = 1} ^ C { [\gcd (i, j, k) = d] d } } } } } \]

枚举 \(i, j, k\) 为 \(d\) 的倍数

\[\begin {aligned} \prod _ {d = 1} ^ {\min (A, B, C)} { \prod _ {i = 1} ^ { \left \lfloor \frac A d \right \rfloor } { \prod _ {j = 1} ^ { \left \lfloor \frac B d \right \rfloor } { (d \gcd (i, j)) ^ { \sum _ {k = 1} ^ { \left \lfloor \frac C d \right \rfloor } { [\gcd (i, j, k) = 1] d } } } } } \end {aligned} \]

对 \([\gcd (i, j, k) = 1]\) 莫比乌斯反演

\[\begin {aligned} \prod _ {d = 1} ^ {\min (A, B, C)} { \prod _ {i = 1} ^ { \left \lfloor \frac A d \right \rfloor } { \prod _ {j = 1} ^ { \left \lfloor \frac B d \right \rfloor } { (d \gcd (i, j)) ^ { \sum _ {k = 1} ^ { \left \lfloor \frac C d \right \rfloor } { \sum _ {t \mid i \wedge t \mid j \wedge t \mid k} \mu (t) d } } } } } \end {aligned} \]

枚举 \(i, j, k\) 为 \(t\) 的倍数，再将 \(\sum _ t\) 拿下来放前面

\[\begin {aligned} \prod _ {d = 1} ^ {\min (A, B, C)} { \prod _ {t = 1} ^ {\min (A, B, C)} { \prod _ {i = 1} ^ { \left \lfloor \frac A {td} \right \rfloor } { \prod _ {j = 1} ^ { \left \lfloor \frac B {td} \right \rfloor } { (td \gcd (i, j)) ^ { d \mu (t) \left \lfloor \frac C {td} \right \rfloor } } } } } \end {aligned} \]

然后经典 \(T = td\)，\(d\) 为因数，废掉 \(t\)，然后 \(\sum _ {d \mid T} d \mu \left (\frac T d \right) = \varphi (T)\)

\[\begin {aligned} \prod _ {T = 1} ^ {\min (A, B, C)} { \prod _ {i = 1} ^ { \left \lfloor \frac A T \right \rfloor } { \prod _ {j = 1} ^ { \left \lfloor \frac B T \right \rfloor } { (T \gcd (i, j)) ^ { \varphi (T) \left \lfloor \frac C T \right \rfloor } } } } \end {aligned} \]

于是又可以把 \(T\) 拿出来，就成了

\[\begin {aligned} \left ( \prod _ {T = 1} ^ {\min (A, B, C)} { T ^ { \varphi (T) \left \lfloor \frac A T \right \rfloor \left \lfloor \frac B T \right \rfloor \left \lfloor \frac C T \right \rfloor } } \right ) \left ( \prod _ {T = 1} ^ {\min (A, B, C)} { \prod _ {i = 1} ^ { \left \lfloor \frac A T \right \rfloor } { \prod _ {j = 1} ^ { \left \lfloor \frac B T \right \rfloor } { \gcd (i, j) ^ { \varphi (T) \left \lfloor \frac C T \right \rfloor } } } } \right ) \end {aligned} \]

这时候有人就发现了，欸这 第一部分 不是和上面的一样的吗？还真是

再仔细思考一下发现上面的那个是在 分子上的，而这个是在 分母上的，两者刚好 抵消了

于是这一部分我们 都 可以 不用考虑，只管 后面的部分，注意到其中

\[\begin {aligned} \prod _ {i = 1} ^ { \left \lfloor \frac A T \right \rfloor } { \prod _ {j = 1} ^ { \left \lfloor \frac B T \right \rfloor } \gcd (i, j) } \end {aligned} \]

其实就是 \(Type ~ 0\) 的式子，我们可以直接用 \(Type ~ 0\) 中的函数 \(O (\sqrt n \log n)\) 的做掉

然后对于 \(T\) 整除分块 即可，时间复杂度 \(O (n ^ {0.75} \log n)\)，于是我们就胜利了！

namespace Type_2_1 {
	int64_t f[MAXN], g[MAXN], rg[MAXN];
	
	inline void Init (const int n = 100000) {
		g[0] = rg[0] = 1; f[0] = 0;
		for (int i = 1; i <= n; ++ i) g[i] = g[i - 1] * i % MOD;
		for (int i = 1; i <= n; ++ i) rg[i] = qpow (g[i]);
	}
	
	inline int64_t Calc1 (const int a, const int b, const int c) {
		int64_t ans = 1, d = 0, v = 0;
		for (int i = 1; i <= std::max ({a, b, c}); ++ i) f[i] = (f[i - 1] + 1ll * phi[i] * (c / i) % (MOD - 1)) % (MOD - 1);
		for (int l = 1, r; l <= std::min (a, b); l = r + 1) {
			r = std::min (a / (a / l), b / (b / l)), d = a / l, v = b / l;
			ans = ans * qpow (qpow (Type_0_2::Calc (d, v)), (f[r] - f[l - 1] + MOD - 1) % (MOD - 1)) % MOD;
		}
		return qpow (ans);
	}
	
	inline int64_t Calc2 (const int a, const int b, const int c) {
		int64_t ans = 1, d = 0, v = 0, k = 0;
		for (int l = 1, r; l <= std::min ({a, b, c}); l = r + 1) {
			r = std::min ({a / (a / l), b / (b / l), c / (c / l)}), d = a / l, v = b / l, k = c / l;
			ans = ans * qpow (g[d], 1ll * (t[r] - t[l - 1] + MOD - 1) % (MOD - 1) * v % (MOD - 1) * k % (MOD - 1)) % MOD;
		}
		return ans;
	}
	
	inline int64_t Solve (const int a, const int b, const int c) {
		return Calc2 (a, b, c) * Calc2 (b, a, c) % MOD * Calc1 (a, b, c) % MOD * Calc1 (a, c, b) % MOD;
	}
	
}

\(f (n) = \prod _ {i = 1} ^ n \varphi (i) \left \lfloor \dfrac c i \right \rfloor\)，\(g (n) = n !\)

前面所有代码中 s 数组是 \(\mu\) 的前缀和，t 数组是 \(\varphi\) 的前缀和

\(O (n ^ {0.75} \log n)\)

虽然但是，我们似乎没有 更快的做法了，但是，这个 \(Type ~ 2\) 的式子确实 过于复杂

有没有什么 更加好的推狮子方法呢？看到这一大堆 \(\prod\)，考虑 \(\ln\) 转化成 \(\sum\)

\[\begin {aligned} \sum _ {i = 1} ^ A { \sum _ {j = 1} ^ B { \sum _ {k = 1} ^ C { \gcd (i, j, k) \ln { \left ( \dfrac { \operatorname {lcm} (i, j) } { \gcd (i, k) } \right ) } } } } \end {aligned} \]

枚举 \(\gcd\)

\[\begin {aligned} \sum _ {d = 1} ^ {\min (A, B, C)} { \sum _ {i = 1} ^ { \left \lfloor \frac A d \right \rfloor } { \sum _ {j = 1} ^ { \left \lfloor \frac B d \right \rfloor } { \sum _ {k = 1} ^ { \left \lfloor \frac C d \right \rfloor } { [\gcd (i, j, k) = 1] d \ln { \left ( \dfrac { \operatorname {lcm} (id, jd) } { \gcd (id, kd) } \right ) } } } } } \end {aligned} \]

莫比乌斯反演，然后进行一套 \(d \to td \to T\)，以及 \(\mu * \operatorname {id} = \varphi\)，最终得到

\[\begin {aligned} \sum _ {T = 1} ^ {\min (A, B, C)} { \varphi (d) \sum _ {i = 1} ^ { \left \lfloor \frac A T \right \rfloor } { \sum _ {j = 1} ^ { \left \lfloor \frac B T \right \rfloor } { \sum _ {k = 1} ^ { \left \lfloor \frac C T \right \rfloor } { \ln \left ( \dfrac { \operatorname {lcm} (i, j) } { \gcd (i, k) } \right ) } } } } \end {aligned} \]

最后我们只需要 \(\exp\) 一下！

\[\begin {aligned} \prod _ {T = 1} ^ {\min (A, B, C)} { \left ( \prod _ {i = 1} ^ { \left \lfloor \frac A T \right \rfloor } { \prod _ {j = 1} ^ { \left \lfloor \frac B T \right \rfloor } { \prod _ {k = 1} ^ { \left \lfloor \frac C T \right \rfloor } { \dfrac { \operatorname {lcm} (i, j) } { \gcd (i, k) } } } } \right ) ^ {\varphi (T)} } \end {aligned} \]

然后发现中间就是 \(Type ~ 0\) 啊，显然 \(\varphi\) 的 前缀和 是容易处理的，整除分块，我们就胜利了

namespace Type_2_2 {
	
	inline int64_t Solve (const int a, const int b, const int c) {
		int64_t ans = 1, d = 0, v = 0, k = 0;
		for (int l = 1, r; l <= std::min ({a, b, c}); l = r + 1) {
			r = std::min ({a / (a / l), b / (b / l), c / (c / l)}), d = a / l, v = b / l, k = c / l;
			ans = ans * qpow (Type_0_2::Solve (d, v, k), (t[r] - t[l - 1] + MOD - 1) % (MOD - 1)) % MOD;
		}
		return ans;
	}
}

其实按照 前面纯纯莫比乌斯反演 的那个做法，我们把指数上的 \(\left \lfloor \dfrac x T \right \rfloor\) 都换成 \(\prod _ {i = 1} ^ {\left \lfloor \frac x T \right \rfloor}\)

然后你就会惊奇的发现 分子部分 与 分母部分 都长成

\[\begin {aligned} \prod _ {T = 1} ^ {\min (A, B, C)} { \prod _ {i = 1} ^ { \left \lfloor \frac A T \right \rfloor } { \prod _ {j = 1} ^ { \left \lfloor \frac B T \right \rfloor } { \prod _ {k = 1} ^ { \left \lfloor \frac C T \right \rfloor } { f (x) ^ {\varphi (T)} } } } } \end {aligned} \]

这样的形式，于是我们观察 \(f (x)\)，分子上是 \(ij\)，分母是 \(\gcd (i, j)\gcd (i, k)\)

于是最后就可以变回 \(\dfrac {\operatorname {lcm} (i, j)} {\gcd (i, k)}\) 这个东西，就和 \(\ln \exp\) 推出来的结果 一模一样 了

出题人只需要考虑怎么写式子就可以了，而选手要考虑的东西就很多了

Luogu P3700 [CQOI2017] 小 Q 的表格

题意比较复杂，建议看 原题面，这个题似乎需要神秘的 注意力

我们先根据 \(2\) 式，有

\[\begin {aligned} \dfrac {f (a, b)} b = \dfrac {f (a, a + b)} {a + b} \end {aligned} \]

你 发现 左边这东西 不大对称，先给它补个 \(a\)

\[\begin {aligned} \dfrac {f (a, b)} {ab} = \dfrac {f (a, a + b)} {a (a + b)} \end {aligned} \]

注意到 右边这个 \(a + b\) 的 \(a\) 就是 凭空加上去的，于是加几个都一样

\[\begin {aligned} \dfrac {f (a, b)} {ab} = \dfrac {f (a, b \bmod a)} {a (b \bmod a)} \end {aligned} \]

考虑 \(a, b \bmod a\) 的性质，容易知道 其 \(\gcd\) 应当相等，于是

\[\begin {aligned} \dfrac {f (a, b)} {ab} = \dfrac { f (\gcd (a, b), \gcd (a, b)) } { \gcd (a, b) ^ 2 } \end {aligned} \]

于是最终得到的 可能更加有效 一些的式子就是

\[\begin {aligned} f (a, b) = f (d, d) \dfrac {ab} {d ^ 2} \end {aligned} \]

这个式子告诉我们一个什么性质？

我们修改 \(f (a, b)\) 可以对应到 \(f (d, d)\) 的修改上，\(d = \gcd (a, b)\)

同样，查询 \(f (a, b)\) 的值也可以通过查询 \(f (d, d)\) 的值得到，这样我们就把范围 缩小了

条件式 好像没办法转化了，这时候来考虑 答案式

\[\begin {aligned} \sum _ {i = 1} ^ k { \sum _ {j = 1} ^ k { f (i, j) } } \end {aligned} \]

套一下 条件式，得到

\[\begin {aligned} \sum _ {i = 1} ^ k { \sum _ {j = 1} ^ k { f (d, d) \dfrac {ij} {d ^ 2} } } \end {aligned} \]

显然 枚举一手 \(d\)，为了让 \(\dfrac {ij} {d ^ 2}\) 是个整数，我们保证 \(\gcd (i, j) = d\) 即可

\[\begin {aligned} \sum _ {d = 1} ^ k { f (d, d) \sum _{i = 1} ^ k { \sum _ {j = 1} ^ k { [\gcd (i, j) = d] \dfrac {ij} {d ^ 2} } } } \end {aligned} \]

然后经典的枚举 \(d\) 的倍数

\[\begin {aligned} \sum _ {d = 1} ^ k { f (d, d) \sum _ {i = 1} ^ { \left \lfloor \frac k d \right \rfloor } { \sum _ {j = 1} ^ { \left \lfloor \frac k d \right \rfloor } { [\gcd (i, j) = 1] ij } } } \end {aligned} \]

显然，上式分成了 \(\sum _ {d = 1} ^ k f (d, d)\) 与 \(\sum _ {i = 1} ...\) 两部分

因为修改在前面部分，故我们考虑 数据结构维护，现在只需要一个 快速求后面部分 的方法

也就是要转化下面这个式子

\[\begin {aligned} \sum _ {i = 1} ^ n { \sum _ {j = 1} ^ n { [\gcd (i, j) = 1] ij } } \end {aligned} \]

直接使用 Trick 3 - ex 即可得到 \(\sum _ {i = 1} ^ n i ^ 2 \varphi (i)\)，这东西 \(O (n)\) 预处理 后即可 \(O (1)\) 查询

如果我们对 \(d\) 整除分块，显然，我们需要维护 \(\sum _ {d = 1} ^ k f(d, d)\) 的 前缀和，同时支持 单点修改

考虑询问次数 \(O (T)\)，单次询问中需要求 \(O (\sqrt n)\) 次前缀和，而修改次数只有 \(O (T)\)

可以想到 \(O (\sqrt n) - O (1)\) 分块来做，但是写起来会复杂一些

所以笔者这里直接 \(O (\log n) - O (\log n)\) 树状数组，也是可以过的，总复杂度 \(O (T \sqrt n \log n)\)

注意 long long 的取

#include <bits/stdc++.h>

const int MAXN = 4000005;
const int MOD  = 1000000007;

int64_t f[MAXN];
int pri[MAXN], phi[MAXN], mu[MAXN], cnt = 0;
bool vis[MAXN];

inline void Sieve (const int n = 4000000) {
	mu[1] = phi[1] = 1;
	for (int i = 2; i <= n; ++ i) {
		if (!vis[i]) pri[++ cnt] = i, phi[i] = i - 1, mu[i] = - 1;
		for (int k = 1; k <= cnt && pri[k] * i <= n; ++ k) {
			vis[pri[k] * i] = 1;
			if (i % pri[k] == 0) {
				phi[i * pri[k]] = phi[i] * pri[k];
				break ;
			}
			mu[i * pri[k]] = - mu[i];
			phi[i * pri[k]] = (pri[k] - 1) * phi[i];
		}
	}
	for (int i = 1; i <= n; ++ i) f[i] = 1ll * i * i % MOD * phi[i];
	for (int i = 1; i <= n; ++ i) f[i] = (f[i - 1] + f[i]) % MOD;
}

int Q, N;

namespace BIT {
	int T[MAXN];
	
	#define lowbit(x) (x & -x)
	
	inline void Add (const int x, const int v) {
		for (int i = x; i <= N; i += lowbit (i)) T[i] = (T[i] + v) % MOD;
	}
	
	inline void Mns (const int x, const int v) {
		for (int i = x; i <= N; i += lowbit (i)) T[i] = (T[i] - v + MOD) % MOD;
	}
	
	inline int  Sum (const int x) {
		int ans = 0;
		for (int i = x; i; i -= lowbit (i)) ans = (ans + T[i]) % MOD;
		return ans;
	}
}

int64_t val[MAXN];

inline void Update (const int x, const int y, const int64_t v) {
	int d = std::__gcd (x, y);
	BIT::Mns (d, val[d] % MOD), val[d] = v / (1ll * (x / d) * (y / d)), BIT::Add (d, val[d] % MOD);
}

inline int Solve (const int k) {
	int ans = 0, d = 0;
	for (int l = 1, r; l <= k; l = r + 1) {
		r = k / (k / l), d = k / l;
		ans = (ans + 1ll * (BIT::Sum (r) - BIT::Sum (l - 1) + MOD) % MOD * f[d] % MOD) % MOD;
	}
	return ans;
}

int a, b, k;
int64_t c;

int main () {
	
	std::cin >> Q >> N, Sieve ();
	
	for (int i = 1; i <= N; ++ i) Update (i, i, 1ll * i * i);
	
	for (int i = 1; i <= Q; ++ i) {
		std::cin >> a >> b >> c >> k, Update (a, b, c);
		std::cout << Solve (k) << '\n';
	}
	
	return 0;
}

5. 引用资料

[0] Number Theory —— H_W_Y