题目描述
设d(x)为x的约数个数,给定N、M,求∑i=1N∑j=1Md(ij)
输入输出格式
输入格式:输入文件包含多组测试数据。第一行,一个整数T,表示测试数据的组数。接下来的T行,每行两个整数N、M。
T行,每行一个整数,表示你所求的答案。
输入输出样例
2 7 4 5 6
110 121
说明
1<=N, M<=50000
1<=T<=50000
看到两个 求和符号,不假思索——莫比乌斯反演!
但是 d(ij) 怎么办呢?
没关系,我们有:
d(ij)=x∣i∑y∣j∑[gcd(x,y)=1]
我们所求的为Ans=i=1∑nj=1∑md(ij)
Ans=i=1∑nj=1∑mx∣i∑y∣j∑[gcd(x,y)=1]
看到这个[gcd(x,y)=1] ,我们就可以根据μ 的性质把它带进去:
Ans=i=1∑nj=1∑mx∣i∑y∣j∑d∣gcd(x,y)∑μ(d)
更换枚举项,由枚举gcd(x,y) 的约数,改为直接枚举d
Ans=i=1∑nj=1∑mx∣i∑y∣j∑d=1∑min(n,m)μ(d)∗[d∣gcd(x,y)]
将μ(d) 可以提出来,因为它与i,j 无关
Ans=d=1∑min(n,m)μ(d)i=1∑nj=1∑mx∣i∑y∣j∑[d∣gcd(x,y)]
接着,由枚举i,j 和它们的约数改变为,枚举它们的约数再直接乘上这些约数的倍数的个数。因为每一个约数都会对它的倍数产生贡献。
Ans=d=1∑min(n,m)μ(d)x=1∑ny=1∑m[d∣gcd(x,y)]⌊xn⌋⌊ym⌋
我们再一次更换枚举项,将枚举x,y 换为枚举dx,dy 。
这样[d∣gcd(x,y)] 这个条件就可以省去。
Ans=d=1∑min(n,m)μ(d)x=1∑⌊dn⌋y=1∑⌊dm⌋⌊dxn⌋⌊dym⌋
我们可以发现⌊dxn⌋ 与yy 无关,所以可以提前。
Ans=d=1∑min(n,m)μ(d)(x=1∑⌊dn⌋⌊dxn⌋)(y=1∑⌊dm⌋⌊dym⌋)
当我们将式子化简成这样的时候,我们已经可以看出,这个式子已经可以做到O(n) 计算了。
然后再使用 分块 的思想,每次询问O(sqrt(n)) 的时间复杂度即可。
附代码:
#include<iostream>
#include<algorithm>
#include<cstdio>
#define MAXN 50010
using namespace std;
int k=0,prime[MAXN],mu[MAXN],sum[MAXN];
long long f[MAXN];
bool np[MAXN];
inline int read(){
int date=0,w=1;char c=0;
while(c<'0'||c>'9'){if(c=='-')w=-1;c=getchar();}
while(c>='0'&&c<='9'){date=date*10+c-'0';c=getchar();}
return date*w;
}
void make(){
int m=MAXN-10;
mu[1]=1;
for(int i=2;i<=m;i++){
if(!np[i]){
prime[++k]=i;
mu[i]=-1;
}
for(int j=1;j<=k&&prime[j]*i<=m;j++){
np[prime[j]*i]=true;
if(i%prime[j]==0)break;
else mu[prime[j]*i]=-mu[i];
}
}
for(int i=1;i<=m;i++)sum[i]=sum[i-1]+mu[i];
for(int i=1;i<=m;i++){
long long s=0;
for(int j=1,last=1;j<=i;j=last+1){
last=i/(i/j);
s+=(long long)(last-j+1)*(i/j);
}
f[i]=s;
}
}
void work(){
int n,m;
long long ans=0;
n=read();m=read();
if(n>m)swap(n,m);
for(int i=1,last=1;i<=n;i=last+1){
last=min(n/(n/i),m/(m/i));
ans+=(long long)(sum[last]-sum[i-1])*f[n/i]*f[m/i];
}
printf("%lld\n",ans);
}
int main(){
int t=read();
make();
while(t--)work();
return 0;
}