这正是 $p$ 在 $\sigma(ij)$ 中的局部贡献。
所有质数的贡献相乘后得到 $\prod_p\sigma\left(p^{\alpha_p+\beta_p}\right)
=\sigma(ij)$。
因此恒等式成立。
:::
代入原式并交换求和：$\operatorname{Ans}(n,m)=\sum_{d\le\min(n,m)}\mu(d)dF_d(n)F_d(m)$，
其中：$F_d(x)=\sum_{u=1}^{\left\lfloor x/d\right\rfloor}\varphi(du)\sigma(u)$。
设所有询问均满足 $n\le m$，最大值分别为 $X,Y$。取阈值 $B\approx\sqrt{\frac{XY}{T}}$。
对于 $d\le B$，预处理 $F_d(x)$ 的前缀和，每次询问直接累加。
对于 $d>B$，按较小端点 $n$ 离线扫描。扫描到 $x$ 时，枚举满足 $d\mid x$ 的大因子，并把 $\mu(d)d\varphi(x)\sigma\left(\frac{x}{d}\right)\varphi(dv)\sigma(v)$ 加入位置 $dv$。这样查询纵坐标前缀 $m$，即可得到所有 $i\le n,j\le m$ 的贡献。单点增加与前缀查询使用根号分块维护。
时间复杂度为 $\mathcal O\left(Y\log X+TB+\frac{XY}{B}+T\sqrt Y\right)$。
```cpp
#include <bits/stdc++.h>
using namespace std;
#define ll long long
const int N = 1e5 + 5, Q = 1e5 + 5, P = 998244353, H = 9e5 + 5;
int T, qn[Q], qm[Q], ans[Q], nx[Q], head[N];
int mu[N], phi[N], sg[N], pr[N], pc, vis[N], cf[N];
int sd[N], sc, off[N], pre[H], pt;
int cnt[N], bg[N], cur[N], dv[H];
int wo[N], ww[H], wt;
int val[N], bs[N], S;
inline void add(int x, int v) {
int y = val[x] + v;
if (y >= P) y -= P;
val[x] = y;
int b = (x - 1) / S;
y = bs[b] + v;
if (y >= P) y -= P;
bs[b] = y;
}
inline int ask(int x) {
int b = (x - 1) / S;
ll s = 0;
for (int i = 0; i < b; i++) {
s += bs[i];
}
for (int i = b * S + 1; i <= x; i++) {
s += val[i];
}
return s % P;
}
int main() {
cin >> T;
int X = 0, Y = 0;
for (int i = 0; i < T; i++) {
cin >> qn[i] >> qm[i];
if (qn[i] > qm[i]) swap(qn[i], qm[i]);
X = max(X, qn[i]);
Y = max(Y, qm[i]);
}
mu[1] = phi[1] = 1;
for (int i = 2; i <= Y; i++) {
if (!vis[i]) {
pr[pc++] = i;
mu[i] = -1;
phi[i] = i - 1;
}
for (int j = 0; j < pc && i * pr[j] <= Y; j++) {
int p = pr[j], x = i * p;
vis[x] = 1;
if (i % p == 0) {
mu[x] = 0;
phi[x] = phi[i] * p;
break;
}
mu[x] = -mu[i];
phi[x] = phi[i] * (p - 1);
}
}
for (int d = 1; d <= Y; d++) {
for (int j = d; j <= Y; j += d) {
int x = sg[j] + d;
if (x >= P) x -= P;
sg[j] = x;
}
}
for (int i = 1; i <= Y; i++) {
if (mu[i]) {
cf[i] = mu[i] > 0 ? i : P - i;
}
}
int B = sqrt((long double)X * Y / T);
if (B < 1) B = 1;
if (B > X) B = X;
for (int d = 1; d <= B; d++) {
if (mu[d]) {
sd[sc++] = d;
off[d] = pt;
pre[pt++] = 0;
for (int u = 1; u <= Y / d; u++) {
pre[pt++] = (pre[pt - 1] + (ll)phi[d * u] * sg[u]) % P;
}
}
}
for (int i = 0; i < T; i++) {
ll s = 0;
for (int j = 0; j < sc && sd[j] <= qn[i]; j++) {
int d = sd[j];
int x = pre[off[d] + qn[i] / d];
int y = pre[off[d] + qm[i] / d];
s += (ll)cf[d] * x % P * y % P;
if (s >= (1LL << 62)) s %= P;
}
ans[i] = s % P;
}
for (int d = B + 1; d <= X; d++) {
if (mu[d]) {
wo[d] = wt;
ww[wt++] = 0;
for (int u = 1; u <= Y / d; u++) {
ww[wt++] = (ll)phi[d * u] * sg[u] % P;
}
for (int x = d; x <= X; x += d) {
cnt[x]++;
}
}
}
bg[1] = 0;
for (int x = 1; x <= X; x++) {
bg[x + 1] = bg[x] + cnt[x];
cur[x] = bg[x];
}
for (int d = B + 1; d <= X; d++) {
if (mu[d]) {
for (int x = d; x <= X; x += d) {
dv[cur[x]++] = d;
}
}
}
memset(head, -1, sizeof(int) * (X + 1));
for (int i = 0; i < T; i++) {
nx[i] = head[qn[i]];
head[qn[i]] = i;
}
S = sqrt(2.0 * Y) + 1;
for (int x = 1; x <= X; x++) {
for (int j = bg[x]; j < bg[x + 1]; j++) {
int d = dv[j], u = x / d, o = wo[d];
int base = (ll)cf[d] * ww[o + u] % P;
for (int v = 1, y = d; v <= Y / d; v++, y += d) {
add(y, (ll)base * ww[o + v] % P);
}
}
for (int i = head[x]; i != -1; i = nx[i]) {
ans[i] += ask(qm[i]);
if (ans[i] >= P) ans[i] -= P;
}
}
for (int i = 0; i < T; i++) {
cout << ans[i] << endl;
}
return 0;
}
```