思路
做这题先要知道一个性质,
[d_{ij}=sum_{x|i}sum_{y|j}[(x,y)=1]
]
然后上莫比乌斯反演颓柿子就好了
[egin{align}&sum_{i=1}^nsum_{j=1}^msum_{x|i}sum_{y|j}[(x,y)=1]\=&sum_{i=1}^nsum_{j=1}^msum_{x|i}sum_{y|j}sum_{d|(x,y)}mu(d)\=&sum_{i=1}^nsum_{j=1}^msum_d^{min(n,m)}mu(d)sum_{x=1}^{lfloorfrac{i}{d}
floor}sum_{y=1}^{lfloorfrac{j}{d}
floor}1\=&sum_d^{min(n,m)}mu(d)sum_{i=1}^nsum_{x=1}^{lfloorfrac{i}{d}
floor}sum_{j=1}^msum_{y=1}^{lfloorfrac{j}{d}
floor}1\=&sum_d^{min(n,m)}mu(d)sum_{x=1}^{lfloorfrac{n}{d}
floor}lfloorfrac{n}{xd}
floorsum_{y=1}^{lfloorfrac{m}{d}
floor}lfloorfrac{m}{dy}
floor\=&sum_d^{min(n,m)}mu(d)sum_{x=1}^{lfloorfrac{n}{d}
floor}lfloorfrac{lfloor frac{n}{d}
floor}{x}
floorsum_{y=1}^{lfloorfrac{m}{d}
floor}lfloorfrac{lfloor frac{m}{d}
floor}{y}
floorend{align}\
]
后面的部分,预处理一个(sum_{i=1}^n lfloor frac{n}{i} floor)就好了,前面上一个整除分块
代码
#include <cstring>
#include <cstdio>
#include <algorithm>
using namespace std;
int isprime[50010],mu[50010],cnt,iprime[50010],n,m,t;
long long sum[50010];
void prime(int n){
mu[1]=1;
isprime[1]=true;
for(int i=2;i<=n;i++){
if(!isprime[i]){
iprime[++cnt]=i;
mu[i]=-1;
}
for(int j=1;j<=cnt&&iprime[j]*i<=n;j++){
isprime[i*iprime[j]]=true;
mu[i*iprime[j]]=-mu[i];
if(i%iprime[j]==0){
mu[i*iprime[j]]=0;
break;
}
}
}
for(int i=1;i<=n;i++)
mu[i]+=mu[i-1];
}
void pre(int n){
for(int i=1;i<=n;i++){
long long ans=0;
for(int l=1,r;l<=i;l=r+1){
r=i/(i/l);
ans=ans+1LL*(r-l+1)*(i/l);
}
sum[i]=ans;
}
}
int main(){
scanf("%d",&t);
prime(50000);
pre(50000);
// printf("ok
");
while(t--){
scanf("%d %d",&n,&m);
long long ans=0;
for(int l=1,r;l<=min(n,m);l=r+1){
r=min(n/(n/l),m/(m/l));
// printf("l=%d
",l);
ans=ans+1LL*(mu[r]-mu[l-1])*sum[n/l]*sum[m/l];
}
printf("%lld
",ans);
}
return 0;
}