公式太长了……我就直接抄一下这位大佬好了……实在懒得打了
首先据说$d(ij)$有个性质$$d(ij)=sum_{x|i}sum_{y|j}[gcd(x,y)=1]$$
我们所求的答案为$$ans=sum_{i=1}^{n}sum_{j=1}^{m}d(ij)$$
$$ans=sum_{i=1}^{n}sum_{j=1}^{m}sum_{x|i}sum_{y|j}[gcd(x,y)=1]$$
考虑一下$gcd(x,y)=1$,我们可以考虑莫比乌斯函数的性质,那么即$sum_{dmid n}mu(d)$与$[n=1]$的结果相同
则有$$ans=sum_{i=1}^{n}sum_{j=1}^{m}sum_{x|i}sum_{y|j}sum_{d|gcd(x,y)}mu(d)$$
然后我们由枚举$gcd(x,y)$的约数改为直接枚举$d$
$$ans=sum_{i=1}^{n}sum_{j=1}^{m}sum_{x|i}sum_{y|j}sum_{d=1}^{min(n,m)}mu(d)*[d|gcd(x,y)]$$
然后把$mu(d)$提取出来
$$ans=sum_{d=1}^{min(n,m)}mu(d)sum_{i=1}^{n}sum_{j=1}^{m}sum_{x|i}sum_{y|j}[d|gcd(x,y)]$$
然后,我们把枚举$i,j$和约数改为直接枚举约数,然后每个约数都会对他所有的倍数产生贡献
$$ans=sum_{d=1}^{min(n,m)}mu(d)sum_{x=1}^{n}sum_{y=1}^{m}[d|gcd(x,y)]lfloorfrac{n}{x} floorlfloorfrac{m}{y} floor$$
然后我们把枚举$x,y$改为枚举$dx,dy$,那么就可以把$[d|gcd(x,y)]$这个条件给消掉
$$ans=sum_{d=1}^{min(n,m)}mu(d)sum_{x=1}^{lfloorfrac{n}{d} floor}sum_{y=1}^{{lfloorfrac{m}{d} floor}}lfloorfrac{n}{dx} floorlfloorfrac{m}{dy} floor$$
$$ans=sum_{d=1}^{min(n,m)}mu(d)(sum_{x=1}^{lfloorfrac{n}{d} floor}lfloorfrac{n}{dx} floor)(sum_{y=1}^{{lfloorfrac{m}{d} floor}}lfloorfrac{m}{dy} floor)$$
然后$sum_{x=1}^{lfloorfrac{n}{d} floor}lfloorfrac{n}{dx} floor$和$sum_{y=1}^{{lfloorfrac{m}{d} floor}}lfloorfrac{m}{dy} floor$的前缀和都可以预处理,直接上整除分块就可以了
1 //minamoto 2 #include<iostream> 3 #include<cstdio> 4 #define ll long long 5 using namespace std; 6 #define getc() (p1==p2&&(p2=(p1=buf)+fread(buf,1,1<<21,stdin),p1==p2)?EOF:*p1++) 7 char buf[1<<21],*p1=buf,*p2=buf; 8 inline int read(){ 9 #define num ch-'0' 10 char ch;bool flag=0;int res; 11 while(!isdigit(ch=getc())) 12 (ch=='-')&&(flag=true); 13 for(res=num;isdigit(ch=getc());res=res*10+num); 14 (flag)&&(res=-res); 15 #undef num 16 return res; 17 } 18 const int N=50005; 19 int vis[N],p[N],mu[N],sum[N],m; 20 ll g[N],ans; 21 void init(int n){ 22 mu[1]=1; 23 for(int i=2;i<=n;++i){ 24 if(!vis[i]) p[++m]=i,mu[i]=-1; 25 for(int j=1;j<=m&&p[j]*i<=n;++j){ 26 vis[i*p[j]]=1; 27 if(i%p[j]==0) break; 28 mu[i*p[j]]=-mu[i]; 29 } 30 } 31 for(int i=1;i<=n;++i) sum[i]=sum[i-1]+mu[i]; 32 for(int i=1;i<=n;++i){ 33 ans=0; 34 for(int l=1,r;l<=i;l=r+1){ 35 r=(i/(i/l)); 36 ans+=1ll*(r-l+1)*(i/l); 37 } 38 g[i]=ans; 39 } 40 } 41 int main(){ 42 // freopen("testdata.in","r",stdin); 43 init(50000); 44 int n,m,T,lim;scanf("%d",&T); 45 while(T--){ 46 scanf("%d%d",&n,&m); 47 lim=min(n,m),ans=0; 48 for(int l=1,r;l<=lim;l=r+1){ 49 r=min(n/(n/l),m/(m/l)); 50 ans+=(sum[r]-sum[l-1])*g[n/l]*g[m/l]; 51 } 52 printf("%lld ",ans); 53 } 54 return 0; 55 }