Link
反演的过程我就简略写写了。
记(n=max(A,B,C))。
(Ans=sumlimits_{i=1}^Asumlimits_{j=1}^Bsumlimits_{k=1}^Csumlimits_{u|i}sumlimits_{v|j}sumlimits_{w|k}epsilon((u,v))epsilon((u,w))epsilon((v,w)))
(Ans=sumlimits_{i=1}^Asumlimits_{j=1}^Bsumlimits_{k=1}^Clfloorfrac Ai
floorlfloorfrac Bj
floorlfloorfrac Ck
floorsumlimits_{u|(i,j)}mu(u)sumlimits_{v|(i,k)}mu(v)sumlimits_{w|(j,k)}mu(w))
(Ans=sumlimits_{u=1}^{min(A,B)}sumlimits_{j=1}^{min(A,C)}sumlimits_{k=1}^{min(B,C)}mu(u)mu(v)mu(w)sumlimits_{operatorname{lcm}(u,v)|i}lfloorfrac Ai
floorsumlimits_{operatorname{lcm}(u,w)|j}lfloorfrac Bj
floorsumlimits_{operatorname{lcm}(v,w)|k}lfloorfrac Ck
floor)
显然(f(x,yin{A,B,C})=sumlimits_{x|d}lfloorfrac yd
floor)可以预处理。
(Ans=sumlimits_{u=1}^{min(A,B)}sumlimits_{j=1}^{min(A,C)}sumlimits_{k=1}^{min(B,C)}mu(u)mu(v)mu(w)f(operatorname{lcm}(u,v),A)f(operatorname{lcm}(u,w),B)f(operatorname{lcm}(v,w),C))
(Ans=sumlimits_{u=1}^nsumlimits_{j=1}^nsumlimits_{k=1}^nmu(u)mu(v)mu(w)f(operatorname{lcm}(u,v),A)f(operatorname{lcm}(u,w),B)f(operatorname{lcm}(v,w),C))
观察到(x>yRightarrow f(x,y)=0)，记(n=max(A,B,C))，那么我们可以建一张无向图(G=([n],{(u,v)|mu(u)
e0wedgemu(v)
e0wedgeoperatorname{lcm}(u,v)le n}))，显然对答案有贡献的((u,v,w))在图中构成一个三元环（注意自环）。
打表发现(n=100000)时(|E|=760741)，那么用(O(|E|^{1.5}))的三元环计数就行了。
然后考虑如何优雅地建边，我们可以枚举(d=gcd(u,v))，再枚举(a=frac ud,b=frac vd)，因为(operatorname{lcm}(u,v)=dab<n)所以枚举的复杂度是(O(nlog^2n))的。
#include<bits/stdc++.h>
using namespace std;
#define LL long long
#define N 100007
#define P 1000000007
int read(){int x;scanf("%d",&x);return x;}
int gcd(int n,int m){return !m? n:gcd(m,n%m);}
int n,tot,pr[N>>3],mu[N],deg[N],flg[N];
LL fa[N],fb[N],fc[N],sa[N],sb[N],sc[N];
vector<pair<int,int> >G[N];
struct edge{int u,v,w;int cmp(){return deg[u]<deg[v]||(deg[u]==deg[v]&&u<v);}}e[N<<4];
void add(int u,int v,int w){G[u].push_back(make_pair(v,w));}
void Init()
{
flg[1]=mu[1]=1;
for(int i=2,j;i<=100000;++i)
{
if(!flg[i]) pr[++tot]=i,mu[i]=-1;
for(j=1;pr[j]*i<=100000&&j<=tot;++j)
{
flg[i*pr[j]]=1;
if(!(i%pr[j])) {mu[i*pr[j]]=0;break;}
mu[i*pr[j]]=-mu[i];
}
}
}
int main()
{
Init();
int n,m,T=read(),a,b,c,i,j,k,u,v,w,w1,w2,t;
LL ans;
for(;T;--T)
{
a=read(),b=read(),c=read(),n=max(a,max(b,c)),m=ans=0;
for(i=1;i<=n;++i) for(j=i;j<=n;j+=i) fa[i]+=a/j,fb[i]+=b/j,fc[i]+=c/j;
for(i=1;i<=n;++i) ans+=mu[i]*mu[i]*mu[i]*fa[i]*fb[i]*fc[i];
for(i=1;i<=n;++i)
for(j=1;j*i<=n;++j)
{
if(!mu[i*j]) continue;
for(k=1;k*i*j<=n;++k)
{
if(!mu[i*k]||k==j||gcd(j,k)!=1) continue;
u=i*j,v=i*k,w=i*j*k,t=mu[u]*mu[u]*mu[v],ans+=t*fa[u]*fb[w]*fc[w],ans+=t*fb[u]*fa[w]*fc[w],ans+=t*fc[u]*fb[w]*fa[w];
if(k<j) continue;
++m,e[m]=edge{u,v,w},++deg[u],++deg[v];
}
}
for(i=1;i<=m;++i)if(e[i].cmp()) add(e[i].u,e[i].v,e[i].w);else add(e[i].v,e[i].u,e[i].w);
for(u=1;u<=n;++u)
{
for(i=0;i<G[u].size();++i) v=G[u][i].first,w=G[u][i].second,sa[v]=fa[w],sb[v]=fb[w],sc[v]=fc[w];
for(i=0;i<G[u].size();++i)
for(v=G[u][i].first,w1=G[u][i].second,j=0;j<G[v].size();++j)
{
w=G[v][j].first,w2=G[v][j].second,t=mu[u]*mu[v]*mu[w];
ans+=t*fa[w1]*fb[w2]*sc[w];
ans+=t*fb[w1]*fa[w2]*sc[w];
ans+=t*fa[w1]*fc[w2]*sb[w];
ans+=t*fc[w1]*fa[w2]*sb[w];
ans+=t*fb[w1]*fc[w2]*sa[w];
ans+=t*fc[w1]*fb[w2]*sa[w];
}
for(i=0;i<G[u].size();++i) v=G[u][i].first,sa[v]=sb[v]=sc[v]=0;
}
for(i=1;i<=n;++i) fa[i]=fb[i]=fc[i]=deg[i]=0,G[i].clear();
printf("%d
",ans%P);
}
}