【BZOJ5332】[SDOI2018]旧试题（数论，三元环计数）

【BZOJ5332】[SDOI2018]旧试题（数论，三元环计数）

题面

BZOJ
洛谷

题解

如果只有一个(sum)，那么我们可以枚举每个答案的出现次数。
首先约数个数这个东西很不爽，就搞一搞，变成(displaystyle sum_{d|i}1)
那么原式就可以写成：(displaystyle sum_{i=1}^Asum_{j=1}^Bsum_{k=1}^Csum_{d=1}^Ad|ijk)。
既然(d|ijk)，意味着(d)可以分别拆成(i)的一个因子，(j)的一个因子，(k)的一个因子的乘积，那么我们接着写，改写成(displaystyle sum_{i=1}^Asum_{j=1}^Bsum_{k=1}^Csum_{x|i}sum_{y|j}sum_{z|k}1)的形式。然而这个样子是会算重的。考虑这里怎么处理。

(displaystyle d(ij)=sum_{x|i}sum_{y|j}[gcd(x,y)=1])。
证明：
令(i*j=prod p_i^{a_i})，令(i=prod p_i^{b_i})，则(j=prod p_i^{a_i-b_i})。
令(x=prod p_i^{c_i},y=prod p_i^{d_i})，如果(gcd(x,y)=1)，那么必定有(c_i)和(d_i)中有一个取值为(0)，那么这个质因子的可以分配的方案数就是((b_i+1)+(a_i-b_i+1)-1=a_i+1)。
发现这个式子与唯一分解定理求解约数个数的式子一样，所以得证。
类似的可以证明(displaystyle d(ijk)=sum_{x|i}sum_{y|j}sum_{z|k}[(x,y)=1][(y,z)=1][(z,x)=1])。

这样子我们就可以把式子写成:$$displaystyle sum_{i=1}^Asum_{j=1}^Bsum_{k=1}^Csum_{x|i}sum_{y|j}sum_{z|k}[(x,y)=1][(y,z)=1][(z,x)=1]$$
然后改变枚举顺序：$$displaystyle sum_{x=1}^Asum_{y=1}^Bsum_{z=1}^C[(x,y)=1][(y,z)=1][(z,x)=1][frac{A}{x}][frac{B}{y}][frac{C}{z}]$$
然后把([(x,y)=1])变成(displaystyle sum_{d|x,d|y}mu(d))。
然后式子就变成了：

[displaystyle sum_{x=1}^Asum_{y=1}^Bsum_{z=1}^Csum_{d1|x,d1|y}mu(d1)sum_{d2|y,d2|z}mu(d2)sum_{d3|z,d3|x}mu(d3)[frac{A}{x}][frac{B}{y}][frac{C}{z}] ]

然后转而枚举(d1,d2,d3)：

[displaystyle sum_{d1=1}^Asum_{d2=1}^Bsum_{d3=1}^Cmu(d1)mu(d2)mu(d3)sum_{lcm(d1,d3)|x}[frac{A}{x}]sum_{lcm(d1,d2)|y}[frac{B}{y}]sum_{lcm(d2,d3)|z}[frac{C}{z}] ]

后面的东西，显然可以在调和级数的复杂度内预处理，并且当(lcm>max(A,B,C))的时候的值就是(0)。

那么考虑枚举两个数，如果它们两个的(lcmle max)的话就在他们之间连上一条边，这样子就得到了一张无向图，把自环处理掉，那么剩下的每一个三元环把他们对应到(d1,d2,d3)上面就是一组解。把自环丢掉的时候也丢掉了有两个数或者三个数相等的情况，所以还需要额外算一下。

然而直接枚举任意两个点是(O(n^2))的。
我们考虑枚举两者的(gcd)，然后再来枚举两个互质的数，以为(mu=0)也是没有贡献的，所以(mu=0)的数也不用考虑，这样子就可以减掉大量没有用的状态。
然后跑个三元环就好了。

