[HAOI2011]Problem b
解法
定义ans(N,M)表示$$sum_{i=1}Nsum_{j=1}M[gcd(i,j)=k]$$
容斥一下答案就是(ans(b,d)-ans(b,c-1)-ans(a-1,d)+ans(a-1,c-1))
我们再考虑单个(ans(N,M))怎么求
先把k提出来
[ans(N,M)=sum_{i=1}^frac{N}{k}sum_{j=1}^frac{M}{k}[gcd(i,j)=1]
]
然后考虑一下莫比乌斯反演
定义
[f(x)=sum_{i=1}^frac{N}{k}sum_{j=1}^frac{M}{k}[gcd(i,j)=x]
]
[g(x)=sum_{x|d}f(d)
]
那么
[g(x)=sum_{x|d}sum_{i=1}^frac{N}{k}sum_{j=1}^frac{M}{k}[gcd(i,j)=d]
]
然后把第一个(sigma)去掉
[g(x)=sum_{i=1}^frac{N}{k}sum_{j=1}^frac{M}{k}[gcd(i,j)\%x=0]
]
然后把(x)提出来去掉(gcd)的影响
[g(x)=sum_{i=1}^frac{N}{kx}sum_{j=1}^frac{M}{kx}[gcd(i,j)\%1=0]
]
显然
[g(x)=frac{N}{kx}frac{M}{kx}
]
根据莫比乌斯反演定理
[f(x)=sum_{x|d}μ(frac{d}{x})g(d)
]
[f(x)=sum_{x|d}μ(frac{d}{x})frac{N}{kd}frac{M}{kd}
]
而(ans(N,M)=f(1))
[ans(N,M)=sum_{d=1}^Nμ(d)frac{N}{kd}frac{M}{kd}
]
前缀和加数论分块(O(sqrt n))算出每一次的结果。
然后为了方便统计答案
我们发现当(kd>N)时，是没有贡献的
于是乎
[ans(N,M)=sum_{d=1}^frac{N}{k}μ(d)frac{N}{kd}frac{M}{kd}
]
所以如果先把(N,M)都除掉(k)那么
[ans(N,M)=sum_{d=1}^Nμ(d)frac{N}{d}frac{M}{d}
]
不要忘记容斥,这还不是最终答案。
完整代码
#include<bits/stdc++.h>
#define il inline
#define ll long long
#define rg register
using namespace std;
const int N=1e7;
ll pre[N],prime[N],isprime[N],s;
void work(){
for(rg int i=2;i<=50000;++i){
if(!isprime[i])pre[i]=-1,prime[++s]=i;
for(rg int j=1;j<=s&&i*prime[j]<=50000;++j){
isprime[i*prime[j]]=1;
pre[i*prime[j]]=-pre[i];
if(i%prime[j]==0){
pre[i*prime[j]]=0;
break;
}
}
}
pre[1]=1;
for(rg int i=1;i<=50000;++i)
pre[i]+=pre[i-1];
}
il ll js(int x,int y){
if(x>y)swap(x,y);
if(x==0)return 0;
ll ans=0;
for(rg int l=1,r;l<=x;l=r+1){
r=min(x/(x/l),y/(y/l));
ans+=(pre[r]-pre[l-1])*(x/l)*(y/l);
}
return ans;
}
int main(){
work();
int T;
cin>>T;
while(T--){
int a,b,c,d,k;
scanf("%d%d%d%d%d",&a,&b,&c,&d,&k);
b/=k,d/=k;
if(a%k)a=a/k+1;
else a/=k;
if(c%k)c=c/k+1;
else c/=k;
printf("%lld
",js(b,d)-js(b,c-1)-js(a-1,d)+js(a-1,c-1));
}
return 0;
}