【LOJ #6076】「2017 山东一轮集训 Day6」三元组（莫比乌斯反演 / 三元环计数）
莫反得到
然后按照旧试题那道题的做法三元环计数即可
#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<<'\n';
}
