莫比乌斯反演

莫比乌斯函数

对于一个正整数$n$，$\mu (n)$为其莫比乌斯函数。

莫比乌斯函数是怎么样定义的呢？

首先，对于任意一个正整数$n$来说，由算数基本定理知，一定可以将其分解成$n=p_1^{a_1}{p}_2^{a_2}...p_k^{a_k}$的形式。

定义$$\mu (n)=\{\begin{array}{l}(-1{)}^k,a_1=a_2=...=a_k=1\\ 0,otherwise\end{array}$$

关于莫比乌斯函数和的一个定理：

$$\sum _{d|n}\mu (d)=\{\begin{array}{l}1,n=1\\ 0,n>1\end{array},(d是n的所有约数)$$

证明：

• 当$n=1$时，显然有$\mu (1)=1$。

• 当$n>1$时，将$n$分解成$n=p_1^{a_1}{p}_2^{a_2}...p_k^{a_k}$，那么$n$的约数$d$一定是这种形式$d=p_1^{b_1}{p}_2^{b_2}...p_k^{b_k},(0<=b_i<=a_i)$。
  由莫比乌斯函数的定义知道，当其中某一个质数的次方$b_i>1$，那么它的莫比乌斯函数$\mu =0$，对答案没有贡献，所以我们不考虑$b_i$大于$1$的情况，于是乎，$b_i$的取值要么是$0$，要么是$1$。
  当这$k$个$b_i$中有奇数个取$1$，那么它的$\mu =-1$，偶数个$1$，$\mu =1$。
  于是，$$\begin{gathered} \sum _{d|n}\mu (d)=C_k^0\mu (d_0)+C_k^1\mu (d_1)+...+C_k^k\mu (d_k) \\ =C_k^0-C_k^1+C_k^2-C_k^3+.... \end{gathered}$$
  由二项式定理，$(a+b{)}^k=C_k^0{a}^0{b}^k+C_k^1{a}^1{b}^{k-1}+C_k^2{a}^2{b}^{k-2}+...+C_k^k{a}^k{b}^0$。
  代入$a=1,b=-1$,有：
  $(1-1{)}^k=C_k^0-C_k^1+C_k^2-C_k^3+....=0$。
  于是得到，$$\sum _{d|n}\mu (d)=0$$

---

先写这么多吧。。。好难

---

莫比乌斯反演

莫比乌斯反演有两种形式。

• 第一种
  $$若F(n)=\sum _{d|n}f(d)，则f(n)=\sum _{d|n}\mu (d)F(\frac{n}{d}),{}_{(d是n的所有约数)}$$

证明：

$f(n)={\sum }_{d|n}\mu (d)F(\frac{n}{d})$

$={\sum }_{d|n}\mu (d){\sum }_{k|\frac{n}{d}}f(k)---{>}_{将F按定义替换，其中k是\frac{n}{d}的所有约数}$

$={\sum }_{k|n}f(k){\sum }_{d|\frac{n}{k}}\mu (d)---{>}_{交换\mu 和f的枚举顺序}$

---

这里为什么可以这样交换呢？

首先，我们之前是先枚举$d$，然后枚举$k$，满足$d|n$且$k|\frac{n}{d}$。
如果我们先枚举$k$，那么我们先要找到$k$所有可能的取值从而再去枚举，从上面的关系可以看出，$k$是$\frac{n}{d}$的所有约数，当$d=1$的时候，$k$此时的所有取值就是$k$的所有取值，因为所有$\frac{n}{d}$的约数一定都是$n$的约数。

所以，第一个$\sum$的条件就是$k|n$，来枚举$k$所有可能的取值，那么接下来就是所有能以该$k$配对的$d$，使得它们仍然满足最上面的两个条件:$d|n$且$k|\frac{n}{d}$，我们可以对第二个条件变形，有$k|\frac{n}{d}\to kd|n\to d|\frac{n}{k}$，所以，第二个$\sum$的条件就是$d|\frac{n}{k}$（用来枚举所有与当前$k$配对的$d$）。

---

继续我们的证明

$f(n)={\sum }_{k|n}f(k){\sum }_{d|\frac{n}{k}}\mu (d)$

当我们写到这里的时候，我们发现第二维的求和正好对应了上面我们证明的莫比乌斯函数和的定理。
该定理告诉了我们，只有当${\sum }_{d|n}\mu (d)$中的$n=1$的时候，这个和为$1$，其余情况全都等于$0$。

所以我们只看$\frac{n}{k}=1$，也就是$k=n$的情况，当$k=n$的时候有：

${\sum }_{k|n}f(k){\sum }_{d|\frac{n}{k}}\mu (d)$

$={\sum }_{1}f(n)\ast 1$

$=f(n)$。

而当$\frac{n}{k}\ne 1$，也就是大于$1$的时候，由于第二个求和的结果是$0$，整个式子乘上$0$结果也是$0$。

所以证明出来了${\sum }_{k|n}f(k){\sum }_{d|\frac{n}{k}}\mu (d)=f(n)+0=f(n)$

至此，莫比乌斯反演的第一种形式得证。

---

再来看莫比乌斯反演的第二种形式

• 第二种

$$若F(n)=\sum _{n|d}f(d),则f(n)=\sum _{n|d}\mu (\frac{d}{n})F(d),{}_{(d是n的所有倍数)}$$

证明：

$f(n)={\sum }_{n|d}\mu (\frac{d}{n})F(d)$

先设$k=\frac{d}{n}$，则$d=nk$，那么式子可以写成

$={\sum }_{n|nk}\mu (k)F(nk)$

$={\sum }_{k=1}^{+\infty }\mu (k)F(nk)---{>}_{化简一下，约掉n，因为是枚举的n的所有倍数，是有无穷多的，k上限是+\infty }$

$={\sum }_{k=1}^{+\infty }\mu (k){\sum }_{nk|t}f(t)---{>}_{将F按定义替换，其中t是nk的所有倍数}$

$={\sum }_{n|t}f(t){\sum }_{k|\frac{t}{n}}\mu (k)---{>}_{交换\mu 和f的枚举顺序}$

---

这里的交换方法和上面的类似。
原本是先枚举$k$，在枚举$t$，且满足$nk|t$。
同样的，我们如果要先枚举$t$，在枚举$k$，就要先找到$t$所有可能的取值，然后就是该取值下所有和t配对的$k$。

很明显，$t$的所有取值就是$n$的所有倍数，所以第一个$\sum$的条件就是$n|t$，然后找到与该$t$配对的所有的$k$，使得$t$和$k$满足，$nk|t$，那么就有$nk|t\to k|\frac{t}{n}$，于是乎，第二个$\sum$的条件就是$k|\frac{t}{n}$。

---

然后，有了
$f(n)={\sum }_{n|t}f(t){\sum }_{k|\frac{t}{n}}\mu (k)$

同样的套路，只有当$\frac{t}{n}=1$的时候，才有值，其他的都是$0$。

那么代入$t=n$，有：

${\sum }_{n|t}f(t){\sum }_{k|\frac{t}{n}}\mu (k)$

$={\sum }_{1}f(n)\ast 1$

$=f(n)$

证毕。

---

感觉学的晕乎乎的，证明写的也不知道对不对

---

模版题目

例题1

解法：

转化一下题目：
将每一个点对$(x,y)$对应到坐标轴上，其实就是求在以点$(a,c)$为左下角$(b,d)$为右上角的区域内，满足$gcd(x,y)=k$的点对的个数。

直接求该区域内的点对个数不好求，我们可以利用类似二维前缀和的思想，定义$S(a,b)$表示以点$(1,1)$为左下角，点$(a,b)$为右上角的区域内点对的个数。那么上图中的区域就可以表示为$S(b,d)-S(a-1,d)-S(b,c-1)+S(a-1,c-1)$。

那么问题就转换成：对于一个以点$(1,1)$为左下角，点$(a,b)$为右上角的区域，我们如何快速的求出该区域内满足$gcd(x,y)=k$的所有点对的个数呢？

设答案个数为$ans$。
那么有$$ans=\sum _{i=1}^a\sum _{j=1}^b[gcd(i,j)=k]{,}_{(如果中括号里的条件成立就是1，否则是0)}$$

接下来在套用莫比乌斯反演公式之前，我们首先要定义出来$F(n)$和$f(n)$分别表示什么。
（我也不知道是怎么想出来这个定义的，还是记住吧QAQ）

定义
$$\begin{gathered} F(n)=\sum _{i=1}^a\sum _{j=1}^b[n|gcd(i,j)] \\ 表示满足gcd(i,j)是n的倍数的点对个数 \end{gathered}$$

定义
$$\begin{gathered} f(n)=\sum _{i=1}^a\sum _{j=1}^b[gcd(i,j)=n] \\ 表示gcd(i,j)恰好等于n的点对的个数 \end{gathered}$$

那么由定义可以推知$F$与$f$之间的关系为:
$$F(n)=\sum _{n|d}f(d)$$

正好满足莫比乌斯反演公式，反演一下有：
$$f(n)=\sum _{n|d}\mu (\frac{d}{n})F(d)$$

首先，我们要求解$f$，就要先把$F$的计算公式推出来。

$F(n)={\sum }_{i=1}^a{\sum }_{j=1}^b[n|gcd(i,j)]$

如果说$gcd(i,j)$是$n$的倍数，那么$i$肯定是$n$的倍数，$j$肯定也是$n$的倍数。
在$1$~$a$范围内，是$n$的倍数的数有$\lfloor \frac{a}{n}\rfloor$个，所以$i$的取值有$\lfloor \frac{a}{n}\rfloor$种，同理，$j$的取值有$\lfloor \frac{b}{n}\rfloor$种，由乘法原理得到，所有这样的$(i,j)$点对有$\lfloor \frac{a}{n}\rfloor \ast \lfloor \frac{b}{n}\rfloor$个。

于是$F(n)={\sum }_{i=1}^a{\sum }_{j=1}^b[n|gcd(i,j)]=\lfloor \frac{a}{n}\rfloor \ast \lfloor \frac{a}{n}\rfloor$

所以，

$$f(n)=\sum _{n|d}\mu (\frac{d}{n})F(d)=\sum _{n|d}\mu (\frac{d}{n})\ast \lfloor \frac{a}{d}\rfloor \ast \lfloor \frac{b}{d}\rfloor$$

再次整理一下，设$t=\frac{d}{n}$，则$d=nt$。

$f(n)={\sum }_{n|d}\mu (\frac{d}{n})\ast \lfloor \frac{a}{d}\rfloor \ast \lfloor \frac{b}{d}\rfloor$
$={\sum }_{t=1}^{+\infty }\mu (t)\ast \lfloor \frac{a}{nt}\rfloor \ast \lfloor \frac{b}{nt}\rfloor$

因为$a、b、n$都是确定的，我们用$a^{{}^{\prime }}=\frac{a}{n}$和$b^{{}^{\prime }}=\frac{b}{n}$来代替。

那么式子变为

$f(n)={\sum }_{t=1}^{+\infty }\mu (t)\ast \lfloor \frac{a^{{}^{\prime }}}{t}\rfloor \ast \lfloor \frac{b^{{}^{\prime }}}{t}\rfloor$

当$t>min(a^{{}^{\prime }},b^{{}^{\prime }})$的时候，式子的值就等于$0$了，所以我们没必要从$1$枚举到$+\infty$，上限到$min(a^{{}^{\prime }},b^{{}^{\prime }})$就够了。

于是乎，我们最终的式子就是!!!
$$f(n)=\sum _{t=1}^{min(a^{{}^{\prime }},b^{{}^{\prime }})}u(t)\ast \lfloor \frac{a^{{}^{\prime }}}{t}\rfloor \ast \lfloor \frac{b^{{}^{\prime }}}{t}\rfloor$$

在计算这个式子的时候，我们可以利用整除分块来降低计算的时间复杂度。

上代码：

#include <iostream>
#include <cstdio>
#include <cmath>
#include <cstring>
#include <set>
#include <map>
#include <queue>
#include <stack>
#include <vector>
#include <string>
#include <algorithm>
#define INF 0x3f3f3f3f

using namespace std;

typedef long long LL;

const int N = 50010;

int n, m;
int cnt, prime[N], mu[N];
bool vis[N];

void get_prime(int x) //利用线性筛法求莫比乌斯函数
{
    mu[1] = 1;
    for(int i = 2; i <= x; i ++)
    {
        if(!vis[i]) prime[cnt ++] = i, mu[i] = -1;
        for(int j = 0; prime[j] <= x / i; j ++)
        {
            vis[prime[j] * i] = true;
            if(i % prime[j] == 0)   break;
            mu[prime[j] * i] = - mu[i];
        }
    }

    for(int i = 1; i <= x; i ++)    mu[i] += mu[i - 1]; //对mu进行一个前缀和
}

LL S(int a, int b, int k) //计算左下角(1, 1)右上角(a, b)的区域内的点对个数
{
    a /= k; b /= k;
    int t = min(a, b);
    
    LL res = 0;
    for(int l = 1, r; l <= t; l = r + 1)
    {
        r = min(t, min(a / (a / l), b / (b / l)));
        res += (LL)(mu[r] - mu[l - 1]) * (a / l) * (b / l);
    }

    return res;
}

int main()
{
#ifdef LOCAL
    freopen("in.in", "r", stdin);
    freopen("out.out", "w", stdout);
#endif

    get_prime(N - 1); //预处理求出所有的莫比乌斯函数

    int T;
    scanf("%d", &T);

    while(T --)
    {
        int a, b, c, d, k;
        scanf("%d%d%d%d%d", &a, &b, &c, &d, &k);
        printf("%lld\n", S(b, d, k) - S(a - 1, d, k) - S(b, c - 1, k) + S(a - 1, c - 1, k));
    }

    return 0;
}

---

（ 趁热打铁，再来一道）
P1447 [NOI2010] 能量采集

（题面就不在这里贴了）

解法：

一共$n\ast m$个植物，每一个植物都有一个$k$值，然后每一个植物都会产生$2k+1$的能量浪费，求总共有多少能量被浪费。

因为这个$2k+1$中的$1$是固定下来的，所以对于n*m个植物，肯定会有$n\ast m$点浪费，先把这个固定的$1$产生的贡献$n\ast m$加上。于是，每一个植物的能量浪费就变成了$2k$，我们只需要统计所有的植物的$k$值的和,$K=k_1+k_2+...+k_{nm}$，然后在乘$2$，最后答案$ans=2K+n\ast m$。

对于每一个点，它的$k$值怎么求?
由题意应该很容易想到对于一个点$(i,j)$，它的$k=gcd(i,j)-1$。

于是，
$$K=\sum _{i=1}^n\sum _{j=1}^m(gcd(i,j)-1)=\sum _{i=1}^n\sum _{j=1}^{m}gcd(i,j)-n\ast m$$

由于$K$的特殊形式，答案ans也可以化简一下
$$ans=2\ast \sum _{i=1}^n\sum _{j=1}^{m}gcd(i,j)-n\ast m$$

问题变成了求${\sum }_{i=1}^n{\sum }_{j=1}^{m}gcd(i,j)$。

我们要把该式子变成类似于上一题的形式。

可以先枚举所有的$gcd(i,j)$的值d，然后统计有多少个$(i,j)$对的$gcd=d$，假设求出来$s$对，那么，这个约数$d$产生的贡献就是$d\ast s$。

于是，式子就变成了这样(假定$n<m$)
$$\sum _{i=1}^n\sum _{j=1}^{m}gcd(i,j)=\sum _{d=1}^{n}d\ast \sum _{i=1}^n\sum _{j=1}^m[gcd(i,j)=d]$$

然后把这个拿出来：${\sum }_{i=1}^n{\sum }_{j=1}^m[gcd(i,j)=d]$

类似上一题的定义方式(不能说是类似，应该说是一模一样 )：
定义
$$\begin{gathered} F(n)=\sum _{i=1}^a\sum _{j=1}^b[n|gcd(i,j)] \\ 表示满足gcd(i,j)是n的倍数的点对个数 \end{gathered}$$

定义
$$\begin{gathered} f(n)=\sum _{i=1}^a\sum _{j=1}^b[gcd(i,j)=n] \\ 表示gcd(i,j)恰好等于n的点对的个数 \end{gathered}$$

接下来就和上一题的推理方式一样了，我就不推了。

反正最后求出来$f(n)$的表达式是这样的。

$$\begin{gathered} f(n)=\sum _{t=1}^{min(a,b)}\mu (t)\ast \lfloor \frac{a}{t}\rfloor \ast \lfloor \frac{b}{t}\rfloor \\ 其中a=\frac{N}{d},b=\frac{M}{d} \\ N和M是题目中的范围，d是当前枚举的约数 \end{gathered}$$

废话少说，上代码!!!

#include <iostream>
#include <cstdio>
#include <cmath>
#include <cstring>
#include <set>
#include <map>
#include <queue>
#include <stack>
#include <vector>
#include <string>
#include <algorithm>
#define INF 0x3f3f3f3f

using namespace std;

typedef long long LL;

const int N = 1e5 + 10;

int n, m;
int cnt, prime[N], mu[N];
bool vis[N];

void get_prime(int x)
{
    mu[1] = 1;
    for(int i = 2; i <= x; i ++)
    {
        if(!vis[i]) prime[cnt ++] = i, mu[i] = -1;
        for(int j = 0; prime[j] <= x / i; j ++)
        {
            vis[prime[j] * i] = true;
            if(i % prime[j] == 0)   break;
            mu[prime[j] * i] = -mu[i];
        }
    }

    for(int i = 1; i <= x; i ++)    mu[i] += mu[i - 1];
}

LL solve(int n, int m, int d)
{
    n /= d; m /= d;
    int t = min(n, m);
    LL res = 0;
    for(int l = 1, r; l <= t; l = r + 1)
    {
        r = min(t, min(n / (n / l), m / (m / l)));
        res += (LL)(mu[r] - mu[l - 1]) * (n / l) * (m / l);
    }
    return res;
}

int main()
{
#ifdef LOCAL
    freopen("in.in", "r", stdin);
    freopen("out.out", "w", stdout);
#endif

    cin >> n >> m;

    int t = min(n, m);
    get_prime(t);

    LL res = 0;
    for(int d = 1; d <= t; d ++) //枚举所有gcd(i, j)的值
        res += d  * solve(n, m, d);

    res = res * 2 - (LL)n * m;

    cout << res << endl;

    return 0;
}

---

update-2021.10.8

进阶题目

P1829 [国家集训队]Crash的数字表格 / JZPTAB
题目大意：
求
$$\sum _{i=1}^N\sum _{j=1}^{M}lcm(i,j)(mod20101009)$$

做法：

• 对式子进行变形：

${\sum }_{i=1}^N{\sum }_{j=1}^{M}lcm(i,j)$

$={\sum }_{i=1}^N{\sum }_{j=1}^M\frac{ij}{gcd(i,j)}$

$={\sum }_{d=1}^T{\sum }_{i=1}^N{\sum }_{j=1}^M\frac{ij}{d}\ast [gcd(i,j)=d]{,}_{(其中T=min(N,M)}$

$={\sum }_{d=1}^T\frac{1}{d}{\sum }_{i=1}^N{\sum }_{j=1}^{M}ij\ast [gcd(i,j)=d]$

$={\sum }_{d=1}^T\frac{1}{d}{\sum }_{i=1}^{\frac{N}{d}}{\sum }_{j=1}^{\frac{M}{d}}ij\ast d^2\ast [gcd(i,j)=1]$

$={\sum }_{d=1}^{T}d{\sum }_{i=1}^{\frac{N}{d}}{\sum }_{j=1}^{\frac{M}{d}}ij\ast [gcd(i,j)=1]$

设函数$g(N,M)$表示${\sum }_{i=1}^N{\sum }_{j=1}^{M}ij\ast [gcd(i,j)=1]$

则上式等于： ${\sum }_{d=1}^{T}d\ast g(\frac{N}{d},\frac{M}{d})$

• 接下来的问题变为如何快速的求出$g(N,M)$.

$$g(N,M)=\sum _{i=1}^N\sum _{j=1}^{M}ij\ast [gcd(i,j)=1]$$

• 这里就可以使用莫比乌斯反演来求：

设$F(n)$表示${\sum }_{i=1}^N{\sum }_{j=1}^{M}ij\ast [n|gcd(i,j)]$。
设$f(n)$表示${\sum }_{i=1}^N{\sum }_{j=1}^{M}ij\ast [gcd(i,j)=n]$。

则有$F(n)={\sum }_{n|d}f(d)$。
反演一下有$f(n)={\sum }_{n|d}\mu (\frac{d}{n})F(d)$。

• 我们要求的是$f(1)$。

$f(1)={\sum }_{d=1}^{+\infty }\mu (d)F(d)$。

• 然后我们需要推导出$F(n)$的表达式。

$F(n)={\sum }_{i=1}^N{\sum }_{j=1}^{M}ij\ast [n|gcd(i,j)]$

$={\sum }_{i=1}^{\frac{N}{n}}{\sum }_{j=1}^{\frac{M}{n}}ij\ast n^2$

$=n^2{\sum }_{i=1}^{\frac{N}{n}}i{\sum }_{j=1}^{\frac{M}{n}}j$

$=n^2\ast \frac{\frac{N}{n}(\frac{N}{n}+1)}{2}\ast \frac{\frac{M}{n}(\frac{M}{n}+1)}{2}{,}_{(注意，此处的除法都是整除)}$

于是得到,$F(n)=n^2\ast \frac{\frac{N}{n}(\frac{N}{n}+1)}{2}\ast \frac{\frac{M}{n}(\frac{M}{n}+1)}{2}$

则$f(1)={\sum }_{d=1}^{+\infty }\mu (d)F(d)$
$={\sum }_{i=1}^{+\infty }\mu (i)F(i){,}_{(换个变量名表示)}$

$={\sum }_{i=1}^{+\infty }\mu (i)\ast i^2\ast \frac{\frac{N}{i}(\frac{N}{i}+1)}{2}\ast \frac{\frac{M}{i}(\frac{M}{i}+1)}{2}$

$={\sum }_{i=1}^{min(N,M)}\mu (i)\ast i^2\ast \frac{\frac{N}{i}(\frac{N}{i}+1)}{2}\ast \frac{\frac{M}{i}(\frac{M}{i}+1)}{2}$

• 最终，我们得到

${\sum }_{i=1}^N{\sum }_{j=1}^{M}lcm(i,j)={\sum }_{d=1}^{T}d\ast g(\frac{N}{d},\frac{M}{d})$

$g(N,M)=f(1)={\sum }_{i=1}^{min(N,M)}\mu (i)\ast i^2\ast \frac{\frac{N}{i}(\frac{N}{i}+1)}{2}\ast \frac{\frac{M}{i}(\frac{M}{i}+1)}{2}$

接下来就可以愉快的做题了！

#include <iostream>
#include <cstdio>
#include <cmath>
#include <cstring>
#include <set>
#include <map>
#include <queue>
#include <stack>
#include <vector>
#include <string>
#include <algorithm>
#define INF 0x3f3f3f3f

using namespace std;

typedef long long LL;

const int N = 1e7 + 10, MOD = 20101009;

int cnt, prime[N], mu[N];
LL sum[N]; //mu[i] * i * i的前缀和
bool vis[N];

void init(int x)
{
    mu[1] = 1;
    for(int i = 2; i <= x; i ++)
    {
        if(!vis[i]) prime[cnt ++] = i, mu[i] = -1;

        for(int j = 0; prime[j] <= x / i; j ++)
        {
            vis[prime[j] * i] = true;
            if(i % prime[j] == 0)   break;
            mu[prime[j] * i] = -mu[i];
        }
    }
    for(int i = 1; i <= x; i ++)
        sum[i] = (sum[i - 1] + (LL)i * i % MOD * (mu[i] + MOD) % MOD) % MOD;
}

LL cal(int x, int y)
{
    return (((LL)x * (x + 1) / 2) % MOD) * (((LL)y * (y + 1) / 2) % MOD) % MOD;
}

LL g(int n, int m)
{
    int t = min(n, m);
    LL res = 0;

    for(int l = 1, r; l <= t; l = r + 1)
    {
        r = min(t, min(n / (n / l), m / (m / l)));
        res = (res + (sum[r] - sum[l - 1] + MOD) % MOD * cal(n / l, m / l)% MOD) % MOD;
    }

    return res;
}

int main()
{
#ifdef LOCAL
    freopen("in.in", "r", stdin);
    freopen("out.out", "w", stdout);
#endif

    int n, m;
    cin >> n >> m;

    int t = min(n, m);
    init(t);

    LL ans = 0;
    for(int d = 1; d <= t; d ++)
        ans = (ans + (LL)d * g(n / d, m / d)) % MOD;
    
    cout << ans << endl;

    return 0;
}

---

AcWing 1358. 约数个数和

做法：

• 首先，$d(ij)={\sum }_{x|i}{\sum }_{y|j}[gcd(x,y)=1]$（根本想不到的转化QAQ）。

证明有点麻烦，不贴了。

所以，原式子等于：
$$\sum _{i=1}^N\sum _{j=1}^{M}d(ij)=\sum _{i=1}^N\sum _{j=1}^M\sum _{x|i}\sum _{y|j}[gcd(x,y)=1]$$

• 对式子进行化简：

${\sum }_{i=1}^N{\sum }_{j=1}^M{\sum }_{x|i}{\sum }_{y|j}[gcd(x,y)=1]$
$={\sum }_{x=1}^N{\sum }_{y=1}^M{\sum }_{x|i}{\sum }_{y|j}[gcd(x,y)=1]{,}_{(交换一下x，y，i，j的枚举顺序)}$

$={\sum }_{x=1}^N{\sum }_{y=1}^M[gcd(x,y)=1]{\sum }_{x|i}{\sum }_{y|j}1$

$={\sum }_{x=1}^N{\sum }_{y=1}^M[gcd(x,y)=1]\ast \lfloor \frac{N}{x}\rfloor \ast \lfloor \frac{M}{y}\rfloor$

$={\sum }_{i=1}^N{\sum }_{j=1}^M[gcd(i,j)=1]\ast \lfloor \frac{N}{i}\rfloor \ast \lfloor \frac{M}{j}\rfloor {,}_{(换一下变量名，x变成i，y变成j)}$

• 到这里后，就可以开始进行莫比乌斯反演了。

设$F(n)={\sum }_{i=1}^N{\sum }_{j=1}^M[n|gcd(i,j)]\ast \lfloor \frac{N}{i}\rfloor \ast \lfloor \frac{M}{j}\rfloor$。

设$f(n)={\sum }_{i=1}^N{\sum }_{j=1}^M[gcd(i,j)=n]\ast \lfloor \frac{N}{i}\rfloor \ast \lfloor \frac{M}{j}\rfloor$。

则有：$F(n)={\sum }_{n|d}f(d)$。

反演一下有：$f(n)={\sum }_{n|d}\mu (\frac{d}{n})F(d)$。

我们要求的是$f(1)={\sum }_{d=1}^{+\infty }\mu (d)F(d)$。

• 接下来推导出$F(n)$的表达式。

$$F(n)=\sum _{i=1}^N\sum _{j=1}^M[n|gcd(i,j)]\ast \lfloor \frac{N}{i}\rfloor \ast \lfloor \frac{M}{j}\rfloor$$

${\sum }_{i=1}^N{\sum }_{j=1}^M[n|gcd(i,j)]\ast \lfloor \frac{N}{i}\rfloor \ast \lfloor \frac{M}{j}\rfloor$

$={\sum }_{i=1}^{\frac{N}{n}}{\sum }_{j=1}^{\frac{M}{n}}\lfloor \frac{N}{ni}\rfloor \ast \lfloor \frac{M}{nj}\rfloor$

$={\sum }_{i=1}^{\frac{N}{n}}\lfloor \frac{N}{ni}\rfloor \ast {\sum }_{j=1}^{\frac{M}{n}}\lfloor \frac{M}{nj}\rfloor$

设$g(n)={\sum }_{i=1}^n\lfloor \frac{n}{i}\rfloor$。

于是,$F(n)=g(\frac{N}{n})\ast g(\frac{M}{n})$。

• 那么，$f(1)$的最终式子是：

$f(1)={\sum }_{i=1}^{+\infty }\mu (i)F(i)$

$={\sum }_{i=1}^{min(N,M)}\mu (i)\ast g(\frac{N}{i})\ast g(\frac{M}{i})$

${}_{(除法都是整除，所以上面的式子就是这个)}$

$={\sum }_{i=1}^{min(N,M)}\mu (i)\ast g(\lfloor \frac{N}{i}\rfloor )\ast g(\lfloor \frac{M}{i}\rfloor )$

${}_{(我们求这个式子的时候也是用整除分块来缩短时间)}$

其中，$g(n)={\sum }_{i=1}^n\lfloor \frac{n}{i}\rfloor$。

因为要用到多次$g(n)$，所以，可以提前把$g(1)$~$g(N)$都预处理出来。

代码：

#include <iostream>
#include <cstdio>
#include <cmath>
#include <cstring>
#include <set>
#include <map>
#include <queue>
#include <stack>
#include <vector>
#include <string>
#include <algorithm>
#define INF 0x3f3f3f3f

using namespace std;

typedef long long LL;

const int N = 5e4 + 10;

int cnt, prime[N], mu[N], g[N];
bool vis[N];

void init(int x)
{
    mu[1] = 1;
    for(int i = 2; i <= x; i ++)
    {
        if(!vis[i]) prime[cnt ++] = i, mu[i] = -1;
        
        for(int j = 0; prime[j] <= x / i; j ++)
        {
            vis[prime[j] * i] = true;
            if(i % prime[j] == 0)   break;
            mu[prime[j] * i] = -mu[i];
        }
    }

    for(int i = 1; i <= x; i ++)    mu[i] += mu[i - 1];
    
    for(int i = 1; i <= x; i ++)
    {
        for(int l = 1, r; l <= i; l = r + 1)
        {
            r = min(i, i / (i / l));
            g[i] += (r - l + 1) * ( i / l);
        }
    }
}

int main()
{
#ifdef LOCAL
    freopen("in.in", "r", stdin);
    freopen("out.out", "w", stdout);
#endif

    init(N - 1);

    int T;  scanf("%d", &T);
    while(T --)
    {
        int n, m;
        scanf("%d%d", &n, &m);
        
        int t = min(n, m);
        LL ans = 0;
        for(int l = 1, r; l <= t; l = r + 1)
        {
            r = min(t, min(n / (n / l), m / (m / l)));
            ans += (LL)(mu[r] - mu[l - 1]) * g[n / l] * g[m / l];
        }

        printf("%lld\n", ans);
    }

    return 0;
}