更好的阅读体验(LaTeX较多，可能加载比较慢)

题面

设 $d(x)$ 为 $x$ 的约数个数，求
$$ \sum_{i=1}^n\sum_{j=1}^md(ij) $$

Solution

引理

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

Proof

一个很显然的想法是，枚举 $i,j$ 中各自的因子，然后相乘得到 $ij$ 的一个因子。但是这样会算重。

考虑如何计数。首先，用算数基本定理表示 $i,j$ ：
$$ i=\prod p_k^{a_k}\\ j=\prod p_k^{b_k}\\ $$
那么 $i\times j$ 就是：
$$ i\times j=\prod p_k^{a_k+b_k} $$
考虑如何枚举 $p_k$ 的所有幂次，设当前枚举的幂次为 $c_k$ .

显然，我们可以 对无序的枚举钦定一个顺序 ，那么有如下情况：

• $c_k\leq a_k$ 我们假定对于每个约数，都优先选择 $a$ 中的部分。那么这里即是枚举 $a$ 中所选的幂次。
• $c_k>a_k$ 这时候必须用到 $b$ ，说明 $a$ 中默认已经选满了，直接枚举的是 $b$ 中的幂次。

显然，在这样的讨论下，无论何时 $a,b$ 都不可能同时 被枚举 。

设当前 $i,j$ 中枚举的因子分别为 $x|i,y|j$ ，那么就相当于满足 $\gcd(x,y)=1$ .

因此，有等式：
$$ d(ij)=\sum_{x|i}\sum_{y|j}[\gcd(x,y)=1] $$
Q.E.D.

由引理，有：
$$ \begin{aligned} \sum_{i=1}^n\sum_{j=1}^md(ij) & =\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[\gcd(x,y)=1]\sum_{i=1}^n\sum_{j=1}^m[x|i,y|j]\\\\ &=\sum_{x=1}^n\sum_{y=1}^m\sum_{d|\gcd(x,y)}\mu(d)\Big\lfloor \frac{n}{x}\Big\rfloor \Big\lfloor \frac{m}{y}\Big\rfloor\\\\ &=\sum_{d=1}^n\mu(d)\sum_{x=1}^n [d|x]\sum_{y=1}^m[d|y]\Big\lfloor \frac{n}{x}\Big\rfloor \Big\lfloor \frac{m}{y}\Big\rfloor\\\\ &=\sum_{d=1}^n\mu(d)\sum_{x=1}^{\lfloor n/d\rfloor}\sum_{y=1}^{\lfloor m/d\rfloor}\Big\lfloor \frac{n}{xd}\Big\rfloor \Big\lfloor \frac{m}{yd}\Big\rfloor\\\\ &=\sum_{d=1}^n\mu(d)\sum_{x=1}^{\lfloor n/d\rfloor}\Big\lfloor \frac{n}{xd}\Big\rfloor\sum_{y=1}^{\lfloor m/d\rfloor} \Big\lfloor \frac{m}{yd}\Big\rfloor\\\\ \end{aligned} $$
令 $T_1=xd,T_2=yd$ ，有
$$ \sum_{d=1}^n\mu(d)\sum_{x=1}^{\lfloor n/d\rfloor}\Big\lfloor \frac{n}{xd}\Big\rfloor\sum_{y=1}^{\lfloor m/d\rfloor} \Big\lfloor \frac{m}{yd}\Big\rfloor=\sum_{d=1}^n\mu(d)\sum_{T_1=1}^n\Big\lfloor \frac{n}{T_1}\Big\rfloor\sum_{T_2=1}^m \Big\lfloor \frac{m}{T_2}\Big\rfloor $$
后面两项都是标准的整除分块式子，第一项就是 $\mu$ 函数的前缀和，直接求就好了。

//Author: RingweEH
const int N=5e4+10;
int n,m,mu[N],pri[N],tot=0,sum[N];
ll f[N];
bool is[N];

void init()
{
    is[1]=mu[1]=1;
    for ( int i=2; i<=(N-10); i++ )
    {
        if ( !is[i] ) { pri[++tot]=i; mu[i]=-1; }
        for ( int j=1; j<=tot && (i*pri[j])<=(N-10); j++ )
        {
            is[i*pri[j]]=1;
            if ( i%pri[j] ) mu[i*pri[j]]=mu[i]*mu[pri[j]];
            else { mu[i*pri[j]]=0; break;}
        }
    }
    sum[0]=0;
    for ( int i=1; i<=(N-10); i++ )
        sum[i]=sum[i-1]+mu[i];
    memset( f,0,sizeof(f) );
    for ( int i=1; i<=(N-10); i++ )
        for ( int l=1,r=0; l<=i; l=r+1 )
        {
            r=i/(i/l);
            f[i]=f[i]+1ll*(r-l+1)*(i/l);
        }
}

int main()
{
    int T=read(); init();
    while ( T-- )
    {
        n=read(); m=read();
        if ( n>m ) swap( n,m );
        ll ans=0;
        for ( int l=1,r=0; l<=n; l=r+1 )
        {
            r=min( n/(n/l),m/(m/l) );
            ans+=(sum[r]-sum[l-1])*f[n/l]*f[m/l];
        }
        printf( "%lld\n",ans );
    }
}