2019.03.07【SDOI2018】【BZOJ5332】【洛谷P4619】旧试题(莫比乌斯反演)(三元环计数)

洛谷传送门

BZOJ传送门


解析:

很明显这是在致敬【SDOI2014】约数个数和。所以才叫"旧试题"

还是先化简式子。首先有二元组的结论,证明在上面那篇博客里面。

d ( i j ) = ∑ k ∣ i ∑ l ∣ j [ gcd ⁡ ( k , l ) = 1 ] d(ij)=\sum_{k\mid i}\sum_{l\mid j}[\gcd(k,l)=1] d(ij)=kilj[gcd(k,l)=1]

其实约数个数的结论可以推广到任意多元组,现在尝试推广到三元组。

d ( i j h ) = ∑ k ∣ i ∑ l t ∣ j h [ gcd ⁡ ( k , l t ) = 1 ] = ∑ k ∣ i ∑ l ∣ j ∑ t ∣ h [ gcd ⁡ ( k , l ) = 1 ] [ gcd ⁡ ( l , t ) = 1 ] [ gcd ⁡ ( k , t ) = 1 ] \begin{aligned} d(ijh)&=&&\sum_{k\mid i}\sum_{lt\mid jh}[\gcd(k,lt)=1]\\ &=&&\sum_{k\mid i}\sum_{l\mid j}\sum_{t\mid h}[\gcd(k,l)=1][\gcd(l,t)=1][\gcd(k,t)=1] \end{aligned} d(ijh)==kiltjh[gcd(k,lt)=1]kiljth[gcd(k,l)=1][gcd(l,t)=1][gcd(k,t)=1]

为了保证因子枚举的不重复性这里加上了 [ gcd ⁡ ( l , t ) = 1 ] [\gcd(l,t)=1] [gcd(l,t)=1]证明动动脑子都能想出来,就不写了

所以开始推式子:

A n s = ∑ i = 1 A ∑ j = 1 B ∑ h = 1 C d ( i j h ) = ∑ i = 1 A ∑ j = 1 B ∑ h = 1 C ∑ k ∣ i ∑ l ∣ j ∑ t ∣ h [ gcd ⁡ ( k , l ) = 1 ] [ gcd ⁡ ( l , t ) = 1 ] [ gcd ⁡ ( k , t ) = 1 ] = ∑ k = 1 A ∑ l = 1 B ∑ t = 1 C ⌊ A k ⌋ ⌊ B l ⌋ ⌊ C t ⌋ [ gcd ⁡ ( k , l ) = 1 ] [ gcd ⁡ ( l , t ) = 1 ] [ gcd ⁡ ( k , t ) = 1 ] = ∑ k = 1 A ∑ l = 1 B ∑ t = 1 C ⌊ A k ⌋ ⌊ B l ⌋ ⌊ C t ⌋ ∑ a ∣ k , a ∣ l μ ( a ) ∑ b ∣ l , b ∣ t μ ( b ) ∑ c ∣ k , c ∣ t μ ( c ) = ∑ a = 1 min ⁡ ( A , B ) ∑ b = 1 min ⁡ ( B , C ) ∑ c = 1 min ⁡ ( A , C ) μ ( a ) μ ( b ) μ ( c ) ∑ lcm ( a , c ) ∣ k ⌊ A k ⌋ ∑ lcm ( a , b ) ∣ l ⌊ B l ⌋ ∑ lcm ( b , c ) ∣ t ⌊ C t ⌋ \begin{aligned} Ans&=&&\sum_{i=1}^A\sum_{j=1}^B\sum_{h=1}^{C}d(ijh)\\ &=&&\sum_{i=1}^A\sum_{j=1}^B\sum_{h=1}^{C}\sum_{k\mid i}\sum_{l\mid j}\sum_{t\mid h}[\gcd(k,l)=1][\gcd(l,t)=1][\gcd(k,t)=1]\\ &=&&\sum_{k=1}^A\sum_{l=1}^B\sum_{t=1}^C\lfloor\frac{A}k\rfloor\lfloor\frac{B}{l}\rfloor\lfloor\frac{C}{t}\rfloor[\gcd(k,l)=1][\gcd(l,t)=1][\gcd(k,t)=1]\\ &=&&\sum_{k=1}^A\sum_{l=1}^B\sum_{t=1}^C\lfloor\frac{A}k\rfloor\lfloor\frac{B}{l}\rfloor\lfloor\frac{C}{t}\rfloor\sum_{a\mid k,a\mid l}\mu(a)\sum_{b\mid l,b\mid t}\mu(b)\sum_{c\mid k,c\mid t}\mu(c)\\ &=&&\sum_{a=1}^{\min(A,B)}\sum_{b=1}^{\min(B,C)}\sum_{c=1}^{\min(A,C)}\mu(a)\mu(b)\mu(c)\sum_{\text{lcm}(a,c)\mid k}\lfloor\frac{A}{k}\rfloor\sum_{\text{lcm}(a,b)\mid l}\lfloor\frac{B}{l}\rfloor\sum_{\text{lcm}(b,c)\mid t}\lfloor\frac{C}{t}\rfloor \end{aligned} Ans=====i=1Aj=1Bh=1Cd(ijh)i=1Aj=1Bh=1Ckiljth[gcd(k,l)=1][gcd(l,t)=1][gcd(k,t)=1]k=1Al=1Bt=1CkAlBtC[gcd(k,l)=1][gcd(l,t)=1][gcd(k,t)=1]k=1Al=1Bt=1CkAlBtCak,alμ(a)bl,btμ(b)ck,ctμ(c)a=1min(A,B)b=1min(B,C)c=1min(A,C)μ(a)μ(b)μ(c)lcm(a,c)kkAlcm(a,b)llBlcm(b,c)ttC

