约数个数和问题

解题思路：

本题要求对于 $i$ 从 $1$ 到 $N$，$j$ 从 $1$ 到 $M$ 中的所有 $ij$ 求一下他们的约数个数的和。

首先先做一步转化，$d(ij)=\sum_{x\vert i}\sum_{y \vert j} [gcd(x,y)=1]$，即 $i$ 的所有约数和 $j$ 的所有约数搭配组成数对，所有最大公因数是 $1$ 的数对个数就是 $ij$ 的约数个数。

这里可以简单证明一下为什么是这样的，我们先将 $i,j$ 分解质因数，$i=p_1^{\alpha_1}p_2^{\alpha_2} … p_k^{\alpha_k}$，$j=p_1^{\beta_1}p_2^{\beta_2}…p_k^{\beta_k}$，则 $ij=p_1^{\alpha_1+\beta_1}p_2^{\alpha_2+\beta_2}…p_k^{\alpha_k+\beta_k}$。因此 $ij$ 的约数个数就是 $(\alpha_1+\beta_1+1)(\alpha_2+\beta_2+1)…(\alpha_k+\beta_k+1)$。然后我们再看一下 $\sum_{x\vert i}\sum_{y \vert j} [gcd(x,y)=1]$ 是多少，我们先看一下只选 $p_1$ 有多少种取值，由于要让 $x,y$ 互质，一种取法就是 $x,y$ 都不取，都是 $1$，然后要么是 $x$ 取要么是 $y$ 取，$x$ 可以取 $1$ 个 $p_1$，也可以取 $2$ 个 $p_1$，以此类推可以取 $\alpha_1$ 个 $p_1$，同理 $y$ 最多可以取到 $\beta_1$ 个 $p_1$，因此对于 $p_1$ 一共有 $\alpha_1+\beta_1+1$ 种取法，以此类推发现这个公式的方案数就是约数个数。

我们原本要求 $\sum_{i=1}^N\sum_{j=1}^M d(i,j)$，现在通过我们的转化，变成了 $\sum_{i=1}^N\sum_{j=1}^M \sum_{x\vert i}\sum_{y \vert j} [gcd(x,y)=1]$，有了这个形式的式子，我们就可以用莫比乌斯反演去求了。

设 $F(n)=\sum_{i=1}^N\sum_{j=1}^M \sum_{x\vert i}\sum_{y \vert j} [n \vert gcd(x,y)]$，也就是最大公约数能整除 $n$ 的所有数对个数。设 $f(n)=\sum_{i=1}^N\sum_{j=1}^M \sum_{x\vert i}\sum_{y \vert j} [gcd(x,y)=n]$，也就是最大公约数恰好等于 $n$ 的所有数对个数。可以发现此时$F(n)$ 和 $f(n)$ 的关系就是 $F(n)=\sum_{n \vert d} f(d)$，这就是莫比乌斯反演的条件，我们就可以得出 $f(n)=\sum_{n \vert d} \mu{(\frac{d}{n})} F(d)$，而我们要求的其实就是 $f(1)$，则 $f(1)=\sum_{d=1}^{min(N,M)} \mu{(d)} F(d)$

接下来我们就需要考虑 $F(d)$ 怎么求，由于 $F(d)$ 中有四个求和，因此我们可以先考虑交换求和次序，我们将 $x,y$ 的求和移到前面来，$x$ 是 $i$ 的约数，$i$ 是 $1 \sim N$，因此 $x$ 也就是 $1 \sim N$，$y$ 同理就是 $1 \sim M$，然后对应的分析一下 $i,j$ 有哪些取值，可以发现 $i,j$ 的取值对条件是没有任何影响的，因此只要 $x,y$ 确定了，那么 $d \vert gcd(x,y)$ 就确定了，而 $i$ 需要是 $x$ 的倍数，$j$ 需要是 $y$ 的倍数，因此只有在 $i$ 是 $x$ 的倍数且 $j$ 是 $y$ 的倍数时，这个情况我们才会枚举到，那么对于此时的 $x,y$，对应枚举到的次数就应该看 $i$ 有多少种 $x$ 的倍数，$j$ 有多少种 $y$ 的倍数，乘在一起就是这对 $x,y$ 总共被枚举到的次数，也就是 $\lfloor \frac{N}{x} \rfloor \lfloor \frac{M}{y} \rfloor$，这样我们要求的就变成了 $\sum_{x=1}^N \sum_{y=1}^M \lfloor \frac{N}{x} \rfloor \lfloor \frac{M}{y} \rfloor [d \vert gcd(x,y)]$。

