【LOJ #6076】「2017 山东一轮集训 Day6」三元组（莫比乌斯反演 / 三元环计数）

传送门

莫反得到$sum_{i=1}^asum_{j=1}^bsum_{k=1}^cfrac{a}{lcm(i,j)}frac b{lcm(j,k)}frac c{lcm(i,k)}$
然后按照旧试题那道题的做法三元环计数即可

#include<bits/stdc++.h>
using namespace std;
#define cs const
#define re register
#define pb push_back
#define pii pair<int,int>
#define ll long long
#define fi first
#define se second
#define bg begin
cs int RLEN=1<<20|1;
inline char gc(){
    static char ibuf[RLEN],*ib,*ob;
    (ib==ob)&&(ob=(ib=ibuf)+fread(ibuf,1,RLEN,stdin));
    return (ib==ob)?EOF:*ib++;
}
inline int read(){
    char ch=gc();
    int res=0;bool f=1;
    while(!isdigit(ch))f^=ch=='-',ch=gc();
    while(isdigit(ch))res=(res+(res<<2)<<1)+(ch^48),ch=gc();
    return f?res:-res;
}
template<class tp>inline void chemx(tp &a,tp b){a<b?a=b:0;}
template<class tp>inline void chemn(tp &a,tp b){a>b?a=b:0;}
cs int mod=1e9+7;
inline int add(int a,int b){return (a+=b)>=mod?(a-mod):a;}
inline int dec(int a,int b){a-=b;return a+(a>>31&mod);}
inline int mul(int a,int b){static ll r;r=1ll*a*b;return (r>=mod)?(r%mod):r;}
inline void Add(int &a,int b){(a+=b)>=mod?(a-=mod):0;}
inline void Dec(int &a,int b){a-=b,a+=a>>31&mod;}
inline void Mul(int &a,int b){static ll r;r=1ll*a*b;a=(r>=mod)?(r%mod):r;}
inline int ksm(int a,int b,int res=1){for(;b;b>>=1,Mul(a,a))(b&1)&&(Mul(res,a),1);return res;}
inline int Inv(int x){return ksm(x,mod-2);}
inline int fix(int x){return (x<0)?x+mod:x;}
inline int gcd(int a,int b){return b?gcd(b,a%b):a;}
cs int N=50005;
int a,b,c,lim,m,deg[N];
struct edge{
	int u,v,w;
	edge(int a=0,int b=0,int c=0):u(a),v(b),w(c){}
}E[N*20+1];
vector<pii> e[N];
int mu[N],pr[N],tot;
bitset<N>vis;
inline void init(cs int len=N-5){
	mu[1]=1;
	for(int i=2;i<=len;i++){
		if(!vis[i])pr[++tot]=i,mu[i]=-1;
		for(int j=1;j<=tot&&i*pr[j]<=len;j++){
			vis[i*pr[j]]=1;
			if(i%pr[j]==0)break;
			mu[i*pr[j]]=-mu[i];
		}
	}
}
inline ll calc1(int u,int lcm){
	return (1ll*(a/lcm)*(b/lcm)*(c/u)+1ll*(a/lcm)*(b/u)*(c/lcm)+1ll*(a/u)*(b/lcm)*(c/lcm));
}
inline ll calc2(int x,int y,int z){
	return (1ll*(a/x)*(b/y)*(c/z)+1ll*(a/x)*(b/z)*(c/y)+1ll*(a/y)*(b/x)*(c/z)+1ll*(a/y)*(b/z)*(c/x)+1ll*(a/z)*(b/x)*(c/y)+1ll*(a/z)*(b/y)*(c/x));
}
int vt[N];
int main(){
	#ifdef Stargazer
	freopen("lx.in","r",stdin);
	#endif
	init();
	a=read(),b=read(),c=read();
	lim=max(max(a,b),c);
	ll res=0;
	for(int g=1;g<=lim;g++)if(mu[g])
	for(int i=1,l=lim/g;i<=l;i++)if(mu[i*g])
	for(int j=i+1,l2=lim/g/i;j<=l2;j++)if(mu[j*g]&&gcd(i,j)==1){
		int u=i*g,v=j*g,lcm=i*j*g;
		deg[u]++,deg[v]++,E[++m]=edge(u,v,lcm);
		res+=calc1(u,lcm)*mu[u]*mu[v]*mu[u]+calc1(v,lcm)*mu[v]*mu[v]*mu[u];
	}
	for(int i=1;i<=lim;i++)res+=1ll*mu[i]*(a/i)*(b/i)*(c/i);
	for(int u,v,i=1;i<=m;i++){
		u=E[i].u,v=E[i].v;
		if(deg[u]>deg[v])e[u].pb(pii(v,E[i].w));
		else e[v].pb(pii(u,E[i].w));
	}
	for(int u=1;u<=lim;u++)if(mu[u]){
		for(pii &x:e[u])vt[x.fi]=x.se;
		for(pii &x:e[u])if(mu[x.fi]){
			int v=x.fi;
			for(pii &z:e[v])if(mu[z.fi]){
				int w=z.fi;if(!vt[w])continue;
				res+=calc2(vt[w],x.se,z.se)*mu[u]*mu[v]*mu[w];
			}
		}
		for(pii &x:e[u])vt[x.fi]=0;
	}
	cout<<(res%mod+mod)%mod<<'
';
}