【LOJ #2565】「SDOI2018」旧试题（莫比乌斯反演/三元环计数）

传送门

首先在$sdoi2015$约数个数和中有结论

$d(ij)=sum_{p|i}sum_{k|j}[gcd(i,j)==1]$

证明很简单，考虑每个质数$p,i=i'*p^a,j=j'*p^b$，在$d(ij)$中的贡献为$a+b+1$
对于后面那个式子的贡献也是$a+b+1$,所以是对的

对于$d(ijk)$类似的$sum_{p|i}sum_{t|j}sum_{l|k}[gcd(p,t)=1][gcd(t,l)=1][gcd(p,l)=1]$

于是可以莫反了

若设$f(x,y)=sum_{x|i}frac y i,L=max(A,B,C)$
则最后推出来就是
$sum_{i=1}^{min(A,B)}sum_{j=1}^{min(A,C)}sum_{k=1}^{min(B,C)}mu(i)mu(j)mu(k)f(lcm(i,j),A)f(lcm(i,k),B)f(lcm(j,k),C)$
循环上界改成$max(A,B,C)$显然是不会有影响的，因为$f(x,y)(x>y)=0$
又$f(x,y)=sum_{x|i}frac y i=sum_{i=1}^{frac y x}frac y{ix}=sum_{i=1}^{frac y x}frac{frac{y}{x}}{i}$
设$f(x)=sum_{i=1}^xfrac x i$
原式可以写作

$sum_{i=1}^{L}sum_{j=1}^{L}sum_{k=1}^{L}mu(i)mu(j)mu(k)f(frac{lcm(i,j)}A)f(frac{lcm(i,k)}B)frac{f(lcm(j,k)}C)$

考虑这东西仍然没法做
考虑把所有$mu$有值的数提出来看做点
如果两个数$i,j$满足$lcm(i,j)le L$就把他们之间连一条边
那么实际上要统计的就是所有三元环

然后你写一发发现$A,B,C=1e5$时加上自环边数是$821535$
于是可以愉快的写三元环计数了
把$i,j,k$中有相等的情况单独算一下

至于建边，可以先枚举$gcd$，然后再枚举$i,j$满足$gcd(i,j)=1$
这样复杂度差不多是$nlog^3n$，可能还要小一点，$n$大概是60000左右

然后就可以$O(nlog^3n+msqrt m)$做完了

#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<typename tp>inline void chemx(tp &a,tp b){a<b?a=b:0;}
template<typename 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);}
cs int N=100005;
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];
ll f[N];
int pr[N],tot,mu[N];
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];
		}
	}
	for(int i=1;i<=len;i++)
	for(int j=i;j<=len;j+=i)f[j]++;
	for(int i=1;i<=len;i++)f[i]+=f[i-1];
}
int vt[N],cnt,a,b,c,lim,deg[N];
inline void clear(){
	cnt=0;
	for(int i=1;i<=lim;i++)deg[i]=0,e[i].clear();
}
inline int gcd(int a,int b){
	return b?gcd(b,a%b):a;
}
inline ll calc1(int u,int lcm){
	return (f[a/lcm]*f[b/lcm]*f[c/u]+f[a/lcm]*f[b/u]*f[c/lcm]+f[a/u]*f[b/lcm]*f[c/lcm]);	
}
inline ll calc2(int l1,int l2,int l3){
	return (f[a/l1]*f[b/l2]*f[c/l3]+f[a/l1]*f[b/l3]*f[c/l2]+f[a/l2]*f[b/l1]*f[c/l3]+f[a/l2]*f[b/l3]*f[c/l1]+f[a/l3]*f[b/l1]*f[c/l2]+f[a/l3]*f[b/l2]*f[c/l1]);
}
inline void solve(){
	a=read(),b=read(),c=read();
	if(a>=b&&a>=c)lim=max(b,c);
	if(b>=a&&b>=c)lim=max(a,c);
	if(c>=a&&c>=b)lim=max(a,b);
	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[++cnt]=edge(u,v,lcm);
				res+=calc1(u,lcm)*mu[u]*mu[u]*mu[v]+calc1(v,lcm)*mu[v]*mu[v]*mu[u];
			}
		}
	}
	for(int i=1;i<=lim;i++)res+=mu[i]*f[a/i]*f[b/i]*f[c/i];
	for(int u,v,i=1;i<=cnt;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<<'
';
	clear();
}
int main(){
	#ifdef Stargazer
	freopen("lx.in","r",stdin);
	#endif
	init();
	int T=read();
	while(T--)solve();
}