「笔记」莫比乌斯反演

目录

• 前置知识
  • 小碎骨
  • 引理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] = egin{cases} 1 & ext{If P is true}\ 0 & ext{Otherwise} end{cases})
此处 (P) 是一个可真可假的命题。

引理1

[forall a,b,cin mathbb{Z},leftlfloordfrac{a}{bc} ight floor = leftlfloor{dfrac{leftlfloordfrac{a}{b} ight floor}{c}} ight floor ]

证明

[dfrac{a}{b} = leftlfloor{dfrac{a}{b}} ight floor + r(0le r < 1) ]

[leftlfloordfrac{a}{bc} ight floor = leftlfloordfrac{a}{b} imesdfrac{1}{c} ight floor = leftlfloor{dfrac{1}{c} imesleft({leftlfloor{dfrac{a}{b} + r} ight floor} ight)} ight floor = leftlfloor{dfrac{leftlfloordfrac{a}{b} ight floor}{c} + dfrac{r}{c}} ight floor = leftlfloor{dfrac{leftlfloor{dfrac{a}{b}} ight floor}{c}} ight floor ]

数论分块

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

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

对于任意一个(i)，最大的满足上式的 (j=leftlfloor{dfrac{n}{leftlfloor{dfrac{n}{i}} ight floor}} ight floor)

---

积性函数

定义

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

性质

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

[egin{aligned} & h(x) = f(x^p)\ & h(x) = f^p(x)\ & h(x) = f(x)g(x)\ & h(x) = sum_{dmid x} f(d)gleft(dfrac{x}{d} ight) 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) = sumlimits_{dmid n} 1)。
• 除数函数 (sigma_{k}(n) = sumlimits_{dmid n} d^k)。
  (k=1) 时为因数和函数，通常简记为 (sigma(n))。
  (k=0) 时为因数个数函数 (sigma_{0}(n))。
• 欧拉函数 (varphi(n) = sumlimits_{i=1}^{n} [gcd(i,n) = 1])。
• 莫比乌斯函数 (mu(n) = egin{cases}1 &n=1\0 &n ext{含有平方因子}\(-1)^k &k ext{为} n ext{的本质不同质因子个数} end{cases})

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

---

莫比乌斯函数

定义

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

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

解释

令 (n = prodlimits_{i=1}^{k} p_{i}^{c_i})，其中(p_i)为质因子，(c_ige 1)。

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

性质

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

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

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

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

反演常用结论

一个反演常用结论:

[[gcd(i,j) = 1] = sumlimits_{dmid gcd(i,j)} {mu (d)} ]

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

证明 2：
暴力展开：([gcd(i,j) = 1] = e(gcd(i,j)) = sumlimits_{dmid 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(fast g) (n) = sum_{dmid n} f(d)gleft(dfrac{n}{d} ight) ]

性质

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

关于单位元 (e)

有：

[e = mu ast 1=sumlimits_{dmid n} mu (d) ]

证明

[egin{aligned} (fast e)(n) = sum_{dmid n} f(d)e(dfrac{n}{d}) = sum_{dmid n} f(d)left[dfrac{n}{d} = 1 ight] end{aligned}]

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

综上，((fast e)(n) = sumlimits_{dmid n} f(d)left[dfrac{n}{d} = 1 ight] = f(n))，满足单位元的性质。
则 (e = mu ast 1) 成立。

除数函数与幂函数

幂函数 (operatorname{Id}_{k}(n) = n^k)。
除数函数 (sigma_{k}(n) = sumlimits_{dmid n} d^k)。

显然有：

[(operatorname{Id}_kast 1)(n) = sum_{dmid n} operatorname{Id_k}(d) = sum_{dmid n} d^k = sigma_k(n) ]

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

[(1ast 1)(n) = sum_{dmid n}1 = sigma_0(n) ]

欧拉函数与恒等函数

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

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

[(varphi ast 1)(p^m) = sum_{dmid 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}) 个，为 (1sim 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 ight) = prodleft(varphi ast 1 ight)(p^m) = prod p^m = n ]

得证。

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

计算

[(fast g) (n) = sum_{dmid n} f(d)gleft(dfrac{n}{d} ight) ]

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

---

莫比乌斯反演

反演是个啥？反演。

公式

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

[large f(n) = (gast 1)(n) = sumlimits_{dmid n}{g(d)} ]

那么有

[large g(n) = (mu ast f)(n)=sumlimits_{dmid n} {mu(d)fleft(dfrac{n}{d} ight)} ]

证明

方法一：运用卷积。

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

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

1. 用 (sumlimits_{dmid n}g(d)) 替换(fleft(dfrac{n}{d} ight))。

   [sum_{dmid n}{mu(d)sum_{kmid frac{n}{d}}g(k)} ]

2. 变换求和顺序。

   [sum_{kmid n}g(k)sum_{dmid frac{n}{k}}{mu(d)} ]

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

举例

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

[e = mu ast 1iffmu = muast e = mu ]

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

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

---

关于一类求和式

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

一般套路

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

[f(n) = gast 1 = sum_{dmid n}g(d) ]

则求和式变为：

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

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

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

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

[sum_{d=1}g(d)leftlfloordfrac{n}{d} ight floorleftlfloordfrac{m}{d} ight floor ]

