「笔记」莫比乌斯反演

目录

• 前置知识
  • 小碎骨
  • 引理1
  • 数论分块
• 积性函数
  • 定义
  • 性质
  • 常见积性函数
• 莫比乌斯函数
  • 定义
  • 性质
  • 反演常用结论
  • 线性筛求莫比乌斯函数
• 狄利克雷(Dirichlet)卷积
  • 性质
  • 关于单位元 \(e\)
  • 除数函数与幂函数
  • 欧拉函数与恒等函数
  • 计算
• 莫比乌斯反演
  • 公式
  • 证明
  • 举例
• 关于一类求和式
  • 一般套路
  • 例 1
  • 例 2
  • 感悟
• 例题
  • [HAOI2011] Problem b
  • [国家集训队]Crash的数字表格
  • SP5971 LCMSUM - LCM Sum
  • [SDOI2015]约数个数和
  • P3768 简单的数学题
  • 「SDOI2017」数字表格
  • P5518 [MtOI2019]幽灵乐团 / 莫比乌斯反演基础练习题
• 写在最后

Updated on 2020.8.6
巨幅更新，对积性函数和狄利克雷卷积部分进行重构。
新增对一类特殊求和式的讲解。

Updated on 2020.9.9
添加了几个例题。

---

前置知识

小碎骨

艾佛森括号 \([P] = \begin{cases} 1 &\text{If P is true}\\ 0 &\text{Otherwise} \end{cases}\)
此处 \(P\) 是一个可真可假的命题。

引理1

\[\forall a,b,c\in \mathbb{Z},\left\lfloor\dfrac{a}{bc}\right\rfloor = \left\lfloor{\dfrac{\left\lfloor\dfrac{a}{b}\right\rfloor}{c}}\right\rfloor \]

证明

\[\dfrac{a}{b} = \left\lfloor{\dfrac{a}{b}}\right\rfloor + r(0\le r < 1) \]

\[\left\lfloor\dfrac{a}{bc}\right\rfloor = \left\lfloor\dfrac{a}{b}\times\dfrac{1}{c}\right\rfloor = \left\lfloor{\dfrac{1}{c}\times\left({\left\lfloor{\dfrac{a}{b} + r}\right\rfloor}\right)}\right\rfloor = \left\lfloor{\dfrac{\left\lfloor\dfrac{a}{b}\right\rfloor}{c} + \dfrac{r}{c}}\right\rfloor = \left\lfloor{\dfrac{\left\lfloor{\dfrac{a}{b}}\right\rfloor}{c}}\right\rfloor \]

数论分块

内容独立了出来，详细内容见 数论分块 - Luckyblock

对于一类含有\(\left\lfloor\frac{n}{i}\right\rfloor\)的求和式 (\(n\) 为常数)，由于\(\left\lfloor\frac{n}{i}\right\rfloor\)单调不增，故存在多个区间\([l,r]\), 使得\(\left\lfloor\frac{n}{i}\right\rfloor = \left\lfloor\frac{n}{j}\right\rfloor(i,j\in [l,r])\) 成立。

对于任意一个\(i\)，最大的满足上式的 \(j=\left\lfloor{\dfrac{n}{\left\lfloor{\dfrac{n}{i}}\right\rfloor}}\right\rfloor\)

---

积性函数

定义

若 \(\gcd(x,y) = 1\) 且 \(f(xy)=f(x)f(y)\)，则 \(f(n)\) 为积性函数。

性质

若 \(f(x)\)，\(g(x)\)均为积性函数，则以下函数也为积性函数：

\[\begin{aligned} & h(x) = f(x^p)\\ & h(x) = f^p(x)\\ & h(x) = f(x)g(x)\\ & h(x) = \sum_{d\mid x} f(d)g\left(\dfrac{x}{d}\right) \end{aligned}\]

常见积性函数

• 单位函数 \(e(n) = [n = 1]\)。
• 幂函数 \(\operatorname{Id}_{k}(n) = n^k\)， \(\operatorname{id}_1(n)\) 通常简记为\(\operatorname{id}(n)\)。
• 常数函数 \(1(n) = 1\)。
• 因数个数 \(\operatorname{d}(n) = \sum\limits_{d\mid n} 1\)。
• 除数函数 \(\sigma_{k}(n) = \sum\limits_{d\mid n} d^k\)。
  \(k=1\) 时为因数和函数，通常简记为 \(\sigma(n)\)。
  \(k=0\) 时为因数个数函数 \(\sigma_{0}(n)\)。
• 欧拉函数 \(\varphi(n) = \sum\limits_{i=1}^{n} [\gcd(i,n) = 1]\)。
• 莫比乌斯函数 \(\mu(n) = \begin{cases}1 &n=1\\0 &n\ \text{含有平方因子}\\(-1)^k &k\text{为}\ n\ \text{的本质不同质因子个数} \end{cases}\)

不是很懂上面写的什么玩意？
不用深究，有个印象继续往下看就好。

---

莫比乌斯函数

定义

\(\mu\) 为莫比乌斯函数，定义为

\[\mu(n) = \begin{cases} 1 &n=1\\0 &n\ \text{含有平方因子}\\(-1)^k &k\text{为}\ n\ \text{的本质不同质因子个数} \end{cases}\]

解释

令 \(n = \prod\limits_{i=1}^{k} p_{i}^{c_i}\)，其中\(p_i\)为质因子，\(c_i\ge 1\)。

1. \(n=1\)时，\(\mu (n) = 1\)。
2. \(n\not ={1}\)时 ，
   • \(\exist i\in [1,k], c_i > 1\) 时，\(\mu (n) = 0\)。
     当某质因子出现次数大于\(1\)时，\(\mu (n) = 0\)。
   • \(\forall i\in [1,k], c_i = 1\) 时，\(\mu (n) = (-1)^k\)。
     当每个质因子只出现一次时，即\(n = \prod\limits_{i=1}^{k}p_i\)，\(\{p_i\}\)中元素唯一。
     \(\mu (n) = (-1)^k\)，此处\(k\)为质因子的种类数。

性质

莫比乌斯函数是积性函数，且具有以下性质

\[\large \sum_{d\mid n} \mu (d) = [n=1] \]

