首先有一个等式是$\varphi(ab)=\frac{\varphi(a)\varphi(b)d}{\varphi(d)}$,其中$d=(a,b)$,这个比较好证,直接按展开式计算可得$\varphi(ab)\varphi(d)=\varphi(a)\varphi(b)d$(等号两边$d$的质因子都出现了两次)
然后推式子,得到$\sum\limits_{T=1}^nf(T)g\left(\left\lfloor\frac nT\right\rfloor,T\right)g\left(\left\lfloor\frac mT\right\rfloor,T\right)$,其中$f(n)=\sum\limits_{d|n}\mu\left(\frac nd\right)\frac d{\varphi(d)},g(n,k)=\sum\limits_{i=1}^n\varphi(ik)$
$f$可以枚举$d$更新倍数预处理,关键在于快速求$g$
我们先对每个$k$预处理$g\left(1\cdots\left\lfloor\frac nk\right\rfloor,k\right)$,这是$O(n\log n)$的
用分段的思想,设一个阈值$C$,对于$\leq\frac nC$的$T$用预处理好的$g$暴力求值,若$T\gt\frac nC$,这时$\left\lfloor\frac nT\right\rfloor\leq C$,我们再预处理$h_{i,j,k}=f(k)\sum\limits_{t=1}^i\varphi(tk)\sum\limits_{t=1}^j\varphi(tk)$,用$h$还有根号分段算剩下的部分
#include<stdio.h> typedef long long ll; const int mod=998244353,T=100000,C=30; void swap(int&a,int&b){a^=b^=a^=b;} int min(int a,int b){return a<b?a:b;} int mul(int a,int b){return a*(ll)b%mod;} int ad(int a,int b){return(a+b)%mod;} int de(int a,int b){return(a-b)%mod;} void inc(int&a,int b){(a+=b)%=mod;} int pow(int a,int b){ int s=1; while(b){ if(b&1)s=mul(s,a); a=mul(a,a); b>>=1; } return s; } int pr[T+10],phi[T+10],mu[T+10]; bool np[T+10]; void sieve(){ int i,j,M=0; phi[1]=1; mu[1]=1; for(i=2;i<=T;i++){ if(!np[i]){ pr[++M]=i; phi[i]=i-1; mu[i]=-1; } for(j=1;j<=M&&i*pr[j]<=T;j++){ np[i*pr[j]]=1; if(i%pr[j]==0){ phi[i*pr[j]]=phi[i]*pr[j]; break; } phi[i*pr[j]]=phi[i]*(pr[j]-1); mu[i*pr[j]]=-mu[i]; } } } int f[T+10],*g[T+10],*h[40][40]; int get(int n,int m){ if(n>m)swap(n,m); int i,nex,s; s=0; for(i=1;i<=min(T/C,n);i++)inc(s,mul(f[i],mul(g[i][n/i],g[i][m/i]))); for(;i<=n;i=nex+1){ nex=min(n/(n/i),m/(m/i)); inc(s,de(h[n/i][m/i][nex-T/C],h[n/i][m/i][i-1-T/C])); } return ad(s,mod); } int main(){ sieve(); int i,j,k,t,x,y,cas,n,m; for(i=1;i<=T;i++){ t=mul(i,pow(phi[i],mod-2)); for(j=i;j<=T;j+=i)inc(f[j],mul(mu[j/i],t)); } for(j=1;j<=T;j++){ g[j]=new int[T/j+1]; g[j][0]=0; for(i=1;i<=T/j;i++)g[j][i]=ad(g[j][i-1],phi[i*j]); } for(i=1;i<=C;i++){ for(j=1;j<=C;j++){ h[i][j]=new int[min(T/i,T/j)-T/C+2]; h[i][j][0]=0; } } for(k=T/C+1;k<=T;k++){ x=0; t=k-T/C; for(i=1;i<=T/k;i++){ inc(x,phi[i*k]); y=0; for(j=1;j<=T/k;j++){ inc(y,phi[j*k]); h[i][j][t]=ad(h[i][j][t-1],mul(mul(x,y),f[k])); } } } scanf("%d",&cas); while(cas--){ scanf("%d%d",&n,&m); printf("%d\n",get(n,m)); } }