有一個結論:\(\varphi(nm)=\frac{\varphi(n)\varphi(m)\gcd(n,m)}{\varphi\left(\gcd(n,m)\right)}\),考慮證實:c++
化簡原式得:git
設 \(f(n)=\sum\limits_{d\mid n}\mu\left(\frac{n}{d}\right)\frac{d}{\varphi(d)},g(T,i)=\sum\limits_{j=1}^{i}\varphi(ij)\),篩出 \(\mu,\varphi\) 後便可預處理 \(f(n)\),發現 \(g(T,i)\) 中的 \(i\) 只有 \(\left\lfloor\frac{n}{T}\right\rfloor\) 種,因而其也能夠預處理。得答案爲:spa
考慮數論分塊,設 \(s(x,y,n)=\sum\limits_{i=1}^nf(i)g\left(x,i\right)g\left(y,i\right)\),那麼每次分塊加上的值就爲 \(s\left(\left\lfloor\frac{n}{l}\right\rfloor,\left\lfloor\frac{m}{l}\right\rfloor,r\right)-s\left(\left\lfloor\frac{n}{l}\right\rfloor,\left\lfloor\frac{m}{l}\right\rfloor,l-1\right)\)。code
發現對於比較大的 \(x,y\) 比較少,所以只預處理比較小的 \(x,y\),大的暴力作便可。設閾值爲 \(S\),得複雜度爲 \(O(n\log n +n S^2 +T(\left\lfloor\frac{n}{S}\right\rfloor+\sqrt n))\)get
#include<bits/stdc++.h> #define maxn 100010 #define maxs 35 #define all 100000 #define p 998244353 using namespace std; typedef long long ll; template<typename T> inline void read(T &x) { x=0;char c=getchar();bool flag=false; while(!isdigit(c)){if(c=='-')flag=true;c=getchar();} while(isdigit(c)){x=(x<<1)+(x<<3)+(c^48);c=getchar();} if(flag)x=-x; } int T,n,m,S=32,tot; int pri[maxn]; ll inv[maxn],mu[maxn],phi[maxn],f[maxn]; bool tag[maxn]; vector<ll> g[maxn],s[maxs][maxs]; int mod(int x) { return x>=p?x-p:x; } void init(int n) { inv[1]=phi[1]=mu[1]=1; for(int i=2;i<=n;++i) { inv[i]=(p-p/i)*inv[p%i]%p; if(!tag[i]) pri[++tot]=i,phi[i]=i-1,mu[i]=p-1; for(int j=1;j<=tot;++j) { int k=i*pri[j]; if(k>n) break; tag[k]=true; if(i%pri[j]) phi[k]=phi[i]*phi[pri[j]],mu[k]=p-mu[i]; else { phi[k]=phi[i]*pri[j]; break; } } } for(int i=1;i<=n;++i) for(int j=i;j<=n;j+=i) f[j]=mod(f[j]+mu[j/i]*i%p*inv[phi[i]]%p); for(int i=1;i<=n;++i) { g[i].resize(n/i+1),g[i][0]=0; for(int j=1;j<=n/i;++j) g[i][j]=mod(g[i][j-1]+phi[i*j]); } for(int i=1;i<=S;++i) { for(int j=1;j<=S;++j) { s[i][j].resize(min(n/i,n/j)+1),s[i][j][0]=1; for(int k=1;k<=min(n/i,n/j);++k) s[i][j][k]=mod(s[i][j][k-1]+f[k]*g[k][i]%p*g[k][j]%p); } } } ll solve(int n,int m) { ll v=0; if(n>m) swap(n,m); for(int i=1;i;++i) { if(n/i<=S&&m/i<=S) { for(int l=i,r;l<=n;l=r+1) { r=min(n/(n/l),m/(m/l)); v=mod(v+mod(s[n/l][m/l][r]-s[n/l][m/l][l-1]+p)); } return v; } v=mod(v+f[i]*g[i][n/i]%p*g[i][m/i]%p); } } int main() { init(all),read(T); while(T--) read(n),read(m),printf("%lld\n",solve(n,m)); return 0; }