#include<iostream>
#include<cstdio>
#include<algorithm>
#include<vector>
using namespace std;
#define ll long long
#define pb push_back
#define MAX 100100
#define MOD 1000000007
inline int read()
{
	int x=0;bool t=false;char ch=getchar();
	while((ch<'0'||ch>'9')&&ch!='-')ch=getchar();
	if(ch=='-')t=true,ch=getchar();
	while(ch<='9'&&ch>='0')x=x*10+ch-48,ch=getchar();
	return t?-x:x;
}
int pri[MAX],tot,mu[MAX];bool zs[MAX];
void Sieve()
{
	mu[1]=1;
	for(int i=2;i<MAX;++i)
	{
		if(!zs[i])pri[++tot]=i,mu[i]=-1;
		for(int j=1;j<=tot&&i*pri[j]<MAX;++j)
		{
			zs[i*pri[j]]=true;
			if(i%pri[j]==0)break;
			mu[i*pri[j]]=-mu[i];
		}
	}
}
int A,B,C,mx,dg[MAX],vis[MAX],book[MAX];
ll ans,fa[MAX],fb[MAX],fc[MAX];
vector<int> U,V,W;
vector<int> e[MAX],w[MAX];
int main()
{
	Sieve();int T=read();
	while(T--)
	{
		A=read();B=read();C=read();mx=max(A,max(B,C));ans=0;
		for(int i=1;i<=A;++i)for(int j=i;j<=A;j+=i)fa[i]+=A/j;
		for(int i=1;i<=B;++i)for(int j=i;j<=B;j+=i)fb[i]+=B/j;
		for(int i=1;i<=C;++i)for(int j=i;j<=C;j+=i)fc[i]+=C/j;
		for(int i=1;i<=A&&i<=B&&i<=C;++i)if(mu[i])ans+=mu[i]*fa[i]*fb[i]*fc[i];
		for(int i=1;i<=mx;++i)
		{
			if(!mu[i])continue;
			for(int j=1;i*j<=mx;++j)
			{
				if(!mu[i*j])continue;
				for(int k=j+1;1ll*i*j*k<=mx;++k)
				{
					if(!mu[i*k])continue;
					if(__gcd(j,k)>1)continue;
					int x=i*j,y=i*k,l=i*j*k;
					ans+=mu[y]*(fa[x]*fb[l]*fc[l]+fa[l]*fb[x]*fc[l]+fa[l]*fb[l]*fc[x]);
					ans+=mu[x]*(fa[y]*fb[l]*fc[l]+fa[l]*fb[y]*fc[l]+fa[l]*fb[l]*fc[y]);
					++dg[x];++dg[y];U.pb(x);V.pb(y);W.pb(l);
				}
			}
		}
		for(int i=0,l=U.size();i<l;++i)
		{
			int u=U[i],v=V[i];
			if(dg[u]<dg[v]||(dg[u]==dg[v]&&u<v))e[v].pb(u),w[v].pb(W[i]);
			else e[u].pb(v),w[u].pb(W[i]);
		}
		for(int i=1;i<=mx;++i)
		{
			for(int j=0,l=e[i].size();j<l;++j)vis[e[i][j]]=i,book[e[i][j]]=w[i][j];
			for(int j=0,l=e[i].size();j<l;++j)
			{
				int u=e[i][j];
				for(int k=0,lk=e[u].size();k<lk;++k)
				{
					int v=e[u][k];if(vis[v]!=i)continue;
					int a=w[i][j],b=book[v],c=w[u][k];
					ans+=mu[i]*mu[u]*mu[v]*(fa[a]*fb[b]*fc[c]+fa[a]*fb[c]*fc[b]+fa[b]*fb[a]*fc[c]+fa[b]*fb[c]*fc[a]+fa[c]*fb[a]*fc[b]+fa[c]*fb[b]*fc[a]);
				}
			}
		}
		printf("%lld
",ans%MOD);
		for(int i=1;i<=mx;++i)fa[i]=fb[i]=fc[i]=dg[i]=vis[i]=0;
		for(int i=1;i<=mx;++i)e[i].clear(),w[i].clear();
		U.clear();V.clear();W.clear();
	}
}