数论分块求解即可。

例 1

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

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

[egin{aligned} &sum_{i=1}^{n}sum_{j=1}^{m}operatorname{Id}(gcd(i,j))\ =& sum_{i=1}^{n}sum_{j=1}^{m}sum_{dmid gcd(i,j)}varphi(d)\ =& sum_{d=1}varphi(d)leftlfloordfrac{n}{d} ight floorleftlfloordfrac{m}{d} ight floor end{aligned}]

例 2

[egin{aligned} &sumlimits_{i=1}^{n}sumlimits_{j=1}^{m}[gcd(i,j) = 1]\ =& sumlimits_{i=1}^{n}sumlimits_{j=1}^{m}e(gcd(i,j))\ =& sumlimits_{i=1}^{n}sumlimits_{j=1}^{m}sum_{dmid gcd(i,j)}mu (d)\ =& sum_{d=1}mu(d)leftlfloordfrac{n}{d} ight floorleftlfloordfrac{m}{d} ight floor end{aligned}]

感悟

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

---

例题

[HAOI2011] Problem b

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

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

(1le n,k,a,b,c,dle 5 imes 10^4)，(ale b,cle d)。

令 (f(n,m) = sumlimits_{i=1}^{n}sumlimits_{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)。

易知原式等价于

[sumlimits_{i=1}^{leftlfloorfrac{n}{k} ight floor}sumlimits_{j=1}^{leftlfloorfrac{m}{k} ight floor}[gcd(i,j) = 1] ]

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

[sumlimits_{i=1}^{leftlfloorfrac{n}{k} ight floor}sumlimits_{j=1}^{leftlfloorfrac{m}{k} ight floor}sum_{dmid gcd(i,j)}{mu(d)} ]

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

[sum_{d=1}mu(d)sum_{i=1}^{leftlfloorfrac{n}{k} ight floor}[dmid i]sum_{j=1}^{leftlfloorfrac{m}{k} ight floor}[dmid j] ]

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

易知(1sim leftlfloordfrac{n}{k} ight floor) 中 (d) 的倍数有 (leftlfloordfrac{leftlfloordfrac{n}{k} ight floor}{d} ight floor = leftlfloordfrac{n}{kd} ight floor) 个（由引理 1 可知），原式变为

[sum_{d=1}mu(d)leftlfloordfrac{n}{kd} ight floorleftlfloordfrac{m}{kd} ight floor ]

预处理 (mu) 后，显然可以数论分块求解，复杂度(O(n + Tsqrt{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
", Solve());
    return 0;
}

---

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

给定 (n,m) , 求：

[sum_{i=1}^nsum_{j=1}^{m} operatorname{lcm}(i,j)mod 20101009 ]

(1le n,mle 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}^nsum_{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}[dmid i]sum_{j=1}^m [dmid j][gcd(dfrac{i}{d}, dfrac{j}{d}=1)] dfrac{ij}{d} ]

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

[egin{aligned} &sum_{d=1}sum_{i=1}^{leftlfloorfrac{n}{d} ight floor}sum_{j=1}^{leftlfloorfrac{m}{d} ight floor}[gcd(i,j)=1]dfrac{id imes jd}{d}\ =& sum_{d=1}sum_{i=1}^{leftlfloorfrac{n}{d} ight floor}sum_{j=1}^{leftlfloorfrac{m}{d} ight floor}[gcd(i,j)=1]ijd\ =& sum_{d=1}d sum_{i=1}^{leftlfloorfrac{n}{d} ight floor}sum_{j=1}^{leftlfloorfrac{m}{d} ight floor}[gcd(i,j)=1]ij end{aligned}]

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

[egin{aligned} f(x,y) =& sumlimits_{i=1}^{x} sumlimits_{j=1}^{y}[gcd(i,j)=1]ij\ =& sum_{i=1}^{x} sum_{j=1}^{y}sum_{dmid gcd(i,j)}mu(d) ij\ =& sum_{d=1}mu(d)sum_{i=1}^{x}[dmid i] sum_{j=1}^{y}[dmid j] ij\ =& sum_{d=1}mu(d) d^2sum_{i=1}^{leftlfloor frac{x}{d} ight floor}sum_{j=1}^{leftlfloorfrac{y}{d} ight floor}ij end{aligned}]

前半部分 (sumlimits_{d=1}mu(d) d^2)，可以考虑筛出 (mu(d)) 后求前缀和。
后半部分是等差数列乘等差数列的形式，设 (g(p,q) = sumlimits_{i=1}^{p} sumlimits_{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) imes p}{2} imes dfrac{(1+q) imes q}{2} ]

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

[f(x,y) = sum_{d=1}mu(d) d^2cdot g(leftlfloor frac{x}{d} ight floor, leftlfloorfrac{y}{d} ight floor) ]

数论分块求解即可。

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

[sum_{d=1}dcdot f(leftlfloorfrac{n}{d} ight floor, leftlfloorfrac{m}{d} ight floor) ]

又是一个可以数论分块求解的形式。
线性筛预处理后 数论分块套数论分块，复杂度 (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<Tle 3 imes 10^5)，(1le nle 10^6)。

法一：无脑暴力

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

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

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

[egin{aligned} &nsum_{i=1}sum_{d|i} [gcd(i,n) = d]dfrac{i}{d}\ =& nsum_{i=1}sum_{d|i} [gcd(dfrac{i}{d},dfrac{n}{d}) = 1]dfrac{i}{d}\ =& nsum_{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})，消去整除的条件，原式等价于：

[nsum_{d|n}sum_{i=1}^{leftlfloorfrac{n}{d} ight floor}[gcd(i,dfrac{n}{d}) = 1]i ]

代入 ([gcd(i,j) = 1] = sumlimits_{dmid gcd(i,j)} {mu (d)})，原式等价于：

[nsum_{d|n}sum_{i=1}^{leftlfloorfrac{n}{d} ight floor} i sum_{t|gcd(i,frac{n}{d})}mu (t) ]

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

[nsum_{d|n}sum_{t|frac{n}{d}} mu(t) sum_{i=1}^{leftlfloorfrac{n}{d} ight floor} [t|i] i ]

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

[nsum_{d|n}sum_{t|frac{n}{d}} mu(t)t sum_{i=1}^{leftlfloorfrac{n}{dt} ight floor} i ]

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

[nsum_{d|n}sum_{t|frac{n}{d}} mu(t)tcdot g(leftlfloordfrac{n}{dt} ight floor) ]

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

[nsum_{T=1}^{n}g(leftlfloordfrac{n}{T} ight floor)sum_{t|T}mu(t)t ]

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

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

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

此时已经可以线性筛 + 数论分块求解，复杂度 (O(n+Tsqrt{n}))，比较菜鸡，时限 500ms 过不了。
考虑筛出 (f) 后再用埃氏筛预处理 (sumlimits_{T=1}^{n}g(leftlfloordfrac{n}{T} ight floor)f(T))，输出时乘上 (n)，复杂度变为 (O(nlog^2 n + n))。

---

法二：

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

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

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

[egin{aligned} &nsum_{d|n}sum_{i=1}^{leftlfloorfrac{n}{d} ight floor}[gcd(id,n) = d]i\ =& nsum_{d|n}sum_{i=1}^{leftlfloorfrac{n}{d} ight floor}[gcd(i,dfrac{n}{d}) = 1]i end{aligned}]

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

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

考虑 (sumlimits_{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}) 对这样的数。
则原式等价于：

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

可以直接预处理答案。
预处理时先线性筛出 (varphi)，再埃氏筛枚举 (i) 的倍数，令它们的答案加上 (frac{varphi(i)i}{2})，最后输出时乘上 (n)。
复杂度 (O(nlog^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
", 1ll * ans[n] * n);
  }
  return 0; 
}

---

[SDOI2015]约数个数和

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

[sum_{i=1}^{n}sum_{j=1}^moperatorname{d}(ij) ]

(1<T,nle 5 imes 10^4)。

一个结论：

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

证明

先考虑 (i = p^a)，(j=p^b(pin 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(pin mathrm{primes})) 时等式成立。

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

---

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

[egin{aligned} operatorname{d}(ij) =& sum_{x|i}sum_{y|j}[gcd(x,y)=1]\ =& sum_{x|i}sum_{y|j}sumlimits_{dmid 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) 提出来，消去整除的条件，原式等价于：

[egin{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}^{leftlfloorfrac{n}{d} ight floor}operatorname{d}(i)sum_{j=1}^{leftlfloorfrac{m}{d} ight floor}operatorname{d}(j) ]

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

[sum_{d=1}{mu (d)}Sleft({leftlfloorfrac{n}{d} ight floor} ight)Sleft({leftlfloorfrac{m}{d} ight floor} ight) ]

线性筛预处理 (mu,operatorname{d})，数论分块求解即可，复杂度 (O(n+Tsqrt{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
", ans);
  }
  return 0;
}
/*
in
1
32 43
*/
/*
out
15420
*/

---

P3768 简单的数学题

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

[left(sum_{i=1}^nsum_{j=1}^nicdot jcdot gcd(i,j) ight)mod p ]

(nleq10^{10})，(5 imes10^8leq pleq1.1 imes10^9)，(pin mathrm{primes})。
时限 4s。

无脑套路暴力。

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

[egin{aligned} &sum_{d=1}dsum_{i=1}^{n}sum_{j=1}^{n}[gcd(i,j)=d]ij\ =& sum_{d=1}dsum_{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})，消去整除的条件，原式等价于：

[egin{aligned} &sum_{d=1}dsum_{i=1}^{leftlfloorfrac{n}{d} ight floor}sum_{j=1}^{leftlfloorfrac{n}{d} ight floor}[gcd(i,j)=1]idcdot jd\ =& sum_{d=1} d^3sum_{i=1}^{leftlfloorfrac{n}{d} ight floor}sum_{j=1}^{leftlfloorfrac{n}{d} ight floor}[gcd(i,j)=1]ij end{aligned}]

代入 ([gcd(i,j) = 1] = sumlimits_{dmid gcd(i,j)} {mu (d)})，原式等价于：

[sum_{d=1} d^3sum_{i=1}^{leftlfloorfrac{n}{d} ight floor}sum_{j=1}^{leftlfloorfrac{n}{d} ight floor}ijsum_{t|gcd(i,j)}mu (t) ]

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

[sum_{d=1} d^3sum_{t=1}mu(t)sum_{i=1}^{leftlfloorfrac{n}{d} ight floor}[t|i] sum_{j=1}^{leftlfloorfrac{n}{d} ight floor}[t|j]ij ]

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

[sum_{d=1} d^3sum_{t=1}mu(t)t^2sum_{i=1}^{leftlfloorfrac{n}{dt} ight floor}sum_{j=1}^{leftlfloorfrac{n}{dt} ight floor}ij ]

发现后面两个求和是等差数列乘等差数列的形式。
设 (g(p,q) = sumlimits_{i=1}^{p} sumlimits_{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) imes p}{2} imes dfrac{(1+q) imes q}{2} ]

代入原式，原式等价于：

[sum_{d=1} d^3sum_{t=1}mu(t)t^2cdot gleft({leftlfloorfrac{n}{dt} ight floor},{leftlfloorfrac{n}{dt} ight floor} ight) ]

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

[egin{aligned} &sum_{T=1}^{n}gleft({leftlfloorfrac{n}{T} ight floor},{leftlfloorfrac{n}{T} ight floor} ight)sum_{d|T}d^3mu{left(dfrac{T}{d} ight)}left(dfrac{T}{d} ight)^2\ =& sum_{T=1}^{n}T^2 gleft({leftlfloorfrac{n}{T} ight floor},{leftlfloorfrac{n}{T} ight floor} ight)sum_{d|T}dcdot mu{left(dfrac{T}{d} ight)} end{aligned}]

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

[sum_{T=1}^{n}T^2 varphi(T) cdot gleft({leftlfloorfrac{n}{T} ight floor},{leftlfloorfrac{n}{T} ight floor} ight) ]

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

有：

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

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

[h(n) = sum_{d|n} d^2varphi(d)cdot gleft(dfrac{n}{d} ight) ]

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

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

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

[h(n) = n^3 ]

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

[egin{aligned} g(1)S(n) &= sum_{i=1}^{n}h(i) - sum_{i=2}^{n}g(i)Sleft(leftlfloordfrac{n}{i} ight floor ight)\ S(n) &= sum_{i=1}^{n} i^3 - sum_{i=2}^{n} i^2cdot Sleft(leftlfloordfrac{n}{i} ight floor ight) end{aligned}]

前一项是自然数的立方和，有 (sumlimits_{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} ]

(1le Tle 10^3)，(1le n,mle 10^6)。
5S，256MB。

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

[largeegin{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(sumlimits_{i=1}^{n}sumlimits_{j=1}^{m}[gcd(i,j)=d] ight)} end{aligned}]

分母套路一波，有：

[egin{aligned} &sumlimits_{i=1}^{n}sumlimits_{j=1}^{m}[gcd(i,j) = d]\ =& sumlimits_{i=1}^{leftlfloorfrac{n}{d} ight floor}sumlimits_{j=1}^{leftlfloorfrac{m}{d} ight floor}[gcd(i,j) = 1]\ =& sumlimits_{i=1}^{leftlfloorfrac{n}{d} ight floor}sumlimits_{j=1}^{leftlfloorfrac{m}{d} ight floor}sum_{kmid gcd(i,j)}mu (k)\ =& sum_{k=1}mu(k)leftlfloordfrac{n}{kd} ight floorleftlfloordfrac{m}{kd} ight floor end{aligned}]

代回原式，原式等于：

[large egin{aligned} &prod_{d = 1}^{n}f_{d}^{left(sumlimits_{k=1}mu(k)leftlfloorfrac{n}{kd} ight floorleftlfloorfrac{m}{kd} ight floor ight)}\ =& prod_{d = 1}^{n}left(f_{d}^{{sumlimits_{k=1}mu(k)}} ight)^{leftlfloorfrac{n}{kd} ight floorleftlfloorfrac{m}{kd} ight floor} end{aligned}]

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

[large egin{aligned} &prod_{d = 1}^{n}left(f_{d}^{{sumlimits_{k=1}mu(k)}} ight)^{leftlfloorfrac{n}{kd} ight floorleftlfloorfrac{m}{kd} ight floor}\ =& prod_{d = 1}^{n}left(prod_{k=1}^{leftlfloorfrac{n}{d} ight floor}f_{d}^{mu(k)leftlfloorfrac{n}{kd} ight floorleftlfloorfrac{m}{kd} ight floor} ight) end{aligned}]

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

[largeprod_{t=1}^{n}left(prod_{d | t}^{n}f_{d}^{{mu(frac{t}{d})}} ight)^{leftlfloorfrac{n}{t} ight floorleftlfloorfrac{m}{t} ight floor} ]

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

设：

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

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

[largeprod_{t=1}^{n}g(t)^{leftlfloorfrac{n}{t} ight floorleftlfloorfrac{m}{t} ight floor} ]

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

总时间复杂度 (O(nlog ^2 n + Tsqrt 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
", 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)} ight)^{f(type)} ]

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

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

(1le A,B,Cle 10^5)，(10^7le ple 1.05 imes 10^9)，(pin mathbb{P})，(T=70)。
2.5S，128MB。

先化下原式，原式等于：

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

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

[egin{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) imes f_1(b,a,c)}{f_2(a,b,c) imes f_2(a,c,b)} ]

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

---

type = 0

[egin{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 ight)^{b imes c} ]

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

---

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

[large egin{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) ight)^c\ =& left(prod_{d=1} d^{left(sumlimits_{i=1}^{a}sumlimits_{j=1}^{b}[gcd(i,j) = d] ight)} ight)^c end{aligned}]

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

[egin{aligned} &sumlimits_{i=1}^{a}sumlimits_{j=1}^{b}[gcd(i,j) = d]\ =& sumlimits_{i=1}^{leftlfloorfrac{a}{d} ight floor}sumlimits_{j=1}^{leftlfloorfrac{b}{d} ight floor}[gcd(i,j) = 1]\ =& sumlimits_{i=1}^{leftlfloorfrac{a}{d} ight floor}sumlimits_{j=1}^{leftlfloorfrac{b}{d} ight floor}sum_{kmid gcd(i,j)}mu (k)\ =& sum_{k=1}mu(k)leftlfloordfrac{a}{kd} ight floorleftlfloordfrac{b}{kd} ight floor end{aligned}]

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

[large egin{aligned} &left(prod_{d=1} d^{left(sumlimits_{k=1}mu(k)leftlfloorfrac{a}{kd} ight floorleftlfloorfrac{b}{kd} ight floor ight)} ight)^c\ =& left(prod_{d=1} left(d^{sumlimits_{k=1}mu(k)} ight)^{leftlfloorfrac{a}{kd} ight floorleftlfloorfrac{b}{kd} ight floor} ight)^c\ =& prod_{d=1} left(prod_{k=1}^{leftlfloorfrac{n}{d} ight floor}d^{left(mu(k)leftlfloorfrac{a}{kd} ight floorleftlfloorfrac{b}{kd} ight floor ight)} ight)^c end{aligned}]

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

[large prod_{t=1}^{n}left(left(prod_{d|t} d^{mu{left(frac{t}{d} ight)}} ight)^{leftlfloorfrac{a}{t} ight floorleftlfloorfrac{b}{t} ight floor} ight)^c ]

设：

[large g_0(t) = prod_{d|t}d^{muleft(frac{t}{d} ight)} ]

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

[large prod_{t=1}^{a}left(g_0(t)^{leftlfloorfrac{a}{t} ight floorleftlfloorfrac{b}{t} ight floor} ight)^c ]

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

---

type = 1

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

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

[large prod_{i=1}^{a}prod_{j=1}^{b}prod_{k=1}^{c} i^{i imes j imes k} = prod_{i=1}^{a}prod_{j=1}^{b}i^{left(i imes j imes sumlimits_{k = 1}^{c} k ight)} ]

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

[left(prod_{i=1}^{a}i^i ight)^{operatorname{sum}(b) imes operatorname{sum}(c)} ]

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

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

---

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

[large egin{aligned} &prod_{i=1}^{a}prod_{j=1}^{b}prod_{k=1}^{c} gcd(i,j)^{i imes j imes k}\ =& left(prod_{i=1}^aprod_{j=1}^bgcd(i,j)^{i imes j} ight)^{operatorname{sum}(c)} end{aligned}]

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

[large left(prod_{d=1}d^{left(sumlimits_{i=1}^a sumlimits_{j=1}^b i imes j[gcd(i,j)=d] ight)} ight)^{operatorname{sum}(c)} ]

大力化简指数，有：

[large egin{aligned} &sumlimits_{i=1}^a sumlimits_{j=1}^b i imes j[gcd(i,j)=d]\ =& d^2 sumlimits_{i=1}^{leftlfloorfrac{a}{d} ight floor} sumlimits_{j=1}^{leftlfloorfrac{b}{d} ight floor} i imes j[gcd(i,j)=1\ =& d^2 sumlimits_{i=1}^{leftlfloorfrac{a}{d} ight floor} i sumlimits_{j=1}^{leftlfloorfrac{b}{d} ight floor} jsumlimits_{t|gcd(i,j)}mu(t)\ =& d^2 sumlimits_{i=1}^{leftlfloorfrac{a}{d} ight floor} i sumlimits_{j=1}^{leftlfloorfrac{b}{d} ight floor} jsumlimits_{k|gcd(i,j)}mu(k)\ =& d^2 sumlimits_{k=1}mu(k)sumlimits_{i=1}^{leftlfloorfrac{a}{d} ight floor} i[k|i] sumlimits_{j=1}^{leftlfloorfrac{b}{d} ight floor} j[k|j]\ =& d^2 sumlimits_{k=1}k^2mu(k)sumlimits_{i=1}^{leftlfloorfrac{a}{kd} ight floor} isumlimits_{j=1}^{leftlfloorfrac{b}{kd} ight floor} j\ =& d^2sumlimits_{k=1}k^2mu(k)operatorname{sum}{left(leftlfloorfrac{a}{kd} ight floor ight)} operatorname{sum}{left(leftlfloorfrac{b}{kd} ight floor ight)}\ end{aligned}]

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

[large left(prod_{d=1}d^{left(d^2sumlimits_{k=1}k^2mu(k)operatorname{sum}{left(leftlfloorfrac{a}{kd} ight floor ight)} operatorname{sum}{left(leftlfloorfrac{b}{kd} ight floor ight)} ight)} ight)^{operatorname{sum}(c)} ]

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

[large egin{aligned} &left(prod_{d=1}left(prod_{k=1}^{leftlfloorfrac{n}{d} ight floor}d^{left(d^2 k^2mu(k) ight)} ight)^{left(operatorname{sum}{left(leftlfloorfrac{a}{kd} ight floor ight)} operatorname{sum}{left(leftlfloorfrac{b}{kd} ight floor ight)} ight)} ight)^{operatorname{sum}(c)}\ =& prod_{t=1}left(left(prod_{d|t}d^{left(d^2left(frac{t}{d} ight)^2muleft(frac{t}{d} ight) ight)} ight)^{operatorname{sum}{left(leftlfloorfrac{a}{t} ight floor ight)} operatorname{sum}{left(leftlfloorfrac{b}{t} ight floor ight)}} ight)^{operatorname{sum}(c)}\ =& prod_{t=1}left(left(prod_{d|t}d^{left(t^2muleft(frac{t}{d} ight) ight)} ight)^{operatorname{sum}{left(leftlfloorfrac{a}{t} ight floor ight)} operatorname{sum}{left(leftlfloorfrac{b}{t} ight floor ight)}} ight)^{operatorname{sum}(c)} end{aligned}]

与二面相同，设：

[large g_1(t) = prod_{d|t}d^{left(t^2muleft(frac{t}{d} ight) ight)} ]

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

[prod_{t=1}left(g_1(t)^{operatorname{sum}{left(leftlfloorfrac{a}{t} ight floor ight)} operatorname{sum}{left(leftlfloorfrac{b}{t} ight floor ight)}} ight)^{operatorname{sum}(c)} ]

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

注意指数的取模。

---

type = 2

[egin{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 egin{aligned} &prod_{i=1}^{a}prod_{j=1}^{b}prod_{k=1}^{c} i^{gcd(i,j,k)}\ =&prod_{d=1}prodlimits_{i=1}^{a}i^{left(sumlimits_{j=1}^{b}sumlimits_{k=1}^{c}[gcd(i,j,k)=d] ight)}\ =& prod_{d=1}prodlimits_{i=1}^{a}i^{left(sumlimits_{j=1}^{leftlfloorfrac{b}{d} ight floor}sumlimits_{k=1}^{leftlfloorfrac{c}{d} ight floor}[gcd(frac{i}{d},j,k)=1] ight)}\ =& prod_{d=1}prodlimits_{i=1}^{leftlfloorfrac{a}{d} ight floor}(id)^{left(dsumlimits_{j=1}^{leftlfloorfrac{b}{d} ight floor}sumlimits_{k=1}^{leftlfloorfrac{c}{d} ight floor}[gcd(i,j,k)=1] ight)}\ =& prod_{d=1}prodlimits_{i=1}^{leftlfloorfrac{a}{d} ight floor}(id)^{left(dsumlimits_{j=1}^{leftlfloorfrac{b}{d} ight floor}sumlimits_{k=1}^{leftlfloorfrac{c}{d} ight floor}sumlimits_{x|gcd(i,j,k)}{mu(x)} ight)}\ =& prod_{d=1}prodlimits_{i=1}^{leftlfloorfrac{a}{d} ight floor}(id)^{left(dsumlimits_{x=1}mu(x)[x|i]sumlimits_{j=1}^{leftlfloorfrac{b}{d} ight floor}[x|j]sumlimits_{k=1}^{leftlfloorfrac{c}{d} ight floor}[x|k] ight)}\ =& prod_{d=1}prod_{x=1}prodlimits_{i=1}^{leftlfloorfrac{a}{d} ight floor}(id)^{left(d imes mu(x)[x|i]sumlimits_{j=1}^{leftlfloorfrac{b}{d} ight floor}[x|j]sumlimits_{k=1}^{leftlfloorfrac{c}{d} ight floor}[x|k] ight)}\ =& prod_{d=1}prod_{x=1}prodlimits_{i=1}^{leftlfloorfrac{a}{xd} ight floor}(ixd)^{left(d imes mu(x){leftlfloorfrac{b}{xd} ight floor}{leftlfloorfrac{c}{xd} ight floor} ight)}\ =& prod_{t = 1}prod_{d|T}prod_{i=1}^{leftlfloorfrac{a}{t} ight floor}(it)^{left(d imes muleft(frac{t}{d} ight){leftlfloorfrac{b}{t} ight floor}{leftlfloorfrac{c}{t} ight floor} ight)}\ =& prod_{t = 1}prod_{d|T}left(t^{leftlfloorfrac{a}{t} ight floor}prod_{i=1}^{leftlfloorfrac{a}{t} ight floor}i ight)^{d imes muleft(frac{t}{d} ight){leftlfloorfrac{b}{t} ight floor}{leftlfloorfrac{c}{t} ight floor}}\ end{aligned}]

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

[large egin{aligned} &prod_{t = 1}prod_{d|t}left(t^{leftlfloorfrac{a}{t} ight floor} imes operatorname{fac}left(leftlfloorfrac{a}{t} ight floor ight) ight)^{left(d imes muleft(frac{t}{d} ight){leftlfloorfrac{b}{t} ight floor}{leftlfloorfrac{c}{t} ight floor} ight)}\ =& prod_{t = 1}left(prod_{d|t}left(t^{leftlfloorfrac{a}{t} ight floor} imes operatorname{fac}left(leftlfloorfrac{a}{t} ight floor ight) ight)^{d imes muleft(frac{t}{d} ight)} ight)^{{leftlfloorfrac{b}{t} ight floor}{leftlfloorfrac{c}{t} ight floor}}\ =& prod_{t = 1}left(left(t^{leftlfloorfrac{a}{t} ight floor} imes operatorname{fac}left(leftlfloorfrac{a}{t} ight floor ight) ight)^{sumlimits_{d|t}d imes muleft(frac{t}{d} ight)} ight)^{{leftlfloorfrac{b}{t} ight floor}{leftlfloorfrac{c}{t} ight floor}} end{aligned}]

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

[large egin{aligned} &prod_{t = 1}left(t^{leftlfloorfrac{a}{t} ight floor} imes operatorname{fac}left(leftlfloorfrac{a}{t} ight floor ight) ight)^{varphi(t){leftlfloorfrac{b}{t} ight floor}{leftlfloorfrac{c}{t} ight floor}}\ =& prod_{t = 1}left(t^{varphi(t)leftlfloorfrac{a}{t} ight floor{leftlfloorfrac{b}{t} ight floor}{leftlfloorfrac{c}{t} ight floor}} imes operatorname{fac}left(leftlfloorfrac{a}{t} ight floor ight)^{varphi(t){leftlfloorfrac{b}{t} ight floor}{leftlfloorfrac{c}{t} ight floor}} ight) end{aligned}]

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

---

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

[large egin{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(sumlimits_{i=1}^{a}sumlimits_{j=1}^{b} [gcd(i,j)=d] ight)} ight)^{sumlimits_{k=1}^{c}gcd(d,k)} end{aligned}]

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

[large egin{aligned} &sumlimits_{k=1}^{c}gcd(d,k)\ =& sumlimits_{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)leftlfloorfrac{c}{x} ight floor end{aligned}]

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

[sumlimits_{i=1}^{a}sumlimits_{j=1}^{b}[gcd(i,j) = d]=sum_{y=1}mu(y)leftlfloordfrac{a}{yd} ight floorleftlfloordfrac{b}{yd} ight floor ]

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

[large egin{aligned} &prod_{d=1} left(d^{left(sumlimits_{y=1}mu(y)leftlfloorfrac{a}{yd} ight floorleftlfloorfrac{b}{yd} ight floor ight)} ight)^{sumlimits_{x|d}varphi(x)leftlfloorfrac{c}{x} ight floor}\ =& prod_{d=1} left(prodlimits_{y=1}d^{left(mu(y)leftlfloorfrac{a}{yd} ight floorleftlfloorfrac{b}{yd} ight floor ight)} ight)^{sumlimits_{x|d}varphi(x)leftlfloorfrac{c}{x} ight floor} end{aligned}]

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

[large egin{aligned} &prod_{d=1} left(prodlimits_{y=1}d^{left(mu(y)leftlfloorfrac{a}{yd} ight floorleftlfloorfrac{b}{yd} ight floor ight)} ight)^{sumlimits_{x|d}varphi(x)leftlfloorfrac{c}{x} ight floor}\ =& prod_{t=1}prod_{d|t} d^{left(muleft(frac{t}{d} ight)leftlfloorfrac{a}{t} ight floorleftlfloorfrac{b}{t} ight floorsumlimits_{x|d}varphi(x)leftlfloorfrac{c}{x} ight floor ight)}\ =& prod_{t=1}left(prod_{d|t} d^{left(muleft(frac{t}{d} ight)sumlimits_{x|d}varphi(x)leftlfloorfrac{c}{x} ight floor ight)} ight)^{leftlfloorfrac{a}{t} ight floorleftlfloorfrac{b}{t} ight floor}\ =& prod_{t=1}left(prod_{d|t} prod_{x|d}d^{left(muleft(frac{t}{d} ight)varphi(x)leftlfloorfrac{c}{x} ight floor ight)} ight)^{leftlfloorfrac{a}{t} ight floorleftlfloorfrac{b}{t} ight floor} end{aligned}]

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

[large egin{aligned} &prod_{t=1}left(prod_{d|t} prod_{x|d}x^{left(muleft(frac{t}{d} ight)varphi(x)leftlfloorfrac{c}{x} ight floor ight)} ight)^{leftlfloorfrac{a}{t} ight floorleftlfloorfrac{b}{t} ight floor}\ &prod_{t=1}left(prod_{d|t} prod_{x|d}{left(frac{d}{x} ight)}^{left(muleft(frac{t}{d} ight)varphi(x)leftlfloorfrac{c}{x} ight floor ight)} ight)^{leftlfloorfrac{a}{t} ight floorleftlfloorfrac{b}{t} ight floor} end{aligned}]

---

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

[large egin{aligned} &prod_{x=1} prod_{t=1}^{operatorname{lim}}[x|t]left(prod_{d|t} [x|d]{x}^{left(muleft(frac{t}{d} ight)varphi(x)leftlfloorfrac{c}{x} ight floor ight)} ight)^{leftlfloorfrac{a}{t} ight floorleftlfloorfrac{b}{t} ight floor}\ =& prod_{x=1} prod_{t=1}^{leftlfloorfrac{operatorname{lim}}{x} ight floor}left(prod_{d|t} {x}^{left(muleft(frac{tx}{dx} ight)varphi(x)leftlfloorfrac{c}{x} ight floor ight)} ight)^{leftlfloorfrac{a}{tx} ight floorleftlfloorfrac{b}{tx} ight floor}\ =& prod_{x=1} prod_{t=1}^{leftlfloorfrac{operatorname{lim}}{x} ight floor}prod_{d|t} {x}^{left(varphi(x)leftlfloorfrac{c}{x} ight floorleftlfloorfrac{a}{tx} ight floorleftlfloorfrac{b}{tx} ight floormuleft(frac{t}{d} ight) ight)} end{aligned}]

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

[large egin{aligned} &prod_{x=1} {x}^{left(sumlimits_{t=1}^{leftlfloorfrac{operatorname{lim}}{x} ight floor}sumlimits_{d|t}varphi(x)leftlfloorfrac{c}{x} ight floorleftlfloorfrac{a}{tx} ight floorleftlfloorfrac{b}{tx} ight floormuleft(frac{t}{d} ight) ight)}\ =& prod_{x=1} {x}^{left(varphi(x)leftlfloorfrac{c}{x} ight floorsumlimits_{t=1}^{leftlfloorfrac{operatorname{lim}}{x} ight floor}leftlfloorfrac{a}{tx} ight floorleftlfloorfrac{b}{tx} ight floorsumlimits_{d|t}muleft(frac{t}{d} ight) ight)} end{aligned}]

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

[large egin{aligned} &prod_{x=1} {x}^{left(varphi(x)leftlfloorfrac{c}{x} ight floorsumlimits_{t=1}^{leftlfloorfrac{operatorname{lim}}{x} ight floor}leftlfloorfrac{a}{tx} ight floorleftlfloorfrac{b}{tx} ight floor[t=1] ight)}\ =& prod_{x=1} {x}^{left(varphi(x)leftlfloorfrac{a}{x} ight floorleftlfloorfrac{b}{x} ight floorleftlfloorfrac{c}{x} ight floor ight)} end{aligned}]

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

---

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

[large egin{aligned} &prod_{x=1} prod_{t=1}^{operatorname{lim}}[x|t]left(prod_{d|t} [x|d]{left(frac{d}{x} ight)}^{left(muleft(frac{t}{d} ight)varphi(x)leftlfloorfrac{c}{x} ight floor ight)} ight)^{leftlfloorfrac{a}{t} ight floorleftlfloorfrac{b}{t} ight floor}\ =& prod_{x=1} prod_{t=1}^{leftlfloorfrac{operatorname{lim}}{x} ight floor}left(prod_{d|tx} [x|d]{left(frac{d}{x} ight)}^{left(muleft(frac{tx}{d} ight)varphi(x)leftlfloorfrac{c}{x} ight floor ight)} ight)^{leftlfloorfrac{a}{tx} ight floorleftlfloorfrac{b}{tx} ight floor}\ =& prod_{x=1} left(prod_{t=1}^{leftlfloorfrac{operatorname{lim}}{x} ight floor}left(prod_{d|tx} [x|d]{left(frac{d}{x} ight)}^{muleft(frac{tx}{d} ight)} ight)^{leftlfloorfrac{a}{tx} ight floorleftlfloorfrac{b}{tx} ight floor} ight)^{varphi(x)leftlfloorfrac{c}{x} ight floor} end{aligned}]

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

[large prod_{x=1} left(prod_{t=1}^{leftlfloorfrac{operatorname{lim}}{x} ight floor}left(prod_{d|t}d^{muleft(frac{t}{d} ight)} ight)^{leftlfloorfrac{a}{tx} ight floorleftlfloorfrac{b}{tx} ight floor} ight)^{varphi(x)leftlfloorfrac{c}{x} ight floor} ]

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

[large prod_{x=1} left(prod_{t=1}^{leftlfloorfrac{operatorname{lim}}{x} ight floor}g_0(t)^{leftlfloorfrac{a}{tx} ight floorleftlfloorfrac{b}{tx} ight floor} ight)^{varphi(x)leftlfloorfrac{c}{x} ight floor} ]

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

[forall a,b,cin mathbb{Z},leftlfloordfrac{a}{bc} ight floor = leftlfloor{dfrac{leftlfloordfrac{a}{b} ight floor}{c}} ight floor ]

则原式等于：

[large prod_{x=1} left(prod_{t=1}^{leftlfloorfrac{operatorname{lim}}{x} ight floor}g_0(t)^{leftlfloorfrac{leftlfloorfrac{a}{x} ight floor}{t} ight floorleftlfloorfrac{leftlfloorfrac{b}{x} ight floor}{t} ight floor} ight)^{varphi(x)leftlfloorfrac{c}{x} ight floor} ]

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

然后就做完了，哈哈哈。

---

一些实现上的小技巧：

• 逆元能预处理就预处理。
• 注意对指数取模，模数为 (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("
");
  }
  return 0;
}

---

写在最后

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

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