【BZOJ1101】【洛谷P3455 】【POI2007】ZAP-Queries(莫比乌斯反演)

传送门

题意:求满足x∈(1,m),y∈(1,n),且gcd(x,y)=dx\in(1,m),y\in(1,n),且gcd(x,y)=dx(1,m),y(1,n),gcd(x,y)=d的数的个数

考虑到我们要求

f(d)=∑i=1m∑j=1n[gcd(i,j)=d]f(d)=\sum_{i=1}^{m}\sum_{j=1}^{n} [gcd(i,j)=d]f(d)=i=1mj=1n[gcd(i,j)=d]

∑i=1⌊md⌋∑j=1⌊nd⌋[gcd(i,j)=1]\sum_{i=1}^{\lfloor \frac m d\rfloor}\sum_{j=1}^{\lfloor \frac n d \rfloor} [gcd(i,j)=1]i=1dmj=1dn[gcd(i,j)=1]
那么ans=f(d)=∑i=1m∑j=1n[gcd(i,j)=d]ans=f(d)=\sum_{i=1}^{m}\sum_{j=1}^{n} [gcd(i,j)=d]ans=f(d)=i=1mj=1n[gcd(i,j)=d]

=∑i=1⌊im⌋∑j=1⌊nd⌋[gcd(i,j)=1]=f(1)(n/=d,m/=d) =\sum_{i=1}^{\lfloor \frac i m\rfloor}\sum_{j=1}^{\lfloor \frac n d \rfloor} [gcd(i,j)=1]=f(1)(n/=d,m/=d) =i=1mij=1dn[gcd(i,j)=1]=f(1)(n/=d,m/=d)

我们发现求[gcd(i,j)=1][gcd(i,j)=1][gcd(i,j)=1]的复杂度是O(n2)O(n^2)O(n2)
怎么办?

考虑在[gcd(i,j)=d]=1[gcd(i,j)=d]=1[gcd(i,j)=d]=1时,我们发现肯定d∣i,d∣jd|i,d|jdi,dj
但是对于很多种ddd的倍数的i,ji,ji,j都只是gcd(i,j)=dp,p≥2,p∈Ngcd(i,j)=dp,p\geq 2,p\in\Ngcd(i,j)=dp,p2,pN

那我们考虑设F(x)=∑i=1m∑j=1n[gcd(i,j)为d的倍数]F(x)=\sum_{i=1}^{m}\sum_{j=1}^{n}[gcd(i,j)为d的倍数]F(x)=i=1mj=1n[gcd(i,j)d]
F(x)=∑i=1m∑j=1n[d∣gcd(i,j)]=⌊md⌋⌊nd⌋F(x)=\sum_{i=1}^{m}\sum_{j=1}^{n}[d|gcd(i,j)]=\lfloor \frac m d \rfloor \lfloor \frac n d \rfloorF(x)=i=1mj=1n[dgcd(i,j)]=dmdn
那么显然有F(d)=∑d∣nf(n)F(d)=\sum_{d|n} f(n)F(d)=dnf(n)

既然F(x)F(x)F(x)f(x)f(x)f(x)之间有规律而且F(x)F(x)F(x)有很好求
我们想到如果能求出f(x)f(x)f(x)F(x)F(x)F(x)的关系

这时就用到了莫比乌斯反演了

f(x)=∑x∣dμ(dx)F(d)=∑x∣dμ(dx)⌊md⌋⌊nd⌋f(x)=\sum_{x|d}\mu(\frac d x)F(d)=\sum_{x|d}\mu(\frac d x)\lfloor \frac m d \rfloor \lfloor \frac n d \rfloorf(x)=xdμ(xd)F(d)=xdμ(xd)dmdn

考虑到我们要求的是f(1)=∑i=1min(nd,md)μ(i)⌊mi⌋⌊ni⌋f(1)=\sum_{i=1}^{min(\frac n d,\frac m d)} \mu(i)\lfloor \frac m i \rfloor \lfloor \frac n i \rfloorf(1)=i=1min(dn,dm)μ(i)imin

这样一次回答的复杂度就是O(n)O(n)O(n)

考虑到多次询问,式子里有整除,可以用整除分块预处理μ\muμ前缀和O(n+m)O(\sqrt n+\sqrt m)O(n+m)求出答案

代码

#include<bits/stdc++.h>
using namespace std;
#define ll long long
inline int read(){
	char ch=getchar();
	int res=0,f=1;
	while(!isdigit(ch)){if(ch=='-')f=-f;ch=getchar();}
	while(isdigit(ch))res=(res+(res<<2)<<1)+(ch^48),ch=getchar();
	return res*f;
}
const int N=50005;
int vis[N],pr[N],mu[N],sum[N],tot;
inline void init(){
	mu[1]=1;
	for(int i=2;i<N;i++){
		if(!vis[i])pr[++tot]=i,mu[i]=-1;
		for(int j=1;j<=tot&&i*pr[j]<N;j++){
			vis[pr[j]*i]=1;
			if(i%pr[j]==0)break;
			mu[i*pr[j]]=-mu[i];
		}
	}
	for(int i=1;i<N;i++)sum[i]=sum[i-1]+mu[i];
}
signed main(){
	int T=read();init();
	for(int cas=1;cas<=T;cas++){
		int b=read(),d=read(),k=read();
        if(k==0){puts("0");continue;}
		ll ans=0;b/=k,d/=k;int p=min(b,d);
		for(int i=1,nxt;i<=p;i=nxt+1){
			nxt=min((b/(b/i)),(d/(d/i)));
			ans+=(1ll*(sum[nxt]-sum[i-1])*(b/i)*(d/i));
		}
		cout<<ans<<'\n';
	}
}

转载于:https://www.cnblogs.com/stargazer-cyk/p/11145679.html

评论
添加红包

请填写红包祝福语或标题

红包个数最小为10个

红包金额最低5元

当前余额3.43前往充值 >
需支付:10.00
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包
实付
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明:

1.余额是钱包充值的虚拟货币,按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载,可以购买VIP、付费专栏及课程。

余额充值