可以发现只有 $gcd(x,y)$ 是 $d$ 的倍数时，这对 $x,y$ 才会是非零项，其他情况都是为零的不用考虑，因此我们可以令 $x’=\frac{x}{d},y’=\frac{y}{d}$，这样我们就不用去枚举 $x,y$，我们直接取枚举倍数，那么上式变成 $\sum_{x’=1}^{\frac{N}{d}}\sum_{y’=1}^{\frac{M}{d}} \lfloor \frac{N}{dx’} \rfloor \lfloor \frac{M}{dy’} \rfloor$，由于我们直接枚举倍数，所以此时 $d\vert gcd(x,y)$ 一定成立，因此这一部分就是答案。

然后我们进一步分析，我们令 $x’=i,y’=j,N’=\frac{N}{d},M’=\frac{M}{d}$，则变成 $\sum_{i=1}^{N’}\sum_{j=1}^{M’} \lfloor \frac{N’}{i} \rfloor \lfloor \frac{M’}{j} \rfloor$，此时我们再将 $\frac{N’}{i}$ 看作 $a_i$，将 $\frac{M’}{j}$ 看作 $b_j$，此时整个式子拆开其实就是 $(a_1b_1+a_1b_2+…+a_1b_{M’})+(a_2b_1+a_2b_2+…a_2b_{M’})+…+(a_{N’}b_1+a_{N’}b_2+…+a_{N’}b_{M’})$，我们通过乘法分配律可以变成 $(a_1+a_2+…+a_{N’})(b_1+b_2+…+b_{M’})$，通过这个思路，我们将原式继续转换得到 $(\sum_{i=1}^{N’} \lfloor \frac{N’}{i} \rfloor) (\lfloor \sum_{j=1}^{M’} \frac{M’}{j} \rfloor)$

我们令 $h(k)=\sum_{i=1}^{k} \lfloor \frac{k}{i} \rfloor$，最终我们要求的就是 $h(N’)h(M’)$。因此 $f(1)=\sum_{d=1}^{min(N,M)}\mu{(d)} h(\frac{N}{d}) h(\frac{M}{d})$。

而 $h(k)$ 我们可以先预处理出来，$h(k)=\sum_{i=1}^{k} \lfloor \frac{k}{i} \rfloor$，这里运用整数分块的原理可知 $\frac{k}{i}$ 最多只有 $O(\sqrt{N})$ 种取值，并且每种取值都是一段连续的区间，因此我们可以直接枚举每段区间，每一段统一的计算即可。这样每个 $h(k)$ 都能用 $O(\sqrt{N})$ 的时间复杂度求出来，我们就可以直接预处理所有的 $h(k)$。

而在求 $f(1)=\sum_{d=1}^{min(N,M)}\mu{(d)} h(\frac{N}{d}) h(\frac{M}{d})$ 时，$d$ 需要从 $1$ 枚举到 $min(N,M)$，而 $h(\frac{N}{d}) h(\frac{M}{d})$ 的取值都取决于 $d$，而这里同样可以用整数分块的思想将求和分成若干个区间来处理，在每个区间中 $h(\frac{N}{d}) h(\frac{M}{d})$ 的值是固定的，因此只需要求一下区间中 $\mu{(d)}$ 的和，这个可以预处理一下前缀和来求，这样我们就能快速的计算出 $f(1)$

C++ 代码

#include <iostream>
#include <cstring>

using namespace std;

typedef long long LL;

const int N = 50010;

int primes[N], cnt; //存储质数
bool st[N]; //记录每个数是否是合数
int mu[N]; //莫比乌斯函数
int sum[N]; //莫比乌斯函数的前缀和
int h[N];

int g(int k, int x) //返回 x 所在的区间的最后一个数
{
    return k / (k / x);
}

void init(int n) //预处理
{
    //预处理莫比乌斯函数
    mu[1] = 1;
    for(int i = 2; i <= n; i++)
    {
        if(!st[i]) primes[cnt++] = i, mu[i] = -1;
        for(int j = 0; primes[j] <= n / i; j++)
        {
            st[primes[j] * i] = true;
            if(i % primes[j] == 0) break;
            mu[primes[j] * i] = -mu[i];
        }
    }
    for(int i = 1; i <= n; i++) sum[i] = sum[i - 1] + mu[i]; //预处理前缀和

    //预处理 h(x) 函数
    for(int i = 1; i <= n; i++)
    {
        for(int l = 1, r; l <= i; l = r + 1)
        {
            r = min(i, g(i, l));
            h[i] += (r - l + 1) * (i / l);
        }
    }
}

int main()
{
    init(N - 1);

    int T;
    scanf("%d", &T);
    while(T--)
    {
        int n, m;
        scanf("%d%d", &n, &m);
        LL res = 0; //记录答案
        int k = min(n, m);
        for(int l = 1, r; l <= k; l = r + 1)
        {
            r = min(k, min(g(n, l), g(m, l)));
            res += (LL)(sum[r] - sum[l - 1]) * h[n / l] * h[m / l];
        }
        printf("%lld\n", res);
    }
    return 0;
}