题目链接
题解
式子不难推，分块打表真的没想到
首先考虑如何拆开(varphi(ij))
考虑公式
[varphi(ij) = ijprodlimits_{p | ij}frac{p - 1}{p}
]
而
[egin{aligned}
varphi(i)varphi(j) &= iprodlimits_{p | i}frac{p - 1}{p} j prodlimits_{p | j}frac{p - 1}{p} \
varphi(i)varphi(j)&= ij prodlimits_{p | ij}frac{p - 1}{p} prodlimits_{p | (i,j)}frac{p - 1}{p} \
varphi(i)varphi(j)(i,j)&= ij prodlimits_{p | ij}frac{p - 1}{p} (i,j)prodlimits_{p | (i,j)}frac{p - 1}{p} \
varphi(i)varphi(j)(i,j)&= varphi(i,j)varphi((i,j)) \
varphi(i,j) &= frac{varphi(i)varphi(j)(i,j)}{varphi((i,j))}
end{aligned}
]
所以我们有
[egin{aligned}
ans &= sumlimits_{i = 1}^{n}sumlimits_{j = 1}^{m} frac{varphi(i)varphi(j)(i,j)}{varphi((i,j))} \
&= sumlimits_{d = 1}^{n}frac{d}{varphi(d)} sumlimits_{d|i}^{n} varphi(i) sumlimits_{d|j}^{m} varphi(j) [(i,j) = d] \
&= sumlimits_{z = 1}^{n}frac{z}{varphi(z)} sumlimits_{z | d} mu(frac{d}{z}) sumlimits_{d|i}^{n} varphi(i) sumlimits_{d|j}^{m} varphi(j) \
&= sumlimits_{z = 1}^{n}frac{z}{varphi(z)} sumlimits_{z | d} mu(frac{d}{z}) sumlimits_{i = 1}^{lfloor frac{n}{d}
floor} varphi(id) sumlimits_{j = 1}^{lfloor frac{m}{d}
floor} varphi(jd) \
&= sumlimits_{d = 1}^{n}sumlimits_{z | d} frac{z}{varphi(z)}mu(frac{d}{z}) sumlimits_{i = 1}^{lfloor frac{n}{d}
floor} varphi(id) sumlimits_{j = 1}^{lfloor frac{m}{d}
floor} varphi(jd)
end{aligned}
]
隐约感觉已经差不多了
令
[F(d) = sumlimits_{z | d} frac{z}{varphi(z)}mu(frac{d}{z})
]
[G(n,d) = sumlimits_{k = 1}^{n}varphi(kd)
]
那么
[ans = sumlimits_{d = 1}^{n}F(d)G(lfloor frac{n}{d}
floor,d)G(lfloor frac{m}{d}
floor,d)
]
注意到(lfloor frac{n}{d}
floor d le n)，(G(n,d))的状态数是(O(n))的，可以预处理
(F(d))可以(O(nlogn))预处理
那么现在我们就可以(O(n))计算辣~
但这还不够，却已经无法从式子上进行优化了
分块！
我们设
[T(n,x,y) = sumlimits_{d = 1}^{n}F(d)G(x,d)G(y,d)
]
我们就可以利用整除分块(T(nxt,lfloor frac{n}{i}
floor,lfloor frac{m}{i}
floor) - T(i - 1,lfloor frac{n}{i}
floor,lfloor frac{m}{i}
floor))进行(O(sqrt{n}))的计算
所以我们只需预处理出(T(n,x,y))
但状态数很多，考虑只打出(B)项(x,y)
考虑到(nx,ny le 10^5)，时间空间复杂度(O(nB))
(y ge B)的(n)只有(lfloor frac{n}{B}
floor)项，暴力计算
所以询问时复杂度(O(sqrt{n} + lfloor frac{n}{B}
floor))
令(Tfrac{n}{B} = nB)，则(B = sqrt{T})
(B)取(100)左右就差不多了
复杂度(O(nlogn + T(sqrt{n} + lfloor frac{n}{B} floor)))
#include<algorithm>
#include<iostream>
#include<cstdlib>
#include<cstring>
#include<cstdio>
#include<vector>
#include<queue>
#include<cmath>
#include<map>
#define LL long long int
#define REP(i,n) for (int i = 1; i <= (n); i++)
#define Redge(u) for (int k = h[u],to; k; k = ed[k].nxt)
#define cls(s,v) memset(s,v,sizeof(s))
#define mp(a,b) make_pair<int,int>(a,b)
#define cp pair<int,int>
using namespace std;
const int maxn = 100005,maxm = 100005,INF = 0x3f3f3f3f;
inline int read(){
int out = 0,flag = 1; char c = getchar();
while (c < 48 || c > 57){if (c == '-') flag = 0; c = getchar();}
while (c >= 48 && c <= 57){out = (out << 1) + (out << 3) + c - 48; c = getchar();}
return flag ? out : -out;
}
const int B = 100,N = 100000,P = 998244353;
int p[maxn],pi,isn[maxn],mu[maxn],phi[maxn],inv[maxn];
int F[maxn],*G[maxn],*T[B + 1][B + 1];
void init(){
inv[0] = inv[1] = 1;
mu[1] = phi[1] = 1;
for (int i = 2; i <= N; i++) inv[i] = 1ll * (P - P / i) * inv[P % i] % P;
for (int i = 2; i <= N; i++){
if (!isn[i]) p[++pi] = i,mu[i] = P - 1,phi[i] = i - 1;
for (int j = 1; j <= pi && i * p[j] <= N; j++){
isn[i * p[j]] = true;
if (i % p[j] == 0){
phi[i * p[j]] = p[j] * phi[i];
break;
}
phi[i * p[j]] = (p[j] - 1) * phi[i];
mu[i * p[j]] = (P - mu[i]) % P;
}
}
for (int i = 1; i <= N; i++)
for (int j = i; j <= N; j += i)
F[j] = (F[j] + 1ll * i * inv[phi[i]] % P * mu[j / i] % P) % P;
for (int d = 1; d <= N; d++){
int len = N / d + 1;
G[d] = new int[len]; G[d][0] = 0;
for (int i = 1; i < len; i++)
G[d][i] = (G[d][i - 1] + phi[i * d]) % P;
}
for (int x = 1; x < B; x++)
for (int y = x; y < B; y++){
int len = N / y + 1;
T[x][y] = new int[len]; T[x][y][0] = 0;
for (int i = 1; i < len; i++)
T[x][y][i] = (T[x][y][i - 1] + 1ll * F[i] * G[i][x] % P * G[i][y] % P) % P;
}
}
void work(){
int n = read(),m = read(),ans = 0;
if (n > m) swap(n,m);
int E = m / B;
for (int i = 1; i <= E; i++)
ans = (ans + 1ll * F[i] * G[i][n / i] % P * G[i][m / i] % P) % P;
for (int i = E + 1,nxt; i <= n; i = nxt + 1){
nxt = min(n / (n / i),m / (m / i));
ans = ((ans + T[n / i][m / i][nxt]) % P + P - T[n / i][m / i][i - 1]) % P;
}
printf("%d
",ans);
}
int main(){
init();
int T = read();
while (T--) work();
return 0;
}