f ( n , k ) = ∑ k ∣ t ⌊ n t ⌋ f(n,k)=\sum_{k\mid t}\lfloor\frac{n}t\rfloor f(n,k)=kttn,显然,当 n n n确定的时候,可以 O ( n log ⁡ n ) O(n\log n) O(nlogn)时间内处理出 f n fn fn数组。

现在最蛋疼的是考虑前面的三个 μ \mu μ

这道题形式上就很反整除分块,我们换一个方向,考虑那些 μ \mu μ值不为 0 0 0,且在 f f f中有值的 a , b , c a,b,c a,b,c

如果 μ ( a ) ! = 0 , μ ( b ) ! = 0 \mu(a)!=0,\mu(b)!=0 μ(a)!=0,μ(b)!=0 f ( max ⁡ ( A , B , C ) , lcm ( a , b ) ) ! = 0 f(\max(A,B,C),\text{lcm}(a,b))!=0 f(max(A,B,C),lcm(a,b))!=0,那么我们在 a , b a,b a,b之间连一条边权为 lcm ( a , b ) \text{lcm}(a,b) lcm(a,b)边(包括自环,不过自环可以单独处理)。(其实这个边权就是为了后面来在 f f f数组里面找东西用的)

其实极限数据边数也只有700多。

现在发现了什么?

我们需要统计的所有 a , b , c a,b,c a,b,c三元组在原图中形成了一个三元环!

现在我们需要做的就是三元环计数了,直接做复杂度可能达到 O ( n m ) O(nm) O(nm),通过一些优化可以优化到 O ( m m ) O(m\sqrt m) O(mm ),具体的看这里:大佬博客

现在考虑一个问题,图怎么建? O ( n 2 ) O(n^2) O(n2)的暴力显然直接 T T T飞了。

我们考虑枚举 gcd ⁡ ( a , b ) , a gcd ⁡ ( a , b ) , b gcd ⁡ ( a , b ) \gcd(a,b),\frac{a}{\gcd(a,b)},\frac{b}{\gcd(a,b)} gcd(a,b),gcd(a,b)a,gcd(a,b)b,显然我们只需要保证 a , b a,b a,b μ \mu μ不为 0 0 0,且 a gcd ⁡ ( a , b ) , b gcd ⁡ ( a , b ) \frac{a}{\gcd(a,b)},\frac{b}{\gcd(a,b)} gcd(a,b)a,gcd(a,b)b互质就行了。

跑得还挺快的。


代码:

#include<bits/stdc++.h>
using namespace std;
#define ll long long
#define re register
#define cs const

cs int P=1e5+5;
int prime[P],pcnt;
bool mark[P];
int mu[P];

