&\sum\limits_{i=1}^{n}\sum\limits_{j=1}^{m}\varphi(ij)\\
=&\sum\limits_{d=1}^{n}\sum\limits_{i=1}^{\lfloor \frac n d\rfloor}\sum\limits_{j=1}^{\lfloor \frac m d\rfloor}\frac{\varphi(i)\varphi(j)d}{\varphi(d)}\cdot[i \bot j]\\
=&\sum\limits_{d=1}^{n}\frac{d}{\varphi(d)}\sum\limits_{k=1}^{\lfloor \frac n d\rfloor }\mu(k)\sum\limits_{i=1}^{\lfloor \frac{n}{dk}\rfloor}\sum\limits_{j=1}^{\lfloor \frac{m}{dk}\rfloor}\varphi(idk)\varphi(jdk)\\
=&\sum\limits_{T=1}^{n}(\sum\limits_{d \mid T}\frac{d}{\varphi(d)}\mu(\frac T d))(\sum\limits_{i=1}^{\lfloor \frac n T \rfloor}\varphi(iT))(\sum\limits_{i=1}^{\lfloor \frac m T \rfloor}\varphi(iT))\\
\end{aligned}
令 f(T)=\sum\limits_{d \mid T}\frac{d}{\varphi(d)}\mu(\frac T d),g(x,y)=\sum\limits_{i=1}^{y}\varphi(ix)，则答案为 \sum\limits_{T=1}^{n}f(T)g(T,\lfloor\frac n T\rfloor)g(T,\lfloor\frac m T\rfloor)。
考虑 f，我们可以筛出 \varphi,\mu，然后枚举倍数 O(n\ln n) 做。对于 g，有性质 xy \le n，所以可以暴力枚举合法的 x,y，复杂度 O(\log n\cdot\sum\limits_{i=1}^{n}\lfloor\frac n i\rfloor)=O(n \log^2 n)。因为利用最开始的性质，对于 1 \le x,y \le n 可以 O(\log n) 复杂度求出 \gcd，并求出 \varphi(xy)。
发现卷积很难做。具体的，不可能预处理出所有 T,n,m。所以只考虑预处理一部分。
选取一个 B，对于 1 \le t \le n,1 \le x,y \le B，求出 s(t,x,y)=\sum\limits_{i=1}^{t}f(i)g(i,x)g(i,y)。这样的话可以代入 T,n,m 并解决 \lfloor\frac{n}{T}\rfloor \le B 的求值。这部分要求 T \ge \lfloor \frac n B \rfloor，可以直接从此处开始进行整除分块，每次差分算贡献 s(r,\lfloor \frac n l \rfloor,\lfloor \frac m l \rfloor)-s(l-1,\lfloor \frac n l \rfloor,\lfloor \frac m l \rfloor)。对于 T \le \lfloor \frac n B \rfloor，可以直接暴力求。
概括的说，原本可以都整除分块，但是为了平衡复杂度，只预处理出固定 n,m 的上界和 T 的下界的一部分，剩下的枚举 T 暴力计算。
时间复杂度 O(n\log n+n\log^2 n+nB^2+q(\frac n B +\sqrt n))=O(n(\log^2 n+B^2+\frac n B+\sqrt n))。
由于 n,q 大致同阶，B 应当取 n^{\frac 1 3}，然而空间是 O(n+n\log n+nB^2)=O(n\log n+nB^2) 的，承受不住，所以代码中 B 取的是 30。
理论上能过，实际也能，但是觉得常数很大所以在 nB^2 的地方循环展开了一下。
```cpp
//g++ c.cpp -o c -g -std=c++14 -O0 -Wall -fsanitize=undefined
#include<iostream>
#include<cstdio>
#include<vector>
#include<algorithm>
#include<chrono>
#define LL long long
using namespace std;
const int maxr=1e6+10,maxn=1e5+10,_maxn=1e5,mod=998244353,maxsqn=33,_maxsqn=32;
int N,sqn,prime[maxn],ptop=0,vis[maxn],phi[maxn],mu[maxn],s[maxn][maxsqn][maxsqn];LL ni[maxn],f[maxn];
vector<int>g[maxn];
int stime(){return chrono::steady_clock::now().time_since_epoch().count()/1000;}
namespace rdin{
char buf[maxr],*l,*r;
inline char gtc(){if(l==r){r=(l=buf)+fread(buf,1,maxr,stdin);}return l==r?EOF:*l++;}
inline int qd(){
int rt=0;char c=gtc();
while(c<'0'||c>'9') c=gtc();
while('0'<=c&&c<='9') rt=(rt<<3)+(rt<<1)+(c^48),c=gtc();
return rt;
}
}
using rdin::qd;
LL piphi(int x,int y){int t=__gcd(x,y);return (LL)phi[x]*phi[y]%mod*t%mod*ni[phi[t]]%mod;}
int main(){
freopen("in.txt","r",stdin);int cntt=stime();
phi[1]=mu[1]=ni[0]=ni[1]=1;for(int i=2;i<=_maxn;i++) ni[i]=mod-mod/i*ni[mod%i]%mod;
for(int i=2;i<=_maxn;i++){
if(!vis[i]) prime[++ptop]=i,mu[i]=-1,phi[i]=i-1;
for(int j=1;j<=ptop&&i*prime[j]<=_maxn;j++){
vis[i*prime[j]]=1;
if(i%prime[j]) phi[i*prime[j]]=phi[i]*phi[prime[j]],mu[i*prime[j]]=-mu[i];
else{phi[i*prime[j]]=phi[i]*prime[j];break;}
}
}
// printf("%d us\n",stime()-cntt);cntt=stime();
for(int i=1;i<=_maxn;i++) for(int j=i;j<=_maxn;j+=i) (f[j]+=i*ni[phi[i]]*mu[j/i])%=mod;
// printf("%d us\n",stime()-cntt);cntt=stime();
for(int i=1;i<=_maxn;i++){//i:fac,j:times;
LL sum=0;g[i].push_back(0);
for(int j=1;i*j<=_maxn;j++) g[i].push_back((sum+=piphi(i,j))%=mod);
}
// printf("%d us\n",stime()-cntt);cntt=stime();
/*
for(int i=1;i<=_maxn;i++){
for(int j=1;j<=_maxsqn&&i*j<=_maxn;j++){
LL pi=f[i]*g[i][j]%mod;
for(int k=1;k<=_maxsqn&&i*k<=_maxn;k++) s[i][j][k]=(s[i-1][j][k]+pi*g[i][k])%mod;
}
}*/
for(int i=1;i<=_maxn;i++){
for(int j=min(_maxsqn,_maxn/i);j;j--){
int k=min(_maxsqn,_maxn/i);LL pi=f[i]*g[i][j]%mod;
for(;k&7;k--) s[i][j][k]=(s[i-1][j][k]+pi*g[i][k])%mod;
// for(;k;k-=8){
// s[i][j][k]=(s[i-1][j][k]+pi*g[i][k])%mod;s[i][j][k-1]=(s[i-1][j][k-1]+pi*g[i][k-1])%mod;
// s[i][j][k-2]=(s[i-1][j][k-2]+pi*g[i][k-2])%mod;s[i][j][k-3]=(s[i-1][j][k-3]+pi*g[i][k-3])%mod;
// s[i][j][k-4]=(s[i-1][j][k-4]+pi*g[i][k-4])%mod;s[i][j][k-5]=(s[i-1][j][k-5]+pi*g[i][k-5])%mod;
// s[i][j][k-6]=(s[i-1][j][k-6]+pi*g[i][k-6])%mod;s[i][j][k-7]=(s[i-1][j][k-7]+pi*g[i][k-7])%mod;
// }
int *p1=&s[i][j][k],*p2=&s[i-1][j][k],*p3=&g[i][k];//we can find that when only use both will be fast
// for(;k;k--,p1--,p2--,p3--) p1[0]=(p2[0]+pi*p3[0])%mod;
for(;k;k-=8,p1-=8,p2-=8,p3-=8){
p1[0]=(p2[0]+pi*p3[0])%mod;p1[-1]=(p2[-1]+pi*p3[-1])%mod;p1[-2]=(p2[-2]+pi*p3[-2])%mod;p1[-3]=(p2[-3]+pi*p3[-3])%mod;
p1[-4]=(p2[-4]+pi*p3[-4])%mod;p1[-5]=(p2[-5]+pi*p3[-5])%mod;p1[-6]=(p2[-6]+pi*p3[-6])%mod;p1[-7]=(p2[-7]+pi*p3[-7])%mod;
}
}
}
// printf("%d us\n",stime()-cntt);cntt=stime();
int T=qd();
while(T--){
int N=qd(),M=qd(),i=1;LL ans=0;if(N<M) swap(N,M);
for(;i<=N/_maxsqn;i++) (ans+=f[i]*g[i][N/i]%mod*g[i][M/i]%mod)%=mod;
for(int l=i,r;l<=M;l=r+1){r=min(N/(N/l),M/(M/l));(ans+=s[r][N/l][M/l]-s[l-1][N/l][M/l])%=mod;}
printf("%lld\n",(ans+mod)%mod);
}
return 0;
}
```