首先有一个重要的式子
ϕ
(
i
j
)
=
ϕ
(
i
)
ϕ
(
j
)
gcd
(
i
,
j
)
ϕ
(
gcd
(
i
,
j
)
\phi(ij)=\frac{\phi(i)\phi(j)\gcd(i,j)}{\phi(\gcd(i,j)}
ϕ(ij)=ϕ(gcd(i,j)ϕ(i)ϕ(j)gcd(i,j)
这题就直接这么推
=
∑
d
=
1
d
ϕ
(
d
)
∑
i
=
1
n
/
d
∑
j
=
1
m
/
d
ϕ
(
i
d
)
ϕ
(
j
d
)
[
gcd
(
i
,
j
)
=
1
]
=\sum\limits_{d=1}\frac{d}{\phi(d)}\sum\limits_{i=1}^{n/d}\sum\limits_{j=1}^{m/d}\phi(id)\phi(jd)[\gcd(i,j)=1]
=d=1∑ϕ(d)di=1∑n/dj=1∑m/dϕ(id)ϕ(jd)[gcd(i,j)=1]
=
∑
d
=
1
d
ϕ
(
d
)
∑
t
=
1
n
/
d
μ
(
t
)
∑
i
=
1
n
/
d
t
∑
j
=
1
m
/
d
t
ϕ
(
i
d
t
)
ϕ
(
j
d
t
)
=\sum\limits_{d=1}\frac{d}{\phi(d)}\sum\limits_{t=1}^{n/d}\mu(t)\sum\limits_{i=1}^{n/dt}\sum\limits_{j=1}^{m/dt}\phi(idt)\phi(jdt)
=d=1∑ϕ(d)dt=1∑n/dμ(t)i=1∑n/dtj=1∑m/dtϕ(idt)ϕ(jdt)
设
T
=
d
t
T=dt
T=dt
=
∑
T
=
1
n
∑
d
∣
T
μ
(
d
/
t
)
d
ϕ
(
d
)
∑
i
=
1
n
/
T
ϕ
(
i
T
)
∑
j
=
1
m
/
T
ϕ
(
j
T
)
=\sum\limits_{T=1}^n\sum\limits_{d|T}\frac{\mu(d/t)d}{\phi(d)}\sum\limits_{i=1}^{n/T}\phi(iT)\sum\limits_{j=1}^{m/T}\phi(jT)
=T=1∑nd∣T∑ϕ(d)μ(d/t)di=1∑n/Tϕ(iT)j=1∑m/Tϕ(jT)
把三部分分别设一下,可以得到
=
∑
T
=
1
n
f
(
T
)
g
(
T
,
n
/
T
)
g
(
T
,
m
/
T
)
=\sum\limits_{T=1}^n f(T)g(T,n/T)g(T,m/T)
=T=1∑nf(T)g(T,n/T)g(T,m/T)
f
f
f很容易预处理
考虑
g
g
g,容易发现
g
(
T
,
n
)
=
g
(
T
,
n
−
1
)
+
ϕ
(
T
n
)
g(T,n)=g(T,n-1)+\phi(Tn)
g(T,n)=g(T,n−1)+ϕ(Tn)
然后再把整条式子设出来,
h ( a , b , n ) = ∑ T = 1 n f ( T ) g ( T , a ) g ( T , b ) h(a,b,n)=\sum\limits_{T=1}^nf(T)g(T,a)g(T,b) h(a,b,n)=T=1∑nf(T)g(T,a)g(T,b)
考虑阈值分治,设一个阈值
B
B
B
预处理
h
(
a
,
b
,
n
)
,
a
≤
B
,
b
≤
B
h(a,b,n),a\le B,b\le B
h(a,b,n),a≤B,b≤B的答案
对于
n
/
T
,
m
/
T
≥
B
n/T,m/T\ge B
n/T,m/T≥B的直接暴力跑
剩下的整除分块
时间复杂度
O
(
n
ln
n
+
n
B
2
+
T
(
n
+
n
B
)
O(n\ln n+nB^2+T(\sqrt{n}+\frac{n}{B})
O(nlnn+nB2+T(n+Bn)
块大小我取的是
B
=
50
B=50
B=50
code:
#include<bits/stdc++.h>
#define N 100005
#define B 50
#define mod 998244353
#define ll long long
using namespace std;
ll qpow(ll x, ll y) {
ll ret = 1;
for(; y; y >>= 1, x = x * x % mod) if(y & 1) ret = ret * x % mod;
return ret;
}
int mu[N], phi[N], f[N], iphi[N];
vector<int> g[N], S[B + 5][B + 5];
int p[N], sz, vis[N];
void init(int n) {
mu[1] = phi[1] = 1;
for(int i = 2; i <= n; i ++) {
if(!vis[i]) mu[i] = (mod - 1), p[++ sz] = i, phi[i] = i - 1;
for(int j = 1; j <= sz && i * p[j] <= n; j ++) {
vis[p[j] * i] = 1;
if(i % p[j] == 0) {
phi[i * p[j]] = 1ll * phi[i] * p[j] % mod;
break;
}
phi[i * p[j]] = 1ll * phi[i] * phi[p[j]] % mod, mu[i * p[j]] = mod - mu[i];
}
}
for(int i = 1; i <= n; i ++) iphi[i] = qpow(phi[i], mod - 2);
for(int i = 1; i <= n; i ++)
for(int j = 1; j <= n / i; j ++)
f[i * j] = (f[i * j] + 1ll * i * mu[j] % mod * iphi[i] % mod) % mod;
for(int T = 1; T <= n; T ++) {
g[T].resize(n / T + 2); g[T][0] = 0;
for(int i = 1; i <= n / T; i ++) {
g[T][i] = (g[T][i - 1] + phi[i * T]) % mod;
}
}
for(int a = 1; a <= B; a ++)
for(int b = 1; b <= B; b ++) { int lim = n / min(a, b);
S[a][b].resize(lim + 2), S[a][b][0] = 0;
for(int T = 1; T <= lim; T ++) {
S[a][b][T] = (S[a][b][T - 1] + 1ll * f[T] * g[T][a] % mod * g[T][b] % mod) % mod;
}
}
}
ll calc(int n, int m) {
ll ans = 0;
for(int i = 1; i <= m / B; i ++) ans = (ans + 1ll * f[i] * g[i][n / i] % mod * g[i][m / i] % mod) % mod;
for(int l = m / B + 1, r = 0; l <= n; l = r + 1) {
r = min(n / (n / l), m / (m / l));
ans = (ans + S[n / l][m / l][r] - S[n / l][m / l][l - 1] + mod) % mod;
}
return ans;
}
int n, m, t;
int main() {
init(N - 5);
scanf("%d", &t);
while(t --) {
scanf("%d%d", &n, &m); if(n > m) swap(n, m);
printf("%lld\n", calc(n, m));
}
return 0;
}