inline void linear_sieves(int len=P-5){
	mu[1]=1;
	for(int re i=2;i<=len;++i){
		if(!mark[i])prime[++pcnt]=i,mu[i]=-1;
		for(int re j=1;i*prime[j]<=len;++j){
			mark[i*prime[j]]=true;
			if(i%prime[j])mu[i*prime[j]]=-mu[i];
			else break;
		}
	}
}

inline int gcd(int a,int b){
	static int tmp;
	while(b){
		tmp=a%b;
		a=b;
		b=tmp;
	}
	return a;
}

cs int mod=1e9+7;
int n,A,B,C;
int fa[P],fb[P],fc[P];
int sa[P],sb[P],sc[P];
ll ans;

struct E{int to,val;};
vector<E> edge[P];

int deg[P];

struct data{
	int x,y,z;
}a[P*10];
int tot;

inline bool check(cs data& a){
	return deg[a.x]==deg[a.y]?a.x<a.y:deg[a.x]<deg[a.y];
}

inline void solve(){
	scanf("%d%d%d",&A,&B,&C);
	n=max(A,max(B,C));tot=ans=0;
	
	for(int re i=1;i<=n;++i)
	for(int re j=i;j<=n;j+=i)fa[i]+=A/j,fb[i]+=B/j,fc[i]+=C/j;
	
	for(int re i=1;i<=n;++i)ans+=(ll)mu[i]*fa[i]*fb[i]*fc[i];
	for(int re i=1;i<=n;++i)if(mu[i])
	for(int re j=1;j*i<=n;++j)if(mu[i*j])
	for(int re k=1;i*j*k<=n;++k)if(k!=j&&mu[i*k]&&gcd(k,j)==1){
		int x=i*j,y=i*k,z=i*j*k,tmp=mu[x]*mu[x]*mu[y];
		ans+=(ll)tmp*fa[x]*fb[z]*fc[z];
		ans+=(ll)tmp*fb[x]*fa[z]*fc[z];
		ans+=(ll)tmp*fc[x]*fa[z]*fb[z];
		if(x>y)a[++tot]=(data){x,y,z},++deg[x],++deg[y];
	}

	for(int re i=1;i<=tot;++i)
	if(check(a[i]))edge[a[i].x].push_back((E){a[i].y,a[i].z});
	else edge[a[i].y].push_back((E){a[i].x,a[i].z});
	
	for(int re x=1;x<=n;++x)if(mu[x]){
		for(int re e=0;e<edge[x].size();++e){
			int y=edge[x][e].to,z=edge[x][e].val;
			sa[y]=fa[z];
			sb[y]=fb[z];
			sc[y]=fc[z];
		}
		for(int re e=0,y,w1;e<edge[x].size();++e){
			y=edge[x][e].to,w1=edge[x][e].val;
			if(mu[y])for(int re k=0,z,w2,tmp;k<edge[y].size();++k){
				z=edge[y][k].to,w2=edge[y][k].val,tmp=mu[x]*mu[y]*mu[z];
				ans+=(ll)tmp*fa[w1]*fb[w2]*sc[z];
				ans+=(ll)tmp*fa[w1]*fc[w2]*sb[z];
				ans+=(ll)tmp*fb[w1]*fa[w2]*sc[z];
				ans+=(ll)tmp*fb[w1]*fc[w2]*sa[z];
				ans+=(ll)tmp*fc[w1]*fa[w2]*sb[z];
				ans+=(ll)tmp*fc[w1]*fb[w2]*sa[z];
			}
		}
		for(int re e=0;e<edge[x].size();++e){
			int y=edge[x][e].to;
			sa[y]=sb[y]=sc[y]=0;
		}
	}
	
	cout<<ans%mod<<"\n";
}

inline void init(){
	for(int re i=1;i<=n;++i)edge[i].clear();
	memset(fa+1,0,sizeof(int)*n);
	memset(fb+1,0,sizeof(int)*n);
	memset(fc+1,0,sizeof(int)*n);
	memset(deg+1,0,sizeof(int)*n);
}

int T;
signed main(){
	linear_sieves();
	scanf("%d",&T);
	while(T--){
		init();
		solve();
	}
	return 0;
}
评论
添加红包

请填写红包祝福语或标题

红包个数最小为10个

红包金额最低5元

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

抵扣说明:

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

余额充值