证明，设 \(n = \prod\limits_{i=1}^{k}{p_i^{c_i}}, n' = \prod\limits_{i=1}^{k}{p_i}\)。

• 根据莫比乌斯函数定义，则有：\(\sum\limits_{d\mid n}{\mu(d)} = \sum\limits_{d\mid n'}{\mu(d)}\)。
• 若 \(n'\) 的某因子 \(d\)，有 \(\mu (d) = (-1)^i\)，则它由 \(i\) 个 本质不同的质因子组成。
  由于质因子总数为 \(k\)，满足上式的因子数为 \(C_{k}^{i}\)。
• 对于原求和式，转为枚举 \(\mu(d)\) 的值。
  则 \(\sum\limits_{d\mid n'}{\mu(d)} = \sum\limits_{i=0}^{k}{C_{k}^{i} \times (-1)^i} = \sum\limits_{i=0}^{k}{C_{k}^{i} \times (-1)^i\times 1^{k-i}}\)
  根据二项式定理，上式 \(= (1+(-1))^k\)，
  易知该式在 \(k=0\)，即 \(n=0\) 时为 \(1\)，否则为 \(0\)。

反演常用结论

一个反演常用结论:

\[[\gcd(i,j) = 1] = \sum\limits_{d\mid \gcd(i,j)} {\mu (d)} \]

证明 1：
设 \(n = \gcd(i,j)\)，则右\(= \sum\limits_{d\mid n} {\mu (d)} = [n = 1] = [\gcd(i,j) = 1]=\) 左。

证明 2：
暴力展开：\([\gcd(i,j) = 1] = e(\gcd(i,j)) = \sum\limits_{d\mid \gcd(i,j)}\mu(d)\)

线性筛求莫比乌斯函数

\(\mu\) 为积性函数，因此可以线性筛莫比乌斯函数。

int pnum, mu[kMaxn], p[kMaxn];
bool vis[kMaxn];

void Euler(int lim_) {
  vis[1] = true, mu[1] = 1ll;
  for (int i = 2; i <= lim_; ++ i) {
    if (! vis[i]) {
      p[++ pnum] = i;
      mu[i] = - 1;
    }
    for (int j = 1; j <= pnum && i * p[j] <= lim_; ++ j) {
      vis[i * p[j]] = true;
      if (i % p[j] == 0) { //勿忘平方因子
        mu[i * p[j]] = 0;
        break;
      }
      mu[i * p[j]] = - mu[i];
    }
  }
}

---

狄利克雷(Dirichlet)卷积

建议阅读 算法学习笔记(35): 狄利克雷卷积 By: Pecco。

定义两个数论函数 \(f,g\) 的狄利克雷卷积为

\[\large(f\ast g) (n) = \sum_{d\mid n} f(d)g\left(\dfrac{n}{d}\right) \]

性质

1. 显然满足 交换律，结合律，分配律。
   • \(f \ast g = g \ast f\)
   • \((f \ast g) \ast h = f\ast (g\ast h)\)
   • \(f\ast (g+h) = f\ast g + f\ast h\)
2. \(e\) 为狄利克雷卷积的单位元，有\((f\ast e)(n) = f(n)\)。
3. 若 \(f,g\) 为积性函数，则 \(f\ast g\) 为积性函数。

关于单位元 \(e\)

有：

\[e = \mu \ast 1=\sum\limits_{d\mid n} \mu (d) \]

证明

\[\begin{aligned} (f\ast e)(n) = \sum_{d\mid n} f(d)e(\dfrac{n}{d}) = \sum_{d\mid n} f(d)\left[\dfrac{n}{d} = 1\right] \end{aligned}\]

• 对于\([\dfrac{n}{d} = 1]\)，当且仅当 \(\dfrac{n}{d}=1\)，即 \(d=n\) 时为 \(1\)，否则为\(0\)。
• 则当 \(d=n\) 时，\(f(d)\left[\dfrac{n}{d} = 1\right] = f(n)\)。
  当 \(d\not ={n}\) 时，\(f(d)\left[\dfrac{n}{d} = 1\right] = 0\)。

综上，\((f\ast e)(n) = \sum\limits_{d\mid n} f(d)\left[\dfrac{n}{d} = 1\right] = f(n)\)，满足单位元的性质。
则 \(e = \mu \ast 1\) 成立。

除数函数与幂函数

幂函数 \(\operatorname{Id}_{k}(n) = n^k\)。
除数函数 \(\sigma_{k}(n) = \sum\limits_{d\mid n} d^k\)。

显然有：

\[(\operatorname{Id}_k\ast 1)(n) = \sum_{d\mid n} \operatorname{Id_k}(d) = \sum_{d\mid n} d^k = \sigma_k(n) \]

当 \(k=0\) 时，\(\operatorname{Id_0}=1\)，\(\sigma_0\) 为因数个数函数，有：

\[(1\ast 1)(n) = \sum_{d\mid n}1 = \sigma_0(n) \]

欧拉函数与恒等函数

\[\begin{aligned} \varphi \ast 1 =& \operatorname{Id}\\ \varphi =& \operatorname{Id}\ast \mu \end{aligned}\]

对于一式，当 \(n=p^m\) 时（\(p\) 为质数），有：

\[(\varphi \ast 1)(p^m) = \sum_{d\mid n}\varphi(d) = \varphi(1) +\sum_{i=1}^{m}\varphi(p^i) = 1 +\sum_{i=1}^{m}(p^i-p^{i-1}) = p^m \]

\(p^i\) 的因子有 \(p^{i-1}\) 个，为 \(1\sim p^{i-1}\)，故 \(\varphi(p^i) = p^i-p^{i-1}\)。

又 \((\varphi \ast 1)(n)\) 为积性函数，则对于任意正整数 \(n\)，有：

\[(\varphi \ast 1)(n) = (\varphi \ast 1)\left(\prod p^m\right) = \prod\left(\varphi \ast 1\right)(p^m) = \prod p^m = n \]

得证。

对于 2 式，在 1 式基础上两侧同时 \(\ast \mu\) 即得。
左侧变为 \(\varphi \ast 1\ast \mu = \varphi \ast e = \varphi\)。

计算

\[(f\ast g) (n) = \sum_{d\mid n} f(d)g\left(\dfrac{n}{d}\right) \]

考虑枚举 \(n\) 的因子，将贡献累加即可。
显然可以使用埃氏筛筛出所有前缀狄利克雷卷积，复杂度 \(O(nk\log n)\)，其中 \(k\) 是计算函数值的复杂度。

---

莫比乌斯反演

反演是个啥？反演。

公式

设\(f(n),g(n)\)为两个数论函数。
如果有

\[\large f(n) = (g\ast 1)(n) = \sum\limits_{d\mid n}{g(d)} \]

那么有

\[\large g(n) = (\mu \ast f)(n)=\sum\limits_{d\mid n} {\mu(d)f\left(\dfrac{n}{d}\right)} \]

证明

方法一：运用卷积。

原问题为：已知 \(f = g\ast 1\)，证明 \(g = f\ast \mu\)。
易知如下转化：\(f\ast \mu = g\ast 1 \ast \mu \Longrightarrow f\ast \mu = g\ast e = g\)。

方法二：对原式进行数论变换。

1. 用 \(\sum\limits_{d\mid n}g(d)\) 替换\(f\left(\dfrac{n}{d}\right)\)。

   \[\sum_{d\mid n}{\mu(d)\sum_{k\mid \frac{n}{d}}g(k)} \]

2. 变换求和顺序。

   \[\sum_{k\mid n}g(k)\sum_{d\mid \frac{n}{k}}{\mu(d)} \]

3. 依据 \(\sum\limits_{d\mid n}{\mu(d)} = [n=1]\)，仅当 \(\dfrac{n}{k} = 1\) 时，\(\sum\limits_{d\mid \frac{n}{k}}{\mu(d)} = 1\)，否则为 \(0\)。
   此时\(k=n\)，故上式等价于 \(\sum\limits_{k\mid n} {[n=k]\cdot g(k)} = g(n)\)。

举例

从 狄利克雷(Dirichlet)卷积 部分可以知道一些积性函数的关系。
尝试对它们进行反演：

\[e = \mu \ast 1\iff\mu = \mu\ast e = \mu \]

\[\sigma_k = \operatorname{Id}_k\ast 1\iff \operatorname{Id}_k = \mu \ast \sigma_k \]

\[\operatorname{Id}=\varphi \ast 1\iff \varphi = \operatorname{Id}\ast \mu \]

---

关于一类求和式

\[\sum_{i=1}^{n}\sum_{j=1}^{m}f(\gcd(i,j)) \]

一般套路

考虑构造出函数 \(g\)，满足下列形式：

\[f(n) = g\ast 1 = \sum_{d\mid n}g(d) \]

则求和式变为：

\[\sum_{i=1}^{n}\sum_{j=1}^{m}\sum_{d\mid \gcd(i,j)}g(d) \]

考虑算术基本定理，发现若满足 \(d\mid \gcd (i,j)\)，则 \(d\mid i\) 且 \(d\mid j\) 成立。
考虑 \(g(d)\) 在何时做出贡献，调整枚举顺序：

\[\sum_{d=1}g(d)\sum_{i=1}^n[d\mid i]\sum_{j=1}^m[d\mid j] \]

\(\sum\limits_{i=1}^{n}[d\mid i]\) 等价于 \(1\sim n\) 中 \(d\) 的倍数的个数，则上式等价于：

\[\sum_{d=1}g(d)\left\lfloor\dfrac{n}{d}\right\rfloor\left\lfloor\dfrac{m}{d}\right\rfloor \]

数论分块求解即可。

例 1

\[\sum_{i=1}^{n}\sum_{j=1}^{m}\gcd(i,j) \]

发现此题的 \(f\) 等价于 \(\operatorname{Id}\)，则上式等价于：

\[\begin{aligned} &\sum_{i=1}^{n}\sum_{j=1}^{m}\operatorname{Id}(\gcd(i,j))\\ =& \sum_{i=1}^{n}\sum_{j=1}^{m}\sum_{d\mid \gcd(i,j)}\varphi(d)\\ =& \sum_{d=1}\varphi(d)\left\lfloor\dfrac{n}{d}\right\rfloor\left\lfloor\dfrac{m}{d}\right\rfloor \end{aligned}\]

例 2

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

感悟

卷点什么东西，把 \(g\) 卷出来。
\(g\) 不一定是特殊意义的函数。

---

例题

[HAOI2011] Problem b

\(n\) 组询问，每次给定参数 \(a,b,c,d,k\)，求：

\[\sum\limits_{x=a}^{b}\sum\limits_{y=c}^{d}[\gcd(x,y) = k] \]

\(1\le n,k,a,b,c,d\le 5\times 10^4\)，\(a\le b,c\le d\)。

令 \(f(n,m) = \sum\limits_{i=1}^{n}\sum\limits_{j=1}^{m}[\gcd(i,j) = k]\)。
根据容斥原理，则原式等价于 \(f(b,d) - f(a-1,d) - f(b,d-1) + f(a-1,d-1)\)。
\(f\) 变成了上述一类求和式的形式，考虑化简 \(f\)。

易知原式等价于

\[\sum\limits_{i=1}^{\left\lfloor\frac{n}{k}\right\rfloor}\sum\limits_{j=1}^{\left\lfloor\frac{m}{k}\right\rfloor}[\gcd(i,j) = 1] \]

代入反演常用结论 \([\gcd(i,j) = 1] = \sum\limits_{d\mid \gcd(i,j)} {\mu (d)}\)，上式化为

\[\sum\limits_{i=1}^{\left\lfloor\frac{n}{k}\right\rfloor}\sum\limits_{j=1}^{\left\lfloor\frac{m}{k}\right\rfloor}\sum_{d\mid \gcd(i,j)}{\mu(d)} \]

变换求和顺序，先枚举\(d\mid \gcd(i,j)\)，可得

\[\sum_{d=1}\mu(d)\sum_{i=1}^{\left\lfloor\frac{n}{k}\right\rfloor}[d\mid i]\sum_{j=1}^{\left\lfloor\frac{m}{k}\right\rfloor}[d\mid j] \]

对于上式后半的解释：当\(d\mid i\) 且 \(d\mid j\) 时，\(d\mid \gcd(i,j)\)。

易知\(1\sim \left\lfloor\dfrac{n}{k}\right\rfloor\) 中 \(d\) 的倍数有 \(\left\lfloor\dfrac{\left\lfloor\dfrac{n}{k}\right\rfloor}{d}\right\rfloor = \left\lfloor\dfrac{n}{kd}\right\rfloor\) 个（由引理 1 可知），原式变为

\[\sum_{d=1}\mu(d)\left\lfloor\dfrac{n}{kd}\right\rfloor\left\lfloor\dfrac{m}{kd}\right\rfloor \]

预处理 \(\mu\) 后，显然可以数论分块求解，复杂度\(O(n + T\sqrt{n})\)。

---

代码

//知识点：莫比乌斯反演
/*
//By:Luckyblock
*/
#include <cstdio>
#include <cctype>
#include <algorithm>
#define ll long long
const int MARX = 6e4 + 10;
//=============================================================
int N, a, b, c, d, k;
int cnt, Mobius[MARX], Prime[MARX], Sum[MARX];
bool vis[MARX];
//=============================================================
inline int read()
{
    int f = 1, w = 0; char ch = getchar();
    for(; !isdigit(ch); ch = getchar()) if(ch == '-') f = -1;
    for(; isdigit(ch); ch = getchar()) w = (w << 3) + (w << 1) + (ch ^ '0');
    return f * w;
}
void GetMobius()
{
    Mobius[1] = 1;
    int MAX = MARX - 10;
    for(int i = 2; i <= MAX; i ++)
    {
      if(! vis[i]) Mobius[i] = - 1, Prime[++ cnt] = i;
      for(int j = 1; j <= cnt && i * Prime[j] <= MAX; j ++)
      {
        vis[i * Prime[j]] = true;
        if(i % Prime[j] == 0) break;
        Mobius[i * Prime[j]] = - Mobius[i];
      }
    }
    for(int i = 1; i <= MAX; i ++) Sum[i] = Sum[i - 1] + Mobius[i];
}
ll Calc(int x, int y)
{
    ll ans = 0ll; int max = std :: min(x, y);
    for(int l = 1, r; l <= max; l = r + 1)
      r = std :: min(x / (x / l), y / (y / l)),
      ans += (1ll * x / (1ll * l * k)) * (1ll * y / (1ll * l * k)) * (Sum[r] - Sum[l - 1]);
    return ans;
}
ll Solve()
{
    a = read(), b = read(), c = read(), d = read(), k = read();
    return Calc(b, d) - Calc(b, c - 1) - Calc(a - 1, d) + Calc(a - 1, c - 1);
}
//=============================================================
int main()
{ 
    N = read(); GetMobius();
    while(N --) printf("%lld\n", Solve());
    return 0;
}

---

[国家集训队]Crash的数字表格

给定 \(n,m\) , 求：

\[\sum_{i=1}^n\sum_{j=1}^{m} \operatorname{lcm}(i,j)\bmod 20101009 \]

\(1\le n,m\le 10^7\)。

易知原式等价于：

\[\sum_{i=1}^{n}\sum_{j=1}^{m}\dfrac{ij}{\gcd (i,j)} \]

考虑枚举 \(\gcd(i,j)\)，设枚举量为 \(d\)。
则 \(d=\gcd(i,j)\) 的充要条件是满足 \(d|i, d|j\) 且 \(\gcd(\dfrac{i}{d},\dfrac{j}{d}) = 1\)，则原式等价于：

\[\sum_{i=1}^n\sum_{j=1}^m \sum_{d=1} \dfrac{ij}{d}[d|i][d|j][\gcd(\frac{i}{d}, \frac{j}{d})=1] \]

先枚举 \(d\)，则原式等价于：

\[\sum_{d=1}\sum_{i=1}^{n}[d\mid i]\sum_{j=1}^m [d\mid j][\gcd(\dfrac{i}{d}, \dfrac{j}{d}=1)] \dfrac{ij}{d} \]

这个 \(d\) 很烦人，把 \(i,j\) 中的 \(d\) 提出来，变为枚举 \(\frac{i}{d}\)，\(\frac{j}{d}\)。
消去 \(d\mid i\)，\(d\mid j\) 的限定条件，则原式等价于：

\[\begin{aligned} &\sum_{d=1}\sum_{i=1}^{\left\lfloor\frac{n}{d}\right\rfloor}\sum_{j=1}^{\left\lfloor\frac{m}{d}\right\rfloor}[\gcd(i,j)=1]\dfrac{id\times jd}{d}\\ =& \sum_{d=1}\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\\ =& \sum_{d=1}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]ij \end{aligned}\]

单独考虑后半部分，设 \(f(x,y) = \sum\limits_{i=1}^{x} \sum\limits_{j=1}^{y}[\gcd(i,j)=1]ij\)。
发现 \(f(x,y)\) 的形式与 [HAOI2011] Problem b 中的式子类似，代入 \([\gcd(i,j) = 1] = \sum\limits_{d\mid \gcd(i,j)} {\mu (d)}\) 套路一波：

\[\begin{aligned} f(x,y) =& \sum\limits_{i=1}^{x} \sum\limits_{j=1}^{y}[\gcd(i,j)=1]ij\\ =& \sum_{i=1}^{x} \sum_{j=1}^{y}\sum_{d\mid \gcd(i,j)}\mu(d) ij\\ =& \sum_{d=1}\mu(d)\sum_{i=1}^{x}[d\mid i] \sum_{j=1}^{y}[d\mid j] ij\\ =& \sum_{d=1}\mu(d) d^2\sum_{i=1}^{\left\lfloor \frac{x}{d}\right\rfloor}\sum_{j=1}^{\left\lfloor\frac{y}{d}\right\rfloor}ij \end{aligned}\]

前半部分 \(\sum\limits_{d=1}\mu(d) d^2\)，可以考虑筛出 \(\mu(d)\) 后求前缀和。
后半部分是等差数列乘等差数列的形式，设 \(g(p,q) = \sum\limits_{i=1}^{p} \sum\limits_{j=1}^{q}ij\)，\(g_{p,q}\) 可以通过下式 \(O(1)\) 计算：

\[g(p,q) = \sum_{i=1}^{p} i \sum_{j=1}^{q}j= \dfrac{(1 + p)\times p}{2}\times \dfrac{(1+q)\times q}{2} \]

则对于 \(f(x,y)\)，有：

\[f(x,y) = \sum_{d=1}\mu(d) d^2\cdot g(\left\lfloor \frac{x}{d}\right\rfloor, \left\lfloor\frac{y}{d}\right\rfloor) \]

数论分块求解即可。

再看回原式，原式等价于：

\[\sum_{d=1}d\cdot f(\left\lfloor\frac{n}{d}\right\rfloor, \left\lfloor\frac{m}{d}\right\rfloor) \]

又是一个可以数论分块求解的形式。
线性筛预处理后 数论分块套数论分块，复杂度 \(O(n + m)\)，瓶颈是线性筛。

---

一些注意的点

处理时会出现求平方的运算，需要特别注意取模问题，ll 都会爆，被坑惨了。

在预处理前缀和的这个地方：

sum[i] = (sum[i - 1] + 1ll * i * i % kMod * (mu[i] + kMod)) % kMod; //仅令 mu + kMod

注意先对平方取模，在乘上 \(\mu\)，否则就会爆掉。
以及可以仅令 \(mu + mod\)。

以及这个地方：

int g(int n_, int m_) {
  return (1ll * n_ * (n_ + 1ll) / 2ll % kMod) * (1ll * m_ * (m_ + 1ll) / 2ll % kMod) % kMod;
}

平方计算，注意随时取模。

代码

//知识点：莫比乌斯反演
/*
By:Luckyblock
*/
#include <algorithm>
#include <cctype>
#include <cmath>
#include <cstdio>
#include <cstring>
#define ll long long
const ll kMod = 20101009;
const int kMaxn = 1e7 + 10;
//=============================================================
int pnum, p[kMaxn];
ll mu[kMaxn], sum[kMaxn];
bool vis[kMaxn];
//=============================================================
inline int read() {
  int f = 1, w = 0;
  char ch = getchar();
  for (; !isdigit(ch); ch = getchar())
    if (ch == '-') f = -1;
  for (; isdigit(ch); ch = getchar()) w = (w << 3) + (w << 1) + (ch ^ '0');
  return f * w;
}
void Getmax(int &fir_, int sec_) {
  if (sec_ > fir_) fir_ = sec_;
}
void Getmin(int &fir_, int sec_) {
  if (sec_ < fir_) fir_ = sec_;
}
void Euler(int lim_) {
  vis[1] = true, mu[1] = 1ll;
  for (int i = 2; i <= lim_; ++ i) {
    if (! vis[i]) {
      p[++ pnum] = i;
      mu[i] = - 1;
    }
    for (int j = 1; j <= pnum && i * p[j] <= lim_; ++ j) {
      vis[i * p[j]] = true;
      if (i % p[j] == 0) { //勿忘平方因子
        mu[i * p[j]] = 0;
        break;
      }
      mu[i * p[j]] = - mu[i];
    }
  }
  sum[1] = 1ll;
  for (int i = 1; i <= lim_; ++ i) {
    sum[i] = (sum[i - 1] + 1ll * i * i % kMod * (mu[i] + kMod)) % kMod; //仅令 mu + kMod
  }
}
int g(int n_, int m_) {
  return (1ll * n_ * (n_ + 1ll) / 2ll % kMod) * (1ll * m_ * (m_ + 1ll) / 2ll % kMod) % kMod;
}
int f(int n_, int m_) {
  int lim = std :: min(n_, m_), ret = 0;
  for (int l = 1, r; l <= lim; l = r + 1) {
    r = std :: min(n_ / (n_ / l), m_ / (m_ / l));
    ret = (ret + 1ll * (sum[r] - sum[l - 1] + kMod) * g(n_ / l, m_ / l) % kMod) % kMod;
  }
  return ret;
}
int Sum(ll l, ll r) {
  return (1ll * (r - l + 1ll) * (l + r) / 2ll) % kMod;
}
//=============================================================
int main() { 
  int n = read(), m = read();
  int lim = std :: min(n, m), ans = 0;
  Euler(lim);
  for (int l = 1, r; l <= lim; l = r + 1) {
    r = std :: min(n / (n / l), m / (m / l));
    ans = (ans + 1ll * Sum(l, r) * f(n / l, m / l) % kMod) % kMod;
  }
  printf("%d", ans);
  return 0;
}
/*
7718820 8445880
*/

---

SP5971 LCMSUM - LCM Sum

\(T\) 次询问，每次询问给定 \(n\)，求：

\[\sum_{i=1}^{n}\operatorname{lcm}(i,n) \]

\(1<T\le 3\times 10^5\)，\(1\le n\le 10^6\)。

法一：无脑暴力

先拆 \(\operatorname{lcm}\)，原式等价于：

\[n\sum_{i=1}^{n}\dfrac{i}{\gcd(i,n)} \]

套路的枚举 \(\gcd(i,n)\)，调换枚举顺序，原式等价于：

\[\begin{aligned} &n\sum_{i=1}\sum_{d|i} [\gcd(i,n) = d]\dfrac{i}{d}\\ =& n\sum_{i=1}\sum_{d|i} [\gcd(\dfrac{i}{d},\dfrac{n}{d}) = 1]\dfrac{i}{d}\\ =& n\sum_{d|n}\sum_{i=1}^{n}[d|i][\gcd(\dfrac{i}{d},\dfrac{n}{d}) = 1]\dfrac{i}{d} \end{aligned}\]

把 \(i,n\) 中的 \(d\) 提出来，变为枚举 \(\frac{i}{d}\)，消去整除的条件，原式等价于：

\[n\sum_{d|n}\sum_{i=1}^{\left\lfloor\frac{n}{d}\right\rfloor}[\gcd(i,\dfrac{n}{d}) = 1]i \]

代入 \([\gcd(i,j) = 1] = \sum\limits_{d\mid \gcd(i,j)} {\mu (d)}\)，原式等价于：

\[n\sum_{d|n}\sum_{i=1}^{\left\lfloor\frac{n}{d}\right\rfloor} i \sum_{t|\gcd(i,\frac{n}{d})}\mu (t) \]

值得注意的是 \(t\) 的上界为 \(\frac{n}{d}\)，\(dt\le n\)。
调换枚举顺序，先枚举 \(t\)，原式等价于：

\[n\sum_{d|n}\sum_{t|\frac{n}{d}} \mu(t) \sum_{i=1}^{\left\lfloor\frac{n}{d}\right\rfloor} [t|i] i \]

套路地消去整除的条件，把 \(i\) 中的 \(t\) 提出来，原式等价于：

\[n\sum_{d|n}\sum_{t|\frac{n}{d}} \mu(t)t \sum_{i=1}^{\left\lfloor\frac{n}{dt}\right\rfloor} i \]

对于最后的一个求和项，设 \(g(x) = \sum\limits_{i=1}^{x}i = \frac{x(x+1)}{2}\)，显然可以 \(O(1)\) 求解，原式等价于：

\[n\sum_{d|n}\sum_{t|\frac{n}{d}} \mu(t)t\cdot g(\left\lfloor\dfrac{n}{dt}\right\rfloor) \]

考虑枚举 \(T = dt\)，显然 \(T\le n\)。
\(\mu(t)t\) 与 \(d\) 无关，可以直接考虑枚举 \(t|T\)，原式等价于：

\[n\sum_{T=1}^{n}g(\left\lfloor\dfrac{n}{T}\right\rfloor)\sum_{t|T}\mu(t)t \]

前半块是一个数论分块的形式，可以 \(O(\sqrt{n})\) 求解。
考虑后半块，设 \(f(n)=\sum\limits_{d|n}\mu(d)d\)，发现它是一个积性函数，可以线性筛筛出，具体地：

\[f(n)= \begin{cases} 1-n &n\in \mathrm{primes} \\ f(\frac{x}{p}) &p^2\mid n\\ f(\frac{x}{p})f(p) &p^2\nmid n \end{cases} \]

其中 \(p\) 为 \(n\) 的最小质因子。

此时已经可以线性筛 + 数论分块求解，复杂度 \(O(n+T\sqrt{n})\)，比较菜鸡，时限 500ms 过不了。
考虑筛出 \(f\) 后再用埃氏筛预处理 \(\sum\limits_{T=1}^{n}g(\left\lfloor\dfrac{n}{T}\right\rfloor)f(T)\)，输出时乘上 \(n\)，复杂度变为 \(O(n\log^2 n + n)\)。

---

法二：

同样先拆 \(\operatorname{lcm}\)，枚举 \(\gcd(i,n)\)，调换枚举顺序，原式等价于：

\[\begin{aligned} &n\sum_{i=1}^{n}\dfrac{i}{\gcd(i,n)}\\ =& n\sum_{i=1}\sum_{d|i} [\gcd(i,n) = d]\dfrac{i}{d}\\ =& n\sum_{d|n}\sum_{i=1}^{n}[d|i][\gcd({i},{n}) = d]\dfrac{i}{d} \end{aligned}\]

把 \(i,n\) 中的 \(d\) 提出来，变为枚举 \(\frac{i}{d}\)，消去整除的条件，原式等价于：

\[\begin{aligned} &n\sum_{d|n}\sum_{i=1}^{\left\lfloor\frac{n}{d}\right\rfloor}[\gcd(id,n) = d]i\\ =& n\sum_{d|n}\sum_{i=1}^{\left\lfloor\frac{n}{d}\right\rfloor}[\gcd(i,\dfrac{n}{d}) = 1]i \end{aligned}\]

调整枚举对象，上式等价于：

\[n\sum_{d|n}\sum_{i=1}^{d}[\gcd(i,d) = 1]i \]

考虑 \(\sum\limits_{i=1}^{d}[\gcd(i,d) = 1]i\) 的实际意义，表示 \([1,d]\) 中与 \(d\) 互质的数的和。
当 \(d>1\) 时，与 \(d\) 互质的数总是成对存在，即若 \(\gcd(i,d)=1\) 成立，则 \(\gcd(d-i,d)=1\) 成立。
每对这样的数的和为 \(d\)，共有 \(\frac{\varphi(d)}{2}\) 对这样的数。
则原式等价于：

\[n\sum_{d|n}\dfrac{\varphi(d)d}{2} \]

可以直接预处理答案。
预处理时先线性筛出 \(\varphi\)，再埃氏筛枚举 \(i\) 的倍数，令它们的答案加上 \(\frac{\varphi(i)i}{2}\)，最后输出时乘上 \(n\)。
复杂度 \(O(n\log^2 n + T)\)。

---

法二代码

//知识点：莫比乌斯反演
/*
By:Luckyblock
*/
#include <algorithm>
#include <cctype>
#include <cstdio>
#include <cstring>
#define ll long long
const int kMaxn = 1e6;
//=============================================================
ll phi[kMaxn + 10], ans[kMaxn + 10];
int pnum, p[kMaxn + 10];
bool flag[kMaxn + 10];
//=============================================================
inline int read() {
  int f = 1, w = 0;
  char ch = getchar();
  for (; !isdigit(ch); ch = getchar())
    if (ch == '-') f = -1;
  for (; isdigit(ch); ch = getchar()) w = (w << 3) + (w << 1) + (ch ^ '0');
  return f * w;
}
void GetMax(int &fir, int sec) {
  if (sec > fir) fir = sec;
}
void GetMin(int &fir, int sec) {
  if (sec < fir) fir = sec;
}
void GetPrime() {
  phi[1] = 1, flag[1] = true; //注意初始化
  for (int i = 2; i <= kMaxn; ++ i) {
    if (! flag[i]) {
      p[++ pnum] = i;
      phi[i] = i - 1ll;
    }
    for (int j = 1; j <= pnum && i * p[j] <= kMaxn; ++ j) {
      flag[i * p[j]] = true;
      if (i % p[j]) {
        phi[i * p[j]] = phi[i] * phi[p[j]];
      } else {
        phi[i * p[j]] = phi[i] * p[j];
        break;
      }
    }
  }
  for (int i = 1; i <= kMaxn; ++ i) {
    for (int j = 1; i * j <= kMaxn; ++ j) {
      ans[i * j] += (i == 1 ? 1 : 1ll * phi[i] * i / 2);
    }
  }
}
//=============================================================
int main() { 
  GetPrime();
  int T = read();
  while (T --) {
    int n = read();
    printf("%lld\n", 1ll * ans[n] * n);
  }
  return 0; 
}

---

[SDOI2015]约数个数和

\(T\) 次询问，每次询问给定 \(n,m\)。
定义 >\(\operatorname{d}(i)\) 为 \(i\) 的约数个数，求：

\[\sum_{i=1}^{n}\sum_{j=1}^m\operatorname{d}(ij) \]

\(1<T,n\le 5\times 10^4\)。

一个结论：

\[\operatorname{d}(ij) = \sum_{x|i}\sum_{y|j}[\gcd(x,y)=1] \]

证明

先考虑 \(i = p^a\)，\(j=p^b(p\in \mathrm{primes})\) 的情况，有：

\[\operatorname{d}(p^{a+b})=\sum_{x|p^a}\sum_{y|p^b}[\gcd(x,y)=1] \]

对于等式左侧，\(p^{a+b}\) 的约数个数为 \(a+b+1\)。
对于等式右侧，在保证 \(\gcd(x,y)=1\) 成立的情况下，有贡献的数对 \((x,y)\) 只能是下列三种形式：

• \(x>0,y-0\)，\(x\) 有 \(a\) 种取值方法。
• \(x=0,y>0\)，\(y\) 有 \(b\) 种取值方法。
• \(x=0,y=0\)。

则等式右侧贡献次数为 \(a+b+1\) 次，等于 \(p^{a+b}\) 的约数个数。
则当 \(i = p^a\)，\(j=p^b(p\in \mathrm{primes})\) 时等式成立。

又不同质因子间相互独立，上述情况可拓展到一般的情况。

---

对 \(\operatorname{d}(i,j)\) 进行化简，代入 \([\gcd(i,j) = 1] = \sum\limits_{d\mid \gcd(i,j)} {\mu (d)}\)，原式等价于：

\[\begin{aligned} \operatorname{d}(ij) =& \sum_{x|i}\sum_{y|j}[\gcd(x,y)=1]\\ =& \sum_{x|i}\sum_{y|j}\sum\limits_{d\mid \gcd(x,y)} {\mu (d)} \end{aligned}\]

调换枚举顺序，先枚举 \(d\)，原式等价于：

\[\sum_{d=1}[d|i][d|j]{\mu (d)}\sum_{x|i}[d|x]\sum_{y|j}[d|y] \]

把各项中的 \(d\) 提出来，消去整除的条件，原式等价于：

\[\begin{aligned} &\sum_{d=1}[d|i][d|j]{\mu (d)}\sum_{x|\frac{i}{d}}\sum_{y|\frac{j}{d}}1\\ =& \sum_{d=1}[d|i][d|j]{\mu (d)}\cdot \operatorname{d}(\dfrac{i}{d})\operatorname{d}(\dfrac{j}{d}) \end{aligned}\]

---

将 \(\operatorname{d}(ij)\) 代回原式，原式等价于：

\[\sum_{i=1}^{n}\sum_{j=1}^m \sum_{d=1}[d|i][d|j]{\mu (d)}\cdot \operatorname{d}(\dfrac{i}{d})\operatorname{d}(\dfrac{j}{d}) \]

调换枚举顺序，先枚举 \(d\)，原式等价于：

\[\sum_{d=1}{\mu (d)}\sum_{i=1}^{n}[d|i]\sum_{j=1}^m [d|j]\cdot \operatorname{d}(\dfrac{i}{d})\operatorname{d}(\dfrac{j}{d}) \]

把 \(i,j\) 中的 \(d\) 提出来，变为枚举 \(\frac{i}{d}, \frac{j}{d}\)，消去整除的条件，原式等价于：

\[\sum_{d=1}{\mu (d)}\sum_{i=1}^{\left\lfloor\frac{n}{d}\right\rfloor}\operatorname{d}(i)\sum_{j=1}^{\left\lfloor\frac{m}{d}\right\rfloor}\operatorname{d}(j) \]

考虑预处理 \(S(x) = \sum\limits_{i=1}^{x}\operatorname{d}(i)\)，则原式等价于：

\[\sum_{d=1}{\mu (d)}S\left({\left\lfloor\frac{n}{d}\right\rfloor}\right)S\left({\left\lfloor\frac{m}{d}\right\rfloor}\right) \]

线性筛预处理 \(\mu,\operatorname{d}\)，数论分块求解即可，复杂度 \(O(n+T\sqrt{n})\)。

---

代码

//知识点：莫比乌斯反演
/*
By:Luckyblock
*/
#include <algorithm>
#include <cctype>
#include <cmath>
#include <cstdio>
#include <cstring>
#define ll long long
const int kMaxn = 5e4 + 10;
//=============================================================
int pnum, p[kMaxn];
ll mu[kMaxn], num[kMaxn], d[kMaxn]; //num 为最小质因子的次数
ll summu[kMaxn], sumd[kMaxn];
bool vis[kMaxn];
//=============================================================
inline int read() {
  int f = 1, w = 0;
  char ch = getchar();
  for (; !isdigit(ch); ch = getchar())
    if (ch == '-') f = -1;
  for (; isdigit(ch); ch = getchar()) w = (w << 3) + (w << 1) + (ch ^ '0');
  return f * w;
}
void Getmax(int &fir_, int sec_) {
  if (sec_ > fir_) fir_ = sec_;
}
void Getmin(int &fir_, int sec_) {
  if (sec_ < fir_) fir_ = sec_;
}
void Euler(int lim_) {
  vis[1] = true;
  mu[1] = d[1] = 1ll;
  for (int i = 2; i <= lim_; ++ i) {
    if (! vis[i]) {
      p[++ pnum] = i;
      mu[i] = - 1;
      num[i] = 1;
      d[i] = 2;
    }
    for (int j = 1; j <= pnum && i * p[j] <= lim_; ++ j) {
      vis[i * p[j]] = true;
      if (i % p[j] == 0) { //勿忘平方因子
        mu[i * p[j]] = 0;
        num[i * p[j]] = num[i] + 1;
        d[i * p[j]] = 1ll * d[i] / num[i * p[j]] * (num[i * p[j]] + 1ll);
        break;
      }
      mu[i * p[j]] = - mu[i];
      num[i * p[j]] = 1;
      d[i * p[j]] = 2ll * d[i]; //

    }
  }
  for (int i = 1; i <= lim_; ++ i) {
    summu[i] = summu[i - 1] + mu[i];
    sumd[i] = sumd[i - 1] + d[i];
  }
}

//=============================================================
int main() { 
  Euler(kMaxn - 10);
  int T = read();
  while (T --) {
    int n = read(), m = read(), lim = std :: min(n, m);
    ll ans = 0ll;
    for (int l = 1, r; l <= lim; l = r + 1) {
      r = std :: min(n / (n / l), m / (m / l));
      ans += 1ll * (summu[r] - summu[l - 1]) * sumd[n / l] * sumd[m / l]; //
    }
    printf("%lld\n", ans);
  }
  return 0;
}
/*
in
1
32 43
*/
/*
out
15420
*/

---

P3768 简单的数学题

给定参数 \(n\)、\(p\)，求：

\[\left(\sum_{i=1}^n\sum_{j=1}^ni\cdot j\cdot \gcd(i,j)\right)\bmod p \]

\(n\leq10^{10}\)，\(5\times10^8\leq p\leq1.1\times10^9\)，\(p\in \mathrm{primes}\)。
时限 4s。

无脑套路暴力。

考虑先枚举 \(\gcd(i,j)\)，原式等价于：

\[\begin{aligned} &\sum_{d=1}d\sum_{i=1}^{n}\sum_{j=1}^{n}[\gcd(i,j)=d]ij\\ =& \sum_{d=1}d\sum_{i=1}^{n}[d|i]\sum_{j=1}^{n}[d|j][\gcd(\dfrac{i}{d},\dfrac{j}{d})=1]ij \end{aligned}\]

提出 \(d\)，变为枚举 \(\frac{i}{d}\) 与 \(\frac{j}{d}\)，消去整除的条件，原式等价于：

\[\begin{aligned} &\sum_{d=1}d\sum_{i=1}^{\left\lfloor\frac{n}{d}\right\rfloor}\sum_{j=1}^{\left\lfloor\frac{n}{d}\right\rfloor}[\gcd(i,j)=1]id\cdot jd\\ =& \sum_{d=1} d^3\sum_{i=1}^{\left\lfloor\frac{n}{d}\right\rfloor}\sum_{j=1}^{\left\lfloor\frac{n}{d}\right\rfloor}[\gcd(i,j)=1]ij \end{aligned}\]

代入 \([\gcd(i,j) = 1] = \sum\limits_{d\mid \gcd(i,j)} {\mu (d)}\)，原式等价于：

\[\sum_{d=1} d^3\sum_{i=1}^{\left\lfloor\frac{n}{d}\right\rfloor}\sum_{j=1}^{\left\lfloor\frac{n}{d}\right\rfloor}ij\sum_{t|\gcd(i,j)}\mu (t) \]

值得注意的是 \(t\) 的上界为 \(\frac{n}{d}\)，\(dt\le n\)。
调换枚举顺序，先枚举 \(t\)，原式等价于：

\[\sum_{d=1} d^3\sum_{t=1}\mu(t)\sum_{i=1}^{\left\lfloor\frac{n}{d}\right\rfloor}[t|i] \sum_{j=1}^{\left\lfloor\frac{n}{d}\right\rfloor}[t|j]ij \]

和上面一样，提出 \(t\)，套路地消去整除的条件，原式等价于：

\[\sum_{d=1} d^3\sum_{t=1}\mu(t)t^2\sum_{i=1}^{\left\lfloor\frac{n}{dt}\right\rfloor}\sum_{j=1}^{\left\lfloor\frac{n}{dt}\right\rfloor}ij \]

发现后面两个求和是等差数列乘等差数列的形式。
设 \(g(p,q) = \sum\limits_{i=1}^{p} \sum\limits_{j=1}^{q}ij\)，\(g_{p,q}\) 可以通过下式 \(O(1)\) 计算：

\[g(p,q) = \sum_{i=1}^{p} i \sum_{j=1}^{q}j= \dfrac{(1 + p)\times p}{2}\times \dfrac{(1+q)\times q}{2} \]

代入原式，原式等价于：

\[\sum_{d=1} d^3\sum_{t=1}\mu(t)t^2\cdot g\left({\left\lfloor\frac{n}{dt}\right\rfloor},{\left\lfloor\frac{n}{dt}\right\rfloor}\right) \]

考虑枚举 \(T = dt\)，显然 \(T\le n\)。
再考虑枚举 \(d|T\)，即可得到 \(t = \frac{T}{d}\)，原式等价于：

\[\begin{aligned} &\sum_{T=1}^{n}g\left({\left\lfloor\frac{n}{T}\right\rfloor},{\left\lfloor\frac{n}{T}\right\rfloor}\right)\sum_{d|T}d^3\mu{\left(\dfrac{T}{d}\right)}\left(\dfrac{T}{d}\right)^2\\ =& \sum_{T=1}^{n}T^2 g\left({\left\lfloor\frac{n}{T}\right\rfloor},{\left\lfloor\frac{n}{T}\right\rfloor}\right)\sum_{d|T}d\cdot \mu{\left(\dfrac{T}{d}\right)} \end{aligned}\]

对于后面这一坨，用 \(\sum\limits_{d|T}d\cdot \mu{\left(\frac{T}{d}\right)} = \operatorname{Id} \ast \mu(T)= \varphi(T)\) 反演，则原式等价于：

\[\sum_{T=1}^{n}T^2 \varphi(T) \cdot g\left({\left\lfloor\frac{n}{T}\right\rfloor},{\left\lfloor\frac{n}{T}\right\rfloor}\right) \]

后半块可以数论分块，考虑前半块。
发现前半段即为 \(\operatorname{Id}^2(T)\varphi(T)\)，又是前缀和形式，考虑杜教筛。

有：

\[f(n) = \operatorname{Id}^2\varphi(n) \]

考虑找到一个函数 \(g\)，构造函数 \(h = f\ast g\) 使其便于求值，有：

\[h(n) = \sum_{d|n} d^2\varphi(d)\cdot g\left(\dfrac{n}{d}\right) \]

看到同时存在 \(d\) 和 \(\frac{n}{d}\)，考虑把 \(d^2\) 消去。
令 \(g = \operatorname{Id}^2\)，有：

\[\begin{aligned} h(n) =& \sum_{d|n} d^2\varphi(d)\cdot \left(\dfrac{n}{d}\right)^2\\ =& n^2\sum_{d|n} \varphi(d)\\ =& n^2 \cdot \varphi \ast 1(n)\\ \end{aligned}\]

又 \(\varphi \ast 1 = \operatorname{Id}\)，则有：

\[h(n) = n^3 \]

找到了合适的 \(g\)，套杜教筛的公式。

\[\begin{aligned} g(1)S(n) &= \sum_{i=1}^{n}h(i) - \sum_{i=2}^{n}g(i)S\left(\left\lfloor\dfrac{n}{i}\right\rfloor\right)\\ S(n) &= \sum_{i=1}^{n} i^3 - \sum_{i=2}^{n} i^2\cdot S\left(\left\lfloor\dfrac{n}{i}\right\rfloor\right) \end{aligned}\]

前一项是自然数的立方和，有 \(\sum\limits_{i=1}^{n} i^3 = (\frac{n(n+1)}{2})^2\)。证明详见：自然数前n项平方和、立方和公式及证明 - 百度文库。
后一项直接等差数列求和 + 数论分块求解即可。

---

「SDOI2017」数字表格

设 \(f_{i}\) 表示斐波那契数列的第 \(i\) 项。
\(T\) 组数据，每次给定参数 \(n,m\)，求：

\[\prod_{i=1}^{n}\prod_{j=1}^{m}f_{\gcd(i,j)} \pmod {10^9 + 7} \]

\(1\le T\le 10^3\)，\(1\le n,m\le 10^6\)。
5S，256MB。

以下钦定 \(n\ge m\)。
大力化式子，先套路地枚举 \(\gcd(i,j)\)，用初中知识把两个 \(\prod\) 化到指数位置，原式等于：

\[\large\begin{aligned} &\prod_{d = 1}^{n}\prod_{i=1}^{n}\prod_{j=1}^{m}f_{d}[\gcd(i,j)=d]\\ =&\prod_{d = 1}^{n}f_{d}^{\left(\sum\limits_{i=1}^{n}\sum\limits_{j=1}^{m}[\gcd(i,j)=d]\right)} \end{aligned}\]

分母套路一波，有：

\[\begin{aligned} &\sum\limits_{i=1}^{n}\sum\limits_{j=1}^{m}[\gcd(i,j) = d]\\ =& \sum\limits_{i=1}^{\left\lfloor\frac{n}{d}\right\rfloor}\sum\limits_{j=1}^{\left\lfloor\frac{m}{d}\right\rfloor}[\gcd(i,j) = 1]\\ =& \sum\limits_{i=1}^{\left\lfloor\frac{n}{d}\right\rfloor}\sum\limits_{j=1}^{\left\lfloor\frac{m}{d}\right\rfloor}\sum_{k\mid \gcd(i,j)}\mu (k)\\ =& \sum_{k=1}\mu(k)\left\lfloor\dfrac{n}{kd}\right\rfloor\left\lfloor\dfrac{m}{kd}\right\rfloor \end{aligned}\]

代回原式，原式等于：

\[\large \begin{aligned} &\prod_{d = 1}^{n}f_{d}^{\left(\sum\limits_{k=1}\mu(k)\left\lfloor\frac{n}{kd}\right\rfloor\left\lfloor\frac{m}{kd}\right\rfloor\right)}\\ =& \prod_{d = 1}^{n}\left(f_{d}^{{\sum\limits_{k=1}\mu(k)}}\right)^{\left\lfloor\frac{n}{kd}\right\rfloor\left\lfloor\frac{m}{kd}\right\rfloor} \end{aligned}\]

考虑再暴力拆一波，原式等于：

\[\large \begin{aligned} &\prod_{d = 1}^{n}\left(f_{d}^{{\sum\limits_{k=1}\mu(k)}}\right)^{\left\lfloor\frac{n}{kd}\right\rfloor\left\lfloor\frac{m}{kd}\right\rfloor}\\ =& \prod_{d = 1}^{n}\left(\prod_{k=1}^{\left\lfloor\frac{n}{d}\right\rfloor}f_{d}^{\mu(k)\left\lfloor\frac{n}{kd}\right\rfloor\left\lfloor\frac{m}{kd}\right\rfloor}\right) \end{aligned}\]

做不动了，但发现变量仅有 \(k,d,kd\)，考虑更换枚举对象改为枚举 \(t = kd\) 与 \(d\)，则原式等于：

\[\large\prod_{t=1}^{n}\left(\prod_{d | t}^{n}f_{d}^{{\mu(\frac{t}{d})}}\right)^{\left\lfloor\frac{n}{t}\right\rfloor\left\lfloor\frac{m}{t}\right\rfloor} \]

枚举对象变成了约数形式。从后面的式子推前面的式子是比较显然的，可以认为这种枚举 \(t=kd\) 的形式是一种逆向操作。

设：

\[\large g(t)=\prod_{d | t}^{n}f_{d}^{{\mu(\frac{t}{d})}} \]

\(g(t)\) 可以用类似埃氏筛的方法 \(O(n\log ^2 n)\) 地预处理出来。再把 \(g\) 代回原式，原式等于：

\[\large\prod_{t=1}^{n}g(t)^{\left\lfloor\frac{n}{t}\right\rfloor\left\lfloor\frac{m}{t}\right\rfloor} \]

可以考虑预处理 \(g(t)\) 的前缀积，数论分块枚举指数求解即可。

总时间复杂度 \(O(n\log ^2 n + T\sqrt n)\)，轻微卡常可以过。

//知识点：莫比乌斯反演 
/*
By:Luckyblock
*/
#include <algorithm>
#include <cctype>
#include <cstdio>
#include <cstring>
#define LL long long
const LL mod = 1e9 + 7;
const int kN = 1e6;
//=============================================================
LL n, m, ans;
int p_num, p[kN + 10];
bool vis[kN + 10];
LL mu[kN + 10], f[kN + 10], g[kN + 10], prod[kN + 10];
LL invf[kN + 10], invp[kN];
//=============================================================
inline int read() {
  int f = 1, w = 0;
  char ch = getchar();
  for (; !isdigit(ch); ch = getchar())
    if (ch == '-') f = -1;
  for (; isdigit(ch); ch = getchar()) {
    w = (w << 3) + (w << 1) + (ch ^ '0');
  }
  return f * w;
}
void Chkmax(int &fir_, int sec_) {
  if (sec_ > fir_) fir_ = sec_;
}
void Chkmin(int &fir_, int sec_) {
  if (sec_ < fir_) fir_ = sec_;
}
LL QPow(LL x_, LL y_) {
  x_ %= mod;
  LL ret = 1;
  for (; y_; y_ >>= 1ll) {
    if (y_ & 1) ret = ret * x_ % mod;
    x_ = x_ * x_ % mod;
  }
  return ret;
}
void Euler() {
  vis[1] = true, mu[1] = 1;
  for (int i = 2; i <= kN; ++ i) {
    if (! vis[i]) {
      p[++ p_num] = i;
      mu[i] = -1;
    }
    for (int j = 1; j <= p_num && i * p[j] <= kN; ++ j) {
      vis[i * p[j]] = true;
      if (i % p[j] == 0) {
        mu[i * p[j]] = 0;
        break;
      }
      mu[i * p[j]] = -mu[i];
    }
  }
}
void Prepare() {
  g[1] = g[2] = 1;
  f[1] = f[2] = 1;
  invf[1] = invf[2] = 1;
  for (int i = 3; i <= kN; ++ i) {
    g[i] = 1;
    f[i] = (f[i - 1] + f[i - 2]) % mod;
    invf[i] = QPow(f[i], mod - 2);
  }
  
  Euler();
  for (int d = 1; d <= kN; ++ d) {
    for (int j = 1; d * j <= kN; ++ j) {
      if (mu[j] == 1) {
        g[d * j] = g[d * j] * f[d] % mod;
      } else if (mu[j] == -1) {
        g[d * j] = g[d * j] * invf[d] % mod;
      }
    }
  }
  invp[0] = prod[0] = 1;
  for (int i = 1; i <= kN; ++ i) {
    prod[i] = prod[i - 1] * g[i] % mod;
    invp[i] = QPow(prod[i], mod - 2);
  }
}
//=============================================================
int main() {
  Prepare();
  int T = read();
  while (T -- ){
    n = read(), m = read(), ans = 1;
    if (n < m) std::swap(n, m);
    for (LL l = 1, r = 1; l <= m; l = r + 1) {
      r = std::min(n / (n / l), m / (m / l));
      ans = (ans * QPow(prod[r] * invp[l - 1] % mod, 1ll * (n / l) * (m / l))) % mod;
    }
    printf("%lld\n", ans);
  }
  return 0;
}

---

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

给定参数 \(p\)，有 \(T\) 组数据，每次给定参数 \(A,B,C\)，求：

\[\prod_{i=1}^{A}\prod_{j=1}^{B}\prod_{k=1}^{C}\left(\dfrac{\operatorname{lcm}(i,j)}{\gcd(i,k)}\right)^{f(type)} \]

其中 \(f(type)\) 的取值如下：

\[f(type) = \begin{cases} 1 &type = 0\\ i\times j\times k &type = 1\\ \gcd(i,j,k) &type = 2 \end{cases}\]

\(1\le A,B,C\le 10^5\)，\(10^7\le p\le 1.05\times 10^9\)，\(p\in \mathbb{P}\)，\(T=70\)。
2.5S，128MB。

先化下原式，原式等于：

\[\prod_{i=1}^{A}\prod_{j=1}^{B}\prod_{k=1}^{C}\left(\dfrac{i\times j }{\gcd(i,j)\times \gcd(i,k)}\right)^{f(type)} \]

发现每一项仅与两个变量有关，设：

\[\begin{aligned} f_1(a,b,c) &= \prod_{i=1}^{a}\prod_{j=1}^{b}\prod_{k=1}^{c} i^{f(type)}\\ f_2(a,b,c) &= \prod_{i=1}^{a}\prod_{j=1}^{b}\prod_{k=1}^{c} \gcd(i,j)^{f(type)} \end{aligned}\]

发现 \(\prod\) 可以随意交换，则原式等价于：

\[\dfrac{f_1(a,b,c)\times f_1(b,a,c)}{f_2(a,b,c)\times f_2(a,c,b)} \]

考虑在 \(type\) 取值不同时，如何快速求得 \(f_1\) 与 \(f_2\)。
一共有 \(6\) 个需要推导的式子，不妨就叫它们 \(1\sim 6\) 面了（

---

type = 0

\[\begin{aligned} f_1(a,b,c) &= \prod_{i=1}^{a}\prod_{j=1}^{b}\prod_{k=1}^{c} i\\ f_2(a,b,c) &= \prod_{i=1}^{a}\prod_{j=1}^{b}\prod_{k=1}^{c} \gcd(i,j) \end{aligned}\]

对于 1 面，显然有：

\[\prod_{i=1}^{a}\prod_{j=1}^{b}\prod_{k=1}^{c} i = \left(\prod_{i=1}^{a}i\right)^{b\times c} \]

预处理阶乘 + 快速幂即可，单次计算时间复杂度 \(O(\log n)\)。

---

再考虑 2 面，套路地枚举 \(\gcd\)，显然有：

\[\large \begin{aligned} &\prod_{i=1}^{a}\prod_{j=1}^{b}\prod_{k=1}^{c} \gcd(i,j)\\ =&\left(\prod_{i=1}^{a}\prod_{j=1}^{b}\gcd(i,j)\right)^c\\ =& \left(\prod_{d=1} d^{\left(\sum\limits_{i=1}^{a}\sum\limits_{j=1}^{b}[\gcd(i,j) = d]\right)}\right)^c \end{aligned}\]

指数是个套路，可以看这里：P3455 [POI2007]ZAP-Queries。于是有：

\[\begin{aligned} &\sum\limits_{i=1}^{a}\sum\limits_{j=1}^{b}[\gcd(i,j) = d]\\ =& \sum\limits_{i=1}^{\left\lfloor\frac{a}{d}\right\rfloor}\sum\limits_{j=1}^{\left\lfloor\frac{b}{d}\right\rfloor}[\gcd(i,j) = 1]\\ =& \sum\limits_{i=1}^{\left\lfloor\frac{a}{d}\right\rfloor}\sum\limits_{j=1}^{\left\lfloor\frac{b}{d}\right\rfloor}\sum_{k\mid \gcd(i,j)}\mu (k)\\ =& \sum_{k=1}\mu(k)\left\lfloor\dfrac{a}{kd}\right\rfloor\left\lfloor\dfrac{b}{kd}\right\rfloor \end{aligned}\]

代回原式，略做处理，则原式等于：

\[\large \begin{aligned} &\left(\prod_{d=1} d^{\left(\sum\limits_{k=1}\mu(k)\left\lfloor\frac{a}{kd}\right\rfloor\left\lfloor\frac{b}{kd}\right\rfloor\right)}\right)^c\\ =& \left(\prod_{d=1} \left(d^{\sum\limits_{k=1}\mu(k)}\right)^{\left\lfloor\frac{a}{kd}\right\rfloor\left\lfloor\frac{b}{kd}\right\rfloor}\right)^c\\ =& \prod_{d=1} \left(\prod_{k=1}^{\left\lfloor\frac{n}{d}\right\rfloor}d^{\left(\mu(k)\left\lfloor\frac{a}{kd}\right\rfloor\left\lfloor\frac{b}{kd}\right\rfloor\right)}\right)^c \end{aligned}\]

像「SDOI2017」数字表格 一样，考虑枚举 \(t=kd\) 和 \(d\)，则原式等于：

\[\large \prod_{t=1}^{n}\left(\left(\prod_{d|t} d^{\mu{\left(\frac{t}{d}\right)}}\right)^{\left\lfloor\frac{a}{t}\right\rfloor\left\lfloor\frac{b}{t}\right\rfloor}\right)^c \]

设：

\[\large g_0(t) = \prod_{d|t}d^{\mu\left(\frac{t}{d}\right)} \]

线性筛预处理 \(\mu\) 后，\(g_0(t)\) 可以用埃氏筛预处理，时间复杂度 \(O(n\log n)\)。再代回原式，原式等于：

\[\large \prod_{t=1}^{a}\left(g_0(t)^{\left\lfloor\frac{a}{t}\right\rfloor\left\lfloor\frac{b}{t}\right\rfloor}\right)^c \]

预处理 \(g_0(t)\) 的前缀积和前缀积的逆元，复杂度 \(O(n\log n)\)。
数论分块 + 快速幂计算即可，单次时间复杂度 \(O(\sqrt n\log n)\)。

---

type = 1

\[\begin{aligned} f_1(a,b,c) &= \prod_{i=1}^{a}\prod_{j=1}^{b}\prod_{k=1}^{c} i^{i\times j\times k}\\ f_2(a,b,c) &= \prod_{i=1}^{a}\prod_{j=1}^{b}\prod_{k=1}^{c} \gcd(i,j)^{i\times j\times k} \end{aligned}\]

考虑 3 面，把 \(\prod k\) 扔到指数位置，有：

\[\large \prod_{i=1}^{a}\prod_{j=1}^{b}\prod_{k=1}^{c} i^{i\times j\times k} = \prod_{i=1}^{a}\prod_{j=1}^{b}i^{\left(i\times j\times \sum\limits_{k = 1}^{c} k\right)} \]

再把 \(\prod j\) 也扔到指数位置，引入 \(\operatorname{sum}(n) = \sum_{i=1}^{n} i = \frac{n(n+1)}{2}\)，原式等于：

\[\left(\prod_{i=1}^{a}i^i\right)^{\operatorname{sum}(b)\times \operatorname{sum}(c)} \]

预处理 \(i^i\) 的前缀积，复杂度 \(O(n\log n)\)。
指数可以 \(O(1)\) 算出，然后快速幂，单次时间复杂度 \(O(\log n)\)。

根据费马小定理，指数需要对 \(p - 1\) 取模。注意 \(p-1\) 不是质数，计算 \(\operatorname{sum}\) 时不能用逆元，但乘不爆 LL，直接算就行。

---

再考虑 4 面，发现 \(k\) 与 \(\gcd\) 无关，则同样把 \(\prod k\) 扔到指数位置，则有：

\[\large \begin{aligned} &\prod_{i=1}^{a}\prod_{j=1}^{b}\prod_{k=1}^{c} \gcd(i,j)^{i\times j\times k}\\ =& \left(\prod_{i=1}^a\prod_{j=1}^b\gcd(i,j)^{i\times j}\right)^{\operatorname{sum}(c)} \end{aligned}\]

套路地枚举 \(\gcd\)，原式等于：

\[\large \left(\prod_{d=1}d^{\left(\sum\limits_{i=1}^a \sum\limits_{j=1}^b i\times j[\gcd(i,j)=d]\right)}\right)^{\operatorname{sum}(c)} \]

大力化简指数，有：

\[\large \begin{aligned} &\sum\limits_{i=1}^a \sum\limits_{j=1}^b i\times j[\gcd(i,j)=d]\\ =& d^2 \sum\limits_{i=1}^{\left\lfloor\frac{a}{d}\right\rfloor} \sum\limits_{j=1}^{\left\lfloor\frac{b}{d}\right\rfloor} i\times j[\gcd(i,j)=1\\ =& d^2 \sum\limits_{i=1}^{\left\lfloor\frac{a}{d}\right\rfloor} i \sum\limits_{j=1}^{\left\lfloor\frac{b}{d}\right\rfloor} j\sum\limits_{t|\gcd(i,j)}\mu(t)\\ =& d^2 \sum\limits_{i=1}^{\left\lfloor\frac{a}{d}\right\rfloor} i \sum\limits_{j=1}^{\left\lfloor\frac{b}{d}\right\rfloor} j\sum\limits_{k|\gcd(i,j)}\mu(k)\\ =& d^2 \sum\limits_{k=1}\mu(k)\sum\limits_{i=1}^{\left\lfloor\frac{a}{d}\right\rfloor} i[k|i] \sum\limits_{j=1}^{\left\lfloor\frac{b}{d}\right\rfloor} j[k|j]\\ =& d^2 \sum\limits_{k=1}k^2\mu(k)\sum\limits_{i=1}^{\left\lfloor\frac{a}{kd}\right\rfloor} i\sum\limits_{j=1}^{\left\lfloor\frac{b}{kd}\right\rfloor} j\\ =& d^2\sum\limits_{k=1}k^2\mu(k)\operatorname{sum}{\left(\left\lfloor\frac{a}{kd}\right\rfloor\right)} \operatorname{sum}{\left(\left\lfloor\frac{b}{kd}\right\rfloor\right)}\\ \end{aligned}\]

指数化不动了，代回原式，原式等于：

\[\large \left(\prod_{d=1}d^{\left(d^2\sum\limits_{k=1}k^2\mu(k)\operatorname{sum}{\left(\left\lfloor\frac{a}{kd}\right\rfloor\right)} \operatorname{sum}{\left(\left\lfloor\frac{b}{kd}\right\rfloor\right)}\right)}\right)^{\operatorname{sum}(c)} \]

同 2 面的情况，先展开一下，再枚举 \(t=kd\) 和 \(d\)，原式等于：

\[\large \begin{aligned} &\left(\prod_{d=1}\left(\prod_{k=1}^{\left\lfloor\frac{n}{d}\right\rfloor}d^{\left(d^2 k^2\mu(k)\right)}\right)^{\left(\operatorname{sum}{\left(\left\lfloor\frac{a}{kd}\right\rfloor\right)} \operatorname{sum}{\left(\left\lfloor\frac{b}{kd}\right\rfloor\right)}\right)}\right)^{\operatorname{sum}(c)}\\ =& \prod_{t=1}\left(\left(\prod_{d|t}d^{\left(d^2\left(\frac{t}{d}\right)^2\mu\left(\frac{t}{d}\right)\right)}\right)^{\operatorname{sum}{\left(\left\lfloor\frac{a}{t}\right\rfloor\right)} \operatorname{sum}{\left(\left\lfloor\frac{b}{t}\right\rfloor\right)}}\right)^{\operatorname{sum}(c)}\\ =& \prod_{t=1}\left(\left(\prod_{d|t}d^{\left(t^2\mu\left(\frac{t}{d}\right)\right)}\right)^{\operatorname{sum}{\left(\left\lfloor\frac{a}{t}\right\rfloor\right)} \operatorname{sum}{\left(\left\lfloor\frac{b}{t}\right\rfloor\right)}}\right)^{\operatorname{sum}(c)} \end{aligned}\]

与二面相同，设：

\[\large g_1(t) = \prod_{d|t}d^{\left(t^2\mu\left(\frac{t}{d}\right)\right)} \]

\(g_1(t)\) 可以用埃氏筛套快速幂筛出，时间复杂度 \(O(n\log^2 n)\)。再代回原式，原式等于：

\[\prod_{t=1}\left(g_1(t)^{\operatorname{sum}{\left(\left\lfloor\frac{a}{t}\right\rfloor\right)} \operatorname{sum}{\left(\left\lfloor\frac{b}{t}\right\rfloor\right)}}\right)^{\operatorname{sum}(c)} \]

同样预处理 \(g_1(t)\) 的前缀积及其逆元，时间复杂度 \(O(n\log n)\)。
整除分块 + 快速幂即可，单次时间复杂度 \(O(\sqrt n\log n)\)。

注意指数的取模。

---

type = 2

\[\begin{aligned} f_1(a,b,c) &= \prod_{i=1}^{a}\prod_{j=1}^{b}\prod_{k=1}^{c} i^{\gcd(i,j,k)}\\ f_2(a,b,c) &= \prod_{i=1}^{a}\prod_{j=1}^{b}\prod_{k=1}^{c} \gcd(i,j)^{\gcd(i,j,k)} \end{aligned}\]

考虑 5 面，手段同上，大力反演化简一波，再调换枚举对象，则有：

\[\large \begin{aligned} &\prod_{i=1}^{a}\prod_{j=1}^{b}\prod_{k=1}^{c} i^{\gcd(i,j,k)}\\ =&\prod_{d=1}\prod\limits_{i=1}^{a}i^{\left(\sum\limits_{j=1}^{b}\sum\limits_{k=1}^{c}[\gcd(i,j,k)=d]\right)}\\ =& \prod_{d=1}\prod\limits_{i=1}^{a}i^{\left(\sum\limits_{j=1}^{\left\lfloor\frac{b}{d}\right\rfloor}\sum\limits_{k=1}^{\left\lfloor\frac{c}{d}\right\rfloor}[\gcd(\frac{i}{d},j,k)=1]\right)}\\ =& \prod_{d=1}\prod\limits_{i=1}^{\left\lfloor\frac{a}{d}\right\rfloor}(id)^{\left(d\sum\limits_{j=1}^{\left\lfloor\frac{b}{d}\right\rfloor}\sum\limits_{k=1}^{\left\lfloor\frac{c}{d}\right\rfloor}[\gcd(i,j,k)=1]\right)}\\ =& \prod_{d=1}\prod\limits_{i=1}^{\left\lfloor\frac{a}{d}\right\rfloor}(id)^{\left(d\sum\limits_{j=1}^{\left\lfloor\frac{b}{d}\right\rfloor}\sum\limits_{k=1}^{\left\lfloor\frac{c}{d}\right\rfloor}\sum\limits_{x|\gcd(i,j,k)}{\mu(x)}\right)}\\ =& \prod_{d=1}\prod\limits_{i=1}^{\left\lfloor\frac{a}{d}\right\rfloor}(id)^{\left(d\sum\limits_{x=1}\mu(x)[x|i]\sum\limits_{j=1}^{\left\lfloor\frac{b}{d}\right\rfloor}[x|j]\sum\limits_{k=1}^{\left\lfloor\frac{c}{d}\right\rfloor}[x|k]\right)}\\ =& \prod_{d=1}\prod_{x=1}\prod\limits_{i=1}^{\left\lfloor\frac{a}{d}\right\rfloor}(id)^{\left(d\times \mu(x)[x|i]\sum\limits_{j=1}^{\left\lfloor\frac{b}{d}\right\rfloor}[x|j]\sum\limits_{k=1}^{\left\lfloor\frac{c}{d}\right\rfloor}[x|k]\right)}\\ =& \prod_{d=1}\prod_{x=1}\prod\limits_{i=1}^{\left\lfloor\frac{a}{xd}\right\rfloor}(ixd)^{\left(d\times \mu(x){\left\lfloor\frac{b}{xd}\right\rfloor}{\left\lfloor\frac{c}{xd}\right\rfloor}\right)}\\ =& \prod_{t = 1}\prod_{d|T}\prod_{i=1}^{\left\lfloor\frac{a}{t}\right\rfloor}(it)^{\left(d\times \mu\left(\frac{t}{d}\right){\left\lfloor\frac{b}{t}\right\rfloor}{\left\lfloor\frac{c}{t}\right\rfloor}\right)}\\ =& \prod_{t = 1}\prod_{d|T}\left(t^{\left\lfloor\frac{a}{t}\right\rfloor}\prod_{i=1}^{\left\lfloor\frac{a}{t}\right\rfloor}i\right)^{d\times \mu\left(\frac{t}{d}\right){\left\lfloor\frac{b}{t}\right\rfloor}{\left\lfloor\frac{c}{t}\right\rfloor}}\\ \end{aligned}\]

引入 \(\operatorname{fac}(n) = \prod_{i=1}^{n} i\)，再根据枚举对象调整一下指数，原式等于：

\[\large \begin{aligned} &\prod_{t = 1}\prod_{d|t}\left(t^{\left\lfloor\frac{a}{t}\right\rfloor}\times \operatorname{fac}\left(\left\lfloor\frac{a}{t}\right\rfloor\right)\right)^{\left(d\times \mu\left(\frac{t}{d}\right){\left\lfloor\frac{b}{t}\right\rfloor}{\left\lfloor\frac{c}{t}\right\rfloor}\right)}\\ =& \prod_{t = 1}\left(\prod_{d|t}\left(t^{\left\lfloor\frac{a}{t}\right\rfloor}\times \operatorname{fac}\left(\left\lfloor\frac{a}{t}\right\rfloor\right)\right)^{d\times \mu\left(\frac{t}{d}\right)}\right)^{{\left\lfloor\frac{b}{t}\right\rfloor}{\left\lfloor\frac{c}{t}\right\rfloor}}\\ =& \prod_{t = 1}\left(\left(t^{\left\lfloor\frac{a}{t}\right\rfloor}\times \operatorname{fac}\left(\left\lfloor\frac{a}{t}\right\rfloor\right)\right)^{\sum\limits_{d|t}d\times \mu\left(\frac{t}{d}\right)}\right)^{{\left\lfloor\frac{b}{t}\right\rfloor}{\left\lfloor\frac{c}{t}\right\rfloor}} \end{aligned}\]

指数中出现了一个经典的狄利克雷卷积的形式，对其进行反演。
将 \((\operatorname{Id}\ast \mu) (n)= \varphi (n)\) 代入原式，原式等于：

\[\large \begin{aligned} &\prod_{t = 1}\left(t^{\left\lfloor\frac{a}{t}\right\rfloor}\times \operatorname{fac}\left(\left\lfloor\frac{a}{t}\right\rfloor\right)\right)^{\varphi(t){\left\lfloor\frac{b}{t}\right\rfloor}{\left\lfloor\frac{c}{t}\right\rfloor}}\\ =& \prod_{t = 1}\left(t^{\varphi(t)\left\lfloor\frac{a}{t}\right\rfloor{\left\lfloor\frac{b}{t}\right\rfloor}{\left\lfloor\frac{c}{t}\right\rfloor}}\times \operatorname{fac}\left(\left\lfloor\frac{a}{t}\right\rfloor\right)^{\varphi(t){\left\lfloor\frac{b}{t}\right\rfloor}{\left\lfloor\frac{c}{t}\right\rfloor}}\right) \end{aligned}\]

预处理 \(t^{\varphi(t)}\) 的前缀积及逆元，阶乘的前缀积及阶乘逆元，\(\pmod {p-1}\) 下的 \(\varphi(t)\) 的前缀和（指数
），时间复杂度 \(O(n\log n)\)。
同样整除分块 + 快速幂即可，单次时间复杂度 \(O(\sqrt n\log n)\)。

---

然后是最掉 sans 的 6 面。有 \(\gcd(i,j,k) = \gcd(\gcd(i,j), k)\)，考虑先枚举 \(\gcd(i,j)\)，然后套路化式子，则有：

\[\large \begin{aligned} &\prod_{i=1}^{a}\prod_{j=1}^{b}\prod_{k=1}^{c} \gcd(i,j)^{\gcd(i,j,k)}\\ =& \prod_{d=1}\prod_{i=1}^{a}\prod_{j=1}^{b}\prod_{k=1}^{c} [\gcd(i,j)=d] d^{\gcd(d,k)}\\ =& \prod_{d=1} \left(d^{\left(\sum\limits_{i=1}^{a}\sum\limits_{j=1}^{b} [\gcd(i,j)=d]\right)}\right)^{\sum\limits_{k=1}^{c}\gcd(d,k)} \end{aligned}\]

先考虑最外面的指数，这也是个套路，可以参考 一个例子。用 \(\operatorname{Id} = \varphi \ast 1\) 反演，显然有：

\[\large \begin{aligned} &\sum\limits_{k=1}^{c}\gcd(d,k)\\ =& \sum\limits_{k=1}^{c}\sum_{x|\gcd(d,k)}\varphi(x)\\ =& \sum_{x=1}\varphi(x)[x|d]\sum_{k=1}^{c}[x|k]\\ =& \sum_{x|d}\varphi(x)\left\lfloor\frac{c}{x}\right\rfloor \end{aligned}\]

再考虑里面的指数，发现这式子在 2 面已经推了一遍了，于是直接拿过来用，有：

\[\sum\limits_{i=1}^{a}\sum\limits_{j=1}^{b}[\gcd(i,j) = d]=\sum_{y=1}\mu(y)\left\lfloor\dfrac{a}{yd}\right\rfloor\left\lfloor\dfrac{b}{yd}\right\rfloor \]

将化简后的两个指数代入原式，原式等于：

\[\large \begin{aligned} &\prod_{d=1} \left(d^{\left(\sum\limits_{y=1}\mu(y)\left\lfloor\frac{a}{yd}\right\rfloor\left\lfloor\frac{b}{yd}\right\rfloor\right)}\right)^{\sum\limits_{x|d}\varphi(x)\left\lfloor\frac{c}{x}\right\rfloor}\\ =& \prod_{d=1} \left(\prod\limits_{y=1}d^{\left(\mu(y)\left\lfloor\frac{a}{yd}\right\rfloor\left\lfloor\frac{b}{yd}\right\rfloor\right)}\right)^{\sum\limits_{x|d}\varphi(x)\left\lfloor\frac{c}{x}\right\rfloor} \end{aligned}\]

与 2、4 面同样套路地，考虑枚举 \(t=yd\) 和 \(d\)，再略作调整，原式等于：

\[\large \begin{aligned} &\prod_{d=1} \left(\prod\limits_{y=1}d^{\left(\mu(y)\left\lfloor\frac{a}{yd}\right\rfloor\left\lfloor\frac{b}{yd}\right\rfloor\right)}\right)^{\sum\limits_{x|d}\varphi(x)\left\lfloor\frac{c}{x}\right\rfloor}\\ =& \prod_{t=1}\prod_{d|t} d^{\left(\mu\left(\frac{t}{d}\right)\left\lfloor\frac{a}{t}\right\rfloor\left\lfloor\frac{b}{t}\right\rfloor\sum\limits_{x|d}\varphi(x)\left\lfloor\frac{c}{x}\right\rfloor\right)}\\ =& \prod_{t=1}\left(\prod_{d|t} d^{\left(\mu\left(\frac{t}{d}\right)\sum\limits_{x|d}\varphi(x)\left\lfloor\frac{c}{x}\right\rfloor\right)}\right)^{\left\lfloor\frac{a}{t}\right\rfloor\left\lfloor\frac{b}{t}\right\rfloor}\\ =& \prod_{t=1}\left(\prod_{d|t} \prod_{x|d}d^{\left(\mu\left(\frac{t}{d}\right)\varphi(x)\left\lfloor\frac{c}{x}\right\rfloor\right)}\right)^{\left\lfloor\frac{a}{t}\right\rfloor\left\lfloor\frac{b}{t}\right\rfloor} \end{aligned}\]

发现要同时枚举 \(d\) 和 \(x\)，化不动了。
从题解里学到一个比较神的技巧，考虑把 \(d\) 拆成 \(x\) 和 \(\frac{d}{x}\) 分别计算贡献再相乘，即分别计算下两式：

\[\large \begin{aligned} &\prod_{t=1}\left(\prod_{d|t} \prod_{x|d}x^{\left(\mu\left(\frac{t}{d}\right)\varphi(x)\left\lfloor\frac{c}{x}\right\rfloor\right)}\right)^{\left\lfloor\frac{a}{t}\right\rfloor\left\lfloor\frac{b}{t}\right\rfloor}\\ &\prod_{t=1}\left(\prod_{d|t} \prod_{x|d}{\left(\frac{d}{x}\right)}^{\left(\mu\left(\frac{t}{d}\right)\varphi(x)\left\lfloor\frac{c}{x}\right\rfloor\right)}\right)^{\left\lfloor\frac{a}{t}\right\rfloor\left\lfloor\frac{b}{t}\right\rfloor} \end{aligned}\]

---

先考虑 \(x\) 的情况，首先把枚举 \(x\) 调整到最外层。设 \(\operatorname{lim}=\max(a,b,c)\)，则原式等于：

\[\large \begin{aligned} &\prod_{x=1} \prod_{t=1}^{\operatorname{lim}}[x|t]\left(\prod_{d|t} [x|d]{x}^{\left(\mu\left(\frac{t}{d}\right)\varphi(x)\left\lfloor\frac{c}{x}\right\rfloor\right)}\right)^{\left\lfloor\frac{a}{t}\right\rfloor\left\lfloor\frac{b}{t}\right\rfloor}\\ =& \prod_{x=1} \prod_{t=1}^{\left\lfloor\frac{\operatorname{lim}}{x}\right\rfloor}\left(\prod_{d|t} {x}^{\left(\mu\left(\frac{tx}{dx}\right)\varphi(x)\left\lfloor\frac{c}{x}\right\rfloor\right)}\right)^{\left\lfloor\frac{a}{tx}\right\rfloor\left\lfloor\frac{b}{tx}\right\rfloor}\\ =& \prod_{x=1} \prod_{t=1}^{\left\lfloor\frac{\operatorname{lim}}{x}\right\rfloor}\prod_{d|t} {x}^{\left(\varphi(x)\left\lfloor\frac{c}{x}\right\rfloor\left\lfloor\frac{a}{tx}\right\rfloor\left\lfloor\frac{b}{tx}\right\rfloor\mu\left(\frac{t}{d}\right)\right)} \end{aligned}\]

把 \(\prod {t}\) 挪到指数位置，原式等于：

\[\large \begin{aligned} &\prod_{x=1} {x}^{\left(\sum\limits_{t=1}^{\left\lfloor\frac{\operatorname{lim}}{x}\right\rfloor}\sum\limits_{d|t}\varphi(x)\left\lfloor\frac{c}{x}\right\rfloor\left\lfloor\frac{a}{tx}\right\rfloor\left\lfloor\frac{b}{tx}\right\rfloor\mu\left(\frac{t}{d}\right)\right)}\\ =& \prod_{x=1} {x}^{\left(\varphi(x)\left\lfloor\frac{c}{x}\right\rfloor\sum\limits_{t=1}^{\left\lfloor\frac{\operatorname{lim}}{x}\right\rfloor}\left\lfloor\frac{a}{tx}\right\rfloor\left\lfloor\frac{b}{tx}\right\rfloor\sum\limits_{d|t}\mu\left(\frac{t}{d}\right)\right)} \end{aligned}\]

指数中又出现了一个经典的狄利克雷卷积的形式，对其进行反演。
将 \((\mu \ast 1) (n)= \epsilon (n)=[n=1]\) 代入原式，原式等于：

\[\large \begin{aligned} &\prod_{x=1} {x}^{\left(\varphi(x)\left\lfloor\frac{c}{x}\right\rfloor\sum\limits_{t=1}^{\left\lfloor\frac{\operatorname{lim}}{x}\right\rfloor}\left\lfloor\frac{a}{tx}\right\rfloor\left\lfloor\frac{b}{tx}\right\rfloor[t=1]\right)}\\ =& \prod_{x=1} {x}^{\left(\varphi(x)\left\lfloor\frac{a}{x}\right\rfloor\left\lfloor\frac{b}{x}\right\rfloor\left\lfloor\frac{c}{x}\right\rfloor\right)} \end{aligned}\]

得到了一个非常优美的式子，而且发现这个式子是 5 面最终答案的一部分。同 5 面的做法，直接整除分块即可。

---

再考虑 \(\frac{d}{x}\) 的情况，同上先把枚举 \(x\) 放到最外层，并调整一下指数，则原式等于：

\[\large \begin{aligned} &\prod_{x=1} \prod_{t=1}^{\operatorname{lim}}[x|t]\left(\prod_{d|t} [x|d]{\left(\frac{d}{x}\right)}^{\left(\mu\left(\frac{t}{d}\right)\varphi(x)\left\lfloor\frac{c}{x}\right\rfloor\right)}\right)^{\left\lfloor\frac{a}{t}\right\rfloor\left\lfloor\frac{b}{t}\right\rfloor}\\ =& \prod_{x=1} \prod_{t=1}^{\left\lfloor\frac{\operatorname{lim}}{x}\right\rfloor}\left(\prod_{d|tx} [x|d]{\left(\frac{d}{x}\right)}^{\left(\mu\left(\frac{tx}{d}\right)\varphi(x)\left\lfloor\frac{c}{x}\right\rfloor\right)}\right)^{\left\lfloor\frac{a}{tx}\right\rfloor\left\lfloor\frac{b}{tx}\right\rfloor}\\ =& \prod_{x=1} \left(\prod_{t=1}^{\left\lfloor\frac{\operatorname{lim}}{x}\right\rfloor}\left(\prod_{d|tx} [x|d]{\left(\frac{d}{x}\right)}^{\mu\left(\frac{tx}{d}\right)}\right)^{\left\lfloor\frac{a}{tx}\right\rfloor\left\lfloor\frac{b}{tx}\right\rfloor}\right)^{\varphi(x)\left\lfloor\frac{c}{x}\right\rfloor} \end{aligned}\]

考虑枚举 \(dx\)，替换原来的 \(d\)，注意一下这里的倍数关系。原式等于：

\[\large \prod_{x=1} \left(\prod_{t=1}^{\left\lfloor\frac{\operatorname{lim}}{x}\right\rfloor}\left(\prod_{d|t}d^{\mu\left(\frac{t}{d}\right)}\right)^{\left\lfloor\frac{a}{tx}\right\rfloor\left\lfloor\frac{b}{tx}\right\rfloor}\right)^{\varphi(x)\left\lfloor\frac{c}{x}\right\rfloor} \]

发现最内层的式子 \(\prod_{d|t}d^{\mu\left(\frac{t}{d}\right)}\)，即为二面处理过的 \(g_0(t)\)。直接代入，原式等于：

\[\large \prod_{x=1} \left(\prod_{t=1}^{\left\lfloor\frac{\operatorname{lim}}{x}\right\rfloor}g_0(t)^{\left\lfloor\frac{a}{tx}\right\rfloor\left\lfloor\frac{b}{tx}\right\rfloor}\right)^{\varphi(x)\left\lfloor\frac{c}{x}\right\rfloor} \]

一个小结论，证明可以看 这里：

\[\forall a,b,c\in \mathbb{Z},\left\lfloor\dfrac{a}{bc}\right\rfloor = \left\lfloor{\dfrac{\left\lfloor\dfrac{a}{b}\right\rfloor}{c}}\right\rfloor \]

则原式等于：

\[\large \prod_{x=1} \left(\prod_{t=1}^{\left\lfloor\frac{\operatorname{lim}}{x}\right\rfloor}g_0(t)^{\left\lfloor\frac{\left\lfloor\frac{a}{x}\right\rfloor}{t}\right\rfloor\left\lfloor\frac{\left\lfloor\frac{b}{x}\right\rfloor}{t}\right\rfloor}\right)^{\varphi(x)\left\lfloor\frac{c}{x}\right\rfloor} \]

于是可以先对外层整除分块，再对内层整除分块。

然后就做完了，哈哈哈。

---

一些实现上的小技巧：

• 逆元能预处理就预处理。
• 注意对指数取模，模数为 \(p-1\)。

//知识点：莫比乌斯反演 
/*
By:Luckyblock
用了比较清晰易懂的变量名，阅读体验应该会比较好。  
vsc 的自动补全真是太好啦！
*/
#include <algorithm>
#include <cctype>
#include <cstdio>
#include <cstring>
using std::min;
using std::max;
#define LL long long
const int Lim = 1e5;
const int kN = 1e5 + 10;
//=============================================================
LL A, B, C, mod, ans;
int T, p_num, p[kN];
bool vis[kN];
LL mu[kN], phi[kN], fac[kN], g[2][kN];
LL sumphi[kN], prodt_phi[kN], prodi_i[kN], prodg[2][kN];
LL inv[kN], inv_fac[kN], inv_prodt_phi[kN], inv_prodg[2][kN];
//=============================================================
inline int read() {
  int f = 1, w = 0;
  char ch = getchar();
  for (; !isdigit(ch); ch = getchar())
    if (ch == '-') f = -1;
  for (; isdigit(ch); ch = getchar()) {
    w = (w << 3) + (w << 1) + (ch ^ '0');
  }
  return f * w;
}
void Chkmax(int &fir_, int sec_) {
  if (sec_ > fir_) fir_ = sec_;
}
void Chkmin(int &fir_, int sec_) {
  if (sec_ < fir_) fir_ = sec_;
}
LL QPow(LL x_, LL y_) {
  x_ %= mod;
  y_ %= mod - 1;
  LL ret = 1;
  for (; y_; y_ >>= 1ll) {
    if (y_ & 1) ret = ret * x_ % mod;
    x_ = x_ * x_ % mod;
  }
  return ret;
}
LL Inv(LL x_) {
  return QPow(x_, mod - 2);
}
LL Sum(LL n_) {
  return ((n_ * (n_ + 1ll)) / 2ll) % (mod - 1);
}
void Euler() {
  vis[1] = true, mu[1] = phi[1] = 1; //初值
  for (int i = 2; i <= Lim; ++ i) {
    if (! vis[i]) {
      p[++ p_num] = i;
      mu[i] = -1;
      phi[i] = i - 1;
    }
    for (int j = 1; j <= p_num && i * p[j] <= Lim; ++ j) {
      vis[i * p[j]] = true;
      if (i % p[j] == 0) {
        mu[i * p[j]] = 0;
        phi[i * p[j]] = phi[i] * p[j];
        break;
      }
      mu[i * p[j]] = -mu[i];
      phi[i * p[j]] = phi[i] * (p[j] - 1);
    }
  }
}
void Prepare() {
  Euler();
  inv[1] = fac[0] = prodt_phi[0] = prodi_i[0] = 1;
  for (int i = 1; i <= Lim; ++ i) {
    g[0][i] = g[1][i] = 1;
    fac[i] = 1ll * fac[i - 1] * i % mod;
    sumphi[i] = (sumphi[i - 1] + phi[i]) % (mod - 1);
    prodi_i[i] = prodi_i[i - 1] * QPow(i, i) % mod;
    if (i > 1) inv[i] = (mod - mod / i) * inv[mod % i] % mod;

    prodt_phi[i] = prodt_phi[i - 1] * QPow(i, phi[i]) % mod;
    inv_prodt_phi[i] = Inv(prodt_phi[i]);
  }

  for (int d = 1; d <= Lim; ++ d) {
    for (int j = 1; d * j <= Lim; ++ j) {
      int t = d * j;
      if (mu[j] == 1) {
        g[0][t] = g[0][t] * d % mod;
        g[1][t] = g[1][t] * QPow(1ll * d, 1ll * t * t) % mod;
      } else if (mu[j] == -1) {
        g[0][t] = g[0][t] * inv[d] % mod;
        g[1][t] = g[1][t] * Inv(QPow(1ll * d, 1ll * t * t)) % mod;
      }
    }
  }
  inv_prodg[0][0] = prodg[0][0] = 1;
  inv_prodg[1][0] = prodg[1][0] = 1;
  inv_prodt_phi[0] = 1;
  for (int i = 1; i <= Lim; ++ i) {
    for (int j = 0; j <= 1; ++ j) {
      prodg[j][i] = prodg[j][i - 1] * g[j][i] % mod;
      inv_prodg[j][i] = Inv(prodg[j][i]);
    }
  }
}
LL f1(LL a_, LL b_, LL c_, int type) {
  if (! type) return QPow(fac[a_], b_ * c_);
  if (type == 1) return QPow(prodi_i[a_], Sum(b_) * Sum(c_));
  LL ret = 1, lim = min(min(a_, b_), c_);
  for (LL l = 1, r = 1; l <= lim; l = r + 1) {
    r = min(min(a_ / (a_ / l), b_ / (b_ / l)), c_ / (c_ / l));
    ret = ret * QPow(prodt_phi[r] * inv_prodt_phi[l - 1], (a_ / l) * (b_ / l) % (mod - 1) * (c_ / l)) % mod;
    ret = ret * QPow(fac[a_ / l], (sumphi[r] - sumphi[l - 1] + mod - 1) % (mod - 1) * (b_ / l) % (mod - 1) * (c_ / l)) % mod;
  }
  return ret;
}
LL f2_2(LL a_, LL b_) { 
  LL ret = 1;
  for (LL l = 1, r = 1; l <= min(a_, b_); l = r + 1) {
    r = min(a_ / (a_ / l), b_ / (b_ / l));
    ret = ret * QPow(prodg[0][r] * inv_prodg[0][l - 1], 1ll * (a_ / l) * (b_ / l)) % mod;
  }
  return ret;
}
LL f2(LL a_, LL b_, LL c_, int type) {
  LL ret = 1;
  if (! type) {
    for (LL l = 1, r = 1; l <= min(a_, b_); l = r + 1) {
      r = min(a_ / (a_ / l), b_ / (b_ / l));
      LL val = QPow(prodg[0][r] * inv_prodg[0][l - 1], 1ll * (a_ / l) * (b_ / l));
      ret = (ret * QPow(val, c_)) % mod;
    }
  } else if (type == 1) {
    for (LL l = 1, r = 1; l <= min(a_, b_); l = r + 1) {
      r = min(a_ / (a_ / l), b_ / (b_ / l));
      LL val = QPow(prodg[1][r] * inv_prodg[1][l - 1], Sum(a_ / l) * Sum(b_ / l));
      ret = (ret * QPow(val, Sum(c_))) % mod;
    }
  } else {
    LL lim = min(min(a_, b_), c_);
    for (LL l = 1, r = 1; l <= lim; l = r + 1) {
      r = min(min(a_ / (a_ / l), b_ / (b_ / l)), c_ / (c_ / l));
      ret = ret * QPow(f2_2(a_ / l, b_ / l), (sumphi[r] - sumphi[l - 1] + mod - 1) % (mod - 1) * (c_ / l)) % mod;
      ret = ret * QPow(prodt_phi[r] * inv_prodt_phi[l - 1], (a_ / l) * (b_ / l) % (mod - 1) * (c_ / l)) % mod;
    }
  }
  return ret;
}
//=============================================================
int main() {
  T = read(), mod = read();
  Prepare();
  while (T -- ) {
    A = read(), B = read(), C = read();
    for (int i = 0; i <= 2; ++ i) {
      ans = f1(A, B, C, i) * f1(B, A, C, i) % mod;
      ans = ans * Inv(f2(A, B, C, i)) % mod * Inv(f2(A, C, B, i)) % mod;
      printf("%lld ", ans);  
    }
    printf("\n");
  }
  return 0;
}

---

写在最后

参考资料：
Oi-Wiki-莫比乌斯反演
算法学习笔记(35): 狄利克雷卷积 By: Pecco
题解 SP5971 【LCMSUM - LCM Sum】 - BJpers2 的博客
题解 SP5971 【LCMSUM - LCM Sum】 - Venus 的博客
题解 P3327 【[SDOI2015]约数个数和】 - suncongbo 的博客

把 Oi-Wiki 上的内容进行了复制 整理扩充。
我是个没有感情的复读机(大雾)