问题
\[\large \sum_{i=1}^n\sum_{j=1}^mlcm(i,j)
\]
结论
结论一：
\[\large \sum_{d=1}^nd\cdot \sum_{k=1}^{\lfloor\frac nd\rfloor}\mu(k)\cdot k^2\cdot\dfrac {{\lfloor\frac n{dk}\rfloor}({\lfloor\frac n{dk}\rfloor}+1)}{2}\dfrac {{\lfloor\frac m{dk}\rfloor}({\lfloor\frac m{dk}\rfloor}+1)}{2}
\]
结论二：（进一步优化）
\[\large \sum_{k=1}^n\lfloor\frac nk\rfloor^2\sum_{d\mid k}\varphi(d)\mu(\frac kd)
\]
分析
先将式子化简：
\[\large
\begin{split}
\sum_{i=1}^n\sum_{j=1}^mlcm(i,j)
&=\sum_{i=1}^n\sum_{j=1}^{m}\sum_{d\vert i,d\vert j,\gcd(\frac id,\frac jd)=1}\frac {ij}{d}\\
&=\sum_{d=1}^nd\sum_{i=1}^{\lfloor\frac nd\rfloor}\sum_{j=1}^{\lfloor\frac md\rfloor}ij[\gcd(i,j)=1]\\
\end{split}
\]
这时我们单独计算后半部分，这里没有使用莫比乌斯反演而是直接用 \(\large \mu*I=\varepsilon\) 的性质：（这部分其实就是上一个问题）
\[\large
\begin{split}
S(n,m)
&=\sum_{i=1}^{n}\sum_{j=1}^{m}ij[\gcd(i,j)=1]\\
&=\sum_{i=1}^{n}\sum_{j=1}^{m}ij\sum_{d\vert \gcd(i,j)}\mu(d)\\
&=\sum_{d=1}^{n}\mu(d)\cdot d^2\sum_{i=1}^{\lfloor\frac nd\rfloor}\sum_{j=1}^{\lfloor\frac md\rfloor}ij\\
&=\sum_{d=1}^{n}\mu(d)\cdot d^2\sum_{i=1}^{\lfloor\frac nd\rfloor}i\sum_{j=1}^{\lfloor\frac md\rfloor}j\\
&=\sum_{d=1}^{n}\mu(d)\cdot d^2\cdot G(\lfloor\frac nd\rfloor)G(\lfloor\frac md\rfloor)
\end{split}
\]
其中，\(\large G(x)=\frac {x(x+1)}{2}\) ，表示等差数列求和。
很容易发现，对于 \(\large S\) 函数，可以通过整除分块来计算，按照 \(\large\lfloor\frac nd\rfloor\) 以及 \(\large \lfloor\frac md\rfloor\) 来分块即可，前面的部分就是预处理一下前缀和即可。
现在代回原式子：
\[\large \sum_{d=1}^nd\cdot S({\lfloor\frac nd\rfloor},{\lfloor\frac md\rfloor})
\]
然后我们可以发现，原式还需要一次除法分块，而第二次除法分块将借助第一次的作为上界进行分块，这就是二次除法分块，这样做的复杂度可以证明是 \(\large O(n^{\frac 34})\) 的。
接下来考虑进一步优化：
\[\large \sum_{d=1}^nd\cdot \sum_{k=1}^{\lfloor\frac nd\rfloor}\mu(k)\cdot k^2\cdot G(\lfloor\frac n{dk}\rfloor)G(\lfloor\frac m{dk}\rfloor)
\]
考虑枚举 \(\large dk\) ：
\[\large \sum_{d=1}^nd\cdot \sum_{d\vert x}\mu(\frac xd)\cdot (\frac xd)^2\cdot G(\lfloor\frac n{x}\rfloor)G(\lfloor\frac m{x}\rfloor)
\]
交换求和符号（先枚举一个数再枚举其倍数等价于先枚举其倍数再枚举倍数的因数）：
\[\large \sum_{x=1}^n G(\lfloor\frac n{x}\rfloor)G(\lfloor\frac m{x}\rfloor)\sum_{d\vert x}\mu(\frac xd)\cdot \frac {x^2}d
\]
考虑先枚举 \(\large \frac xd\) ：
\[\large \sum_{x=1}^n G(\lfloor\frac n{x}\rfloor)G(\lfloor\frac m{x}\rfloor)\sum_{k\vert x}kx\cdot \mu(k)
\]
直接把 \(\large x\) 提出来：
\[\large \sum_{x=1}^n G(\lfloor\frac n{x}\rfloor)G(\lfloor\frac m{x}\rfloor)x\sum_{k\vert x}k\cdot \mu(k)
\]
于是后面那部分就是 \(\large F(x)=\sum\limits_{d\vert x}id(x)\mu(x)\) ，因为 \(id\) 和 \(\mu\) 都是积性函数，所以 \(F\) 也是积性函数。
所以可以考虑线性筛出来，前面部分就是一个整除分块，直接做即可。
时间复杂度是 \(\large O(V+m\sqrt{V})\) ，其中 \(V\) 是值域，\(m\) 是询问次数，如果使用杜教筛可以优化到 \(\large O(V^{\frac 23}+m\sqrt{V})\) 。
代码
#include<bits/stdc++.h>
using namespace std;
template <typename T>
inline void read(T &x){
x=0;bool f=false;char ch=getchar();
while(!isdigit(ch)){f|=ch=='-';ch=getchar();}
while(isdigit(ch)){x=x*10+(ch^48);ch=getchar();}
x=f?-x:x;
return ;
}
template <typename T>
inline void write(T x){
if(x<0) x=-x,putchar('-');
if(x>9) write(x/10);
putchar(x%10^48);
return ;
}
#define int long long
#define ll long long
#define ull unsigned long long
#define rep(i,x,y) for(int i=(x);i<=(y);i++)
#define dep(i,y,x) for(int i=(y);i>=(x);i--)
const int N=1e7+5,M=2e5+5,MOD=1e8+9;
int n,m,k;
int prime[N],f[N],cnt,pre[N],low[N];
ll val[N],Ans;
bool isprime[N];
inline int inc(int x,int y){x+=y;return x>=MOD?x-MOD:x;}
inline void incc(int &x,int y){x+=y;if(x>=MOD) x-=MOD;}
bool vis[N];
inline void GetPrimes(int d){
vis[1]=true;f[1]=low[1]=pre[1]=1;//对1进行定义
for(ll i=2;i<=d;i++){
if(!vis[i]) prime[++cnt]=low[i]=i,f[i]=1-i;//对质数进行定义
pre[i]=inc(pre[i-1],i*f[i]%MOD);
for(ll j=1;j<=cnt&&i*prime[j]<=d;j++){
vis[i*prime[j]]=true;
if(i%prime[j]==0){
low[i*prime[j]]=low[i]*prime[j];
if (low[i]==i) f[i*prime[j]]=1-prime[j];//对质数的若干次幂进行定义(一般由f[i]递推);
else f[i*prime[j]]=f[i/low[i]]*f[low[i]*prime[j]];//把当前质因子的所有项剃掉，然后把新的加回来
break;
}
low[i*prime[j]]=prime[j];
f[i*prime[j]]=f[i]*f[prime[j]];//通过 f(n),f(m) 计算 f(nm)
}
}
return ;
}
inline int GetSum(int x){return (1ll*x*(x+1)/2)%MOD;}
inline int Sum(int x,int y){
int res=0;
if(x>y) swap(x,y);
for(int l=1,r;l<=x;l=r+1){
r=min(x/(x/l),y/(y/l));
incc(res,1ll*(inc((pre[r]-pre[l-1])%MOD,MOD))*GetSum(x/l)%MOD*GetSum(y/l)%MOD);
}
return res;
}
inline int Calc(int x,int y){
int res=0;
for(int l=1,r;l<=x;l=r+1){
r=min(x/(x/l),y/(y/l));
incc(res,1ll*inc((pre[r]-pre[l-1])%MOD,MOD)%MOD*GetSum(x/l)%MOD*GetSum(y/l)%MOD);
}
return res;
}
signed main(){
GetPrimes(1e7);
int t;read(t);
while(t--){
read(n),read(m);
if(n>m) swap(n,m);
write(Calc(n,m)),putchar('\n');
}
return 0;
}