LOJ#6076. 「2017 山东一轮集训 Day6」三元组

Description

给定三个整数 $a,b,c$，一个三元组 $(i,j,k)$ 是合法的，当且仅当满足：

$i,j,k$ 均为整数
$1 leq i leq a,1 leq j leq b,1 leq k leq c$
$gcd(i,j)=gcd(i,k)=gcd(j,k)=1$
请求出合法的三元组数量对 $10^9+7$ 取模的值。

Solution

egin{align}
&sum _{i=1}^a sum _{j=1}^bsum _{k=1}^c [(i,j)=1][(j,k)=1][(i,k)=1]\
= & sum _{i=1}^a sum _{j=1}^bsum _{k=1}^c sum _{x|(i,j)}mu (x)sum _{y|(j,k)}mu (j)sum _{z|(i,k)}mu (z)\
= & sum _{x=1}sum _{y=1}sum _{z=1}mu (x)mu (y)mu (z)frac{a}{lcm(x,z)}frac{b}{lcm(x,y)}frac{c}{lcm(y,z)}
end{align}

发现值不为0的三元组数量有限，所以枚举所有值不为0的三元组，做三元环计数，特判有两个或者三个数相同的情况

#include<iostream>
#include<cstdio>
#include<cmath>
using namespace std;
int tot,prime[50005],lim,cnt,d[50005],head[50005],vis[50005];
const int mod=1e9+7;
long long ans,mu[50005],a,b,c;
bool vst[50005];
struct Node{
    int u,v,w;
}node[1000005];
struct Edge{
    int to,nxt,w;
}edge[1000005];
inline int read(){
    int f=1,w=0;
    char ch=0;
    while(ch<'0'||ch>'9'){if(ch=='-')f=-1;ch=getchar();}
    while(ch>='0'&&ch<='9')w=(w<<1)+(w<<3)+ch-'0',ch=getchar();
    return f*w;
}
int gcd(int x,int y){return y?gcd(y,x%y):x;}
inline long long calc1(long long x,long long y){return (a/x)*(b/y)*(c/y)+(a/y)*(b/x)*(c/y)+(a/y)*(b/y)*(c/x);}
inline long long calc2(long long x,long long y,long long z){return (a/x)*(b/y)*(c/z)+(a/x)*(b/z)*(c/y)+(a/y)*(b/x)*(c/z)+(a/y)*(b/z)*(c/x)+(a/z)*(b/x)*(c/y)+(a/z)*(b/y)*(c/x);}
inline void addedge(int u,int v,int w){edge[++tot]=(Edge){v,head[u],w},head[u]=tot;}
int main(){
    mu[1]=1;
    for(int i=2;i<=50000;i++){
        if(!vst[i])prime[++tot]=i,mu[i]=-1;
        for(int j=1;j<=tot&&i*prime[j]<=50000;j++){
            vst[i*prime[j]]=true;
            if(!(i%prime[j]))break;
            else mu[i*prime[j]]=-mu[i];
        }
    }
    a=read(),b=read(),c=read(),lim=max(a,max(b,c)),tot=0;
    for(int g=1;g<=lim;g++)if(mu[g])for(int i=1;i<=lim/g;i++)if(mu[i*g])for(int j=i+1;j<=lim/g/i;j++)if(mu[j*g]&&gcd(i,j)==1){
        int u=i*g,v=j*g,lcm=i*j*g;
        node[++cnt]=(Node){u,v,lcm},++d[u],++d[v],ans+=mu[u]*mu[u]*mu[v]*calc1(u,lcm)+mu[u]*mu[v]*mu[v]*calc1(v,lcm);
    }
    for(int i=1;i<=lim;i++)if(mu[i])ans+=mu[i]*(a/i)*(b/i)*(c/i);
    for(int i=1;i<=cnt;i++)d[node[i].u]>d[node[i].v]?addedge(node[i].u,node[i].v,node[i].w):addedge(node[i].v,node[i].u,node[i].w);
    for(int i=1;i<=lim;i++)if(mu[i]){
        for(int j=head[i];j;j=edge[j].nxt)vis[edge[j].to]=edge[j].w;
        for(int j=head[i];j;j=edge[j].nxt)if(mu[edge[j].to])for(int k=head[edge[j].to];k;k=edge[k].nxt)if(mu[edge[k].to]&&vis[edge[k].to])ans+=mu[i]*mu[edge[j].to]*mu[edge[k].to]*calc2(edge[j].w,edge[k].w,vis[edge[k].to]);
        for(int j=head[i];j;j=edge[j].nxt)vis[edge[j].to]=0;
    }
    printf("%lld
",(ans%mod+mod)%mod);
    return 0;
}