首页 > 技术文章 > 三元环计数

REKonib 2022-01-20 23:33 原文

一个挺有意思的科技

模板题:洛谷P1989 无向图三元环计数

给定一个 \(n\) 个点 \(m\) 条边的简单无向图,求其三元环个数。
\(1\le n\le 10^5,1\le m\le 2\times10^5\)

考虑给所有边规定一个方向,以度数较小的点为起点。两点度数相等则以编号较小的点为起点。
这样建出来的新图显然是有向无环图
我们只需在新图中对形如 \(<u\to v>,<v\to w>,<u\to w>\) 的子图计数

枚举 \(u\) ,枚举 \(u\) 的出边 \(<u\to v>\) ,枚举 \(v\) 的出边 \(<v\to w>\)
在枚举 \(u\) 时给 \(u\) 的所有出边都打上标记,容易 \(O(1)\) 检验 \(u\)\(w\) 间是否有边
然后就过了,因为这个算法的时间复杂度是 \(O(m\sqrt m)\)
下面给出证明

新图上有一个性质:每个点的出度都不超过 \(\sqrt m\)
如果存在 \(d_u>\sqrt m\) ,注意到 \(u\) 的任意出边 \(<u\to v>\) 都满足 \(d_v\ge d_u>\sqrt m\)
这样的 \(v\) 共有 \(d_u\) 个,故 \(\sum\limits_{u\to v} d_v\ge d_u\sqrt m>m\)
\(\sum\limits_{v=1}^n d_v=m\) ,显然矛盾,故原结论得证

考虑每条边 \(<u\to v>\) ,它对复杂度的贡献为 \(d_v\) (枚举 \(v\) 的出边) 。
共有 \(m\) 条边,所以总贡献不超过 \(m\sqrt m\)
即时间复杂度为 \(O(m\sqrt m)\)

Click to view the code
#include<stdio.h>
const int N=100010,M=200010;
int n,m,x,y,d[N],last[N],v[N]; long long s;
struct node { int y,pre; }a[M],E[M];
int main() {
    scanf("%d%d",&n,&m);
    for (int i=1; i<=m; ++i)
        scanf("%d%d",&x,&y),
        ++d[x],++d[y],a[i]={x,y};
    for (int i=1; i<=m; ++i) { // 定向
        x=a[i].y,y=a[i].pre;
        if (d[x]<d[y]||(d[x]==d[y]&&x<y))
            E[i]={y,last[x]},last[x]=i; // x->y
        else E[i]={x,last[y]},last[y]=i; // y->x
    }
    for (int i=1; i<=n; ++i) {
        for (int j=last[i]; j; j=E[j].pre)
            v[E[j].y]=i; // 打标记
        for (int j=last[i]; j; j=E[j].pre) //
            for (int k=last[E[j].y]; k; k=E[k].pre)
                s+=(v[E[k].y]==i); // 计数
    }
    printf("%lld\n",s);
    return 0;
}

依然是模板题:洛谷P6815 [PA2009]Cakes

一个有 \(n\) 个点 \(m\) 条边的无向图,每个点有一个点权 \(a\)
对于一个三元环 \((i,j,k)(i<j<k)\),它的贡献是 \(\max(a_i,a_j,a_k)\)
求所有三元环的贡献和。

有点卡常
把链式前向星换成 vector ,这样内存访问连续一些,代码会变快
然后还可以加个快读
或者直接开 O2 也能过

Click to view the code
#include<stdio.h>
#include<vector>
#include<ctype.h>
#define pb push_back
const int N=100010,M=250010; long long ans;
int n,m,a[N],d[N],last[N],v[N],x[M],y[M];
std::vector<int>E[N];
inline int max(int x,int y) { return x>y?x:y; }
int read() {
    int x=0; char ch=getchar();
    while (!isdigit(ch)) ch=getchar();
    while (isdigit(ch)) x=x*10+(ch^48),ch=getchar();
    return x;
}

int main() {
    n=read(),m=read();
    for (int i=1; i<=n; ++i) a[i]=read();
    for (int i=1; i<=m; ++i)
        ++d[x[i]=read()],++d[y[i]=read()];
    for (int i=1; i<=m; ++i) {
        int X=x[i],Y=y[i];
        if (d[X]<d[Y]||(d[X]==d[Y]&&X<Y)) E[X].pb(Y);
        else E[Y].pb(X);
    }
    for (int i=1; i<=n; ++i) {
        for (auto j:E[i]) v[j]=i;
        for (auto j:E[i]) for (auto k:E[j])
            v[k]==i&&(ans+=max(max(a[i],a[j]),a[k]));
    }
    printf("%lld\n",ans);
    return 0;
}


给出一个 \(n\) 个点 \(m\) 条边的无向图,在这个图的 \(m\) 条边中任选 3 条边构成子图,求能得到连通子图的方案数。
\(3\le n\le 10^5,3\le m\le 3\times10^5\)

考虑如上三种情况即可

链形的可以枚举中间的边统计。树形的可以枚举中心点(指图中的 5 ),环形的就是跑一个三元环计数。

Click to view the code
#include<stdio.h>
const int N=100010,M=300010;
int n,m,d[N],v[N],last[N],x[M],y[M]; long long s,s2;
struct node { int y,pre; }E[M];
int main() {
    scanf("%d%d",&n,&m);
    for (int i=1; i<=m; ++i)
        scanf("%d%d",x+i,y+i),
        ++d[x[i]],++d[y[i]];
    for (int i=1; i<=m; ++i) {
        int tx=x[i],ty=y[i];
        s+=1ll*(d[tx]-1)*(d[ty]-1); 
        // 对链计数,并且每个三元环也会被数到 3 次
        // 因此后面要减去三元环个数的 2 倍来去重
        if (d[tx]<d[ty]||(d[tx]==d[ty]&&tx<ty))
            E[i]={ty,last[tx]},last[tx]=i;
        else E[i]={tx,last[ty]},last[ty]=i;
    }
    for (int i=1; i<=n; ++i) {
        s+=1ll*d[i]*(d[i]-1)*(d[i]-2)/6; // 对树计数
        for (int j=last[i]; j; j=E[j].pre) v[E[j].y]=i;
        for (int j=last[i]; j; j=E[j].pre)
            for (int k=last[E[j].y]; k; k=E[k].pre)
                s2+=(v[E[k].y]==i); // 对环计数
    }
    printf("%lld\n",s-(s2<<1));
    return 0;
}

[SDOI2018] 旧试题
给出 \(T\)\(A,B,C\) ,求出下面式子的值(对 \(10^9+7\) 取模)

\[\sum_{i=1}^A\sum_{j=1}^B\sum_{k=1}^Cd(ijk) \]

其中 \(d(ijk)\) 表示 \(i\times j\times k\) 的正因数个数
\(1\le T\le 10,1\le A,B,C\le 10^5,1\le \sum\max\{A,B,C\}\le 2\times10^5\)

数论 + 图论

对于乘积的约束个数我们有一个套路

(约定 \(x\perp y\) 表示 \(\gcd(x,y)=1\)

\[d(ijk)=\sum_{x\mid i}\sum_{y\mid j}\sum_{z\mid k}[x\perp y][y\perp z][z\perp x] \]

\(i,j,k\) 分解质因数即可证明,这里不再赘述

然后题目所求就是

\[\sum_{i=1}^A\sum_{j=1}^B\sum_{k=1}^C\sum_{x\mid i}\sum_{y\mid j}\sum_{z\mid k}[x\perp y][y\perp z][z\perp x] \]

推一下式子

\[\newcommand{\div}[2]{\left\lfloor\dfrac{#1}{#2}\right\rfloor} \sum_{x=1}^A\sum_{y=1}^B\sum_{z=1}^C[x\perp y][y\perp z][z\perp x]\div Ax \div By \div Cz\\ \\ \sum_{x=1}^A\div Ax\sum_{y=1}^B\div By\sum_{z=1}^C\div Cz\sum_{u\mid\gcd(x,y)}\mu(u)\sum_{v\mid\gcd(y,z)}\mu(v)\sum_{w\mid\gcd(z,x)}\mu(w)\\ \\ \sum_{u=1}^A\sum_{v=1}^B\sum_{w=1}^C\mu(u)\mu(v)\mu(w)\sum_{{\rm lcm}(w,u)\mid x}\div Ax\sum_{{\rm lcm}(u,v)\mid y}\div By\sum_{{\rm lcm}(v,w)\mid z}\div Cz\]

\(f(n,x)=\sum\limits_{x\mid i}\left\lfloor\dfrac ni\right\rfloor\) ,则上式即为

\[\sum_{u=1}^A\sum_{v=1}^B\sum_{w=1}^C\mu(u)\mu(v)\mu(w)f(A,{\rm lcm}(w,u))f(B,{\rm lcm}(u,v))f(C,{\rm lcm}(v,w)) \]

\(f\) 可以简单预处理一下,\(O(n\log n)\) 就够了

然而你发现推到这里依然只能 \(O(TABC)\) 爆算,啥用没有

我们考虑一个图论的建模,对任意两点 \(u,v\) ,连一条权值为 \({\rm lcm}(u,v)\) 的边

然后就变成了找三元环的问题(但要注意这题允许出现自环)

\(\mu(u)=0\) 的点和 \({\rm lcm}(u,v)>\max\{A,B,C\}\) 的边对答案的贡献都是 0 ,可以直接剪去

我们通过看题解写程序跑极限数据发现,处理之后的新图上的边数很少,\(O(m\sqrt m)\) 的算法可以被接受

但到这里还没做完,因为我们还没有考虑怎么建边

\(O(n^2)\) 暴力建边显然不可取,我们需要一个足够快的做法(最好是近似 \(O(m)\)

枚举 \(\rm lcm(u,v)\) 。根据前面的剪枝优化,可以令 \(1\le {\rm lcm}\le \max\{A,B,C\}\)\(\mu({\rm lcm})\ne 0\)
\(\mu\) 值不为 0 说明没有平方因子,容易想到将 \(\rm lcm\) 分解质因数
每个质因子都有三种方案:仅 \(u\) 包含,仅 \(v\) 包含,\(u\)\(v\) 均包含
然后枚举就完事了,反正连出的每条边都是未被剪去的有效边,时间复杂度差不多就是 \(O(m)\)

下面贴出代码

重边和自环好烦(),建议单独计算

Click to view the code
#include<stdio.h>
#include<string.h>
#include<vector>
#define clr(a) memset(a,0,sizeof a)
#define pb push_back
const int P=1e9+7,N=200010,M=1000010;
int T,A,B,C,n,m,cnt,ans,t,p[20000];
int v[N],mu[N],fA[N],fB[N],fC[N],d[N];
struct edge { int x,y,v; }E[M];
struct node { int y,v; }v2[N];
std::vector<node>a[N];
const double eps=1e-12; double inv[N];
inline int max(int x,int y) { return x>y?x:y; }
inline void add(int f,int x) { // 即 ans+=f*x
    if (f==1) (ans+=x)>=P&&(ans-=P);
    else (ans-=x)<0&&(ans+=P);
}

void calc_mu() { // 筛 mu ,顺便预处理一波单位分数,目的是用浮点数乘法优化整除
    mu[1]=1,inv[1]=1;
    for (int i=2; i<=n; ++i) {
        inv[i]=1./i+eps;
        v[i]||(p[++cnt]=i,mu[i]=-1);
        for (int j=1; j<=cnt&&(t=p[j]*i)<=n; ++j) {
            v[t]=1;
            if (i%p[j]==0) break;
            mu[t]=-mu[i];
        }
    }
}
void calc_f() { // 预处理 f
    for (int i=1; i<=n; ++i)
        for (int j=i; j<=n; j+=i) {
            double t=inv[j];
            fA[i]+=A*t,fB[i]+=B*t,fC[i]+=C*t; // 即加上 A/j,B/j,C/j
        }
}
void build(int x,int lcm,int u,int v) { // 建边
    if (u<v) E[++m]={u,v,lcm},++d[u],++d[v];
    for (int i=x+1; i<=cnt; ++i) {
        if (1ll*lcm*p[i]>n) break;
        int t=lcm*p[i];
        build(i,t,u*p[i],v);
        build(i,t,u,v*p[i]);
        build(i,t,u*p[i],v*p[i]);
    }
}

int count(int x,int y,int z) { // 计数
    int s=1ll*fA[z]*fB[x]*fC[y]%P;
    s=(1ll*fA[x]*fB[y]*fC[z]+s)%P;
    s=(1ll*fA[y]*fB[z]*fC[x]+s)%P;
    s=(1ll*fA[x]*fB[z]*fC[y]+s)%P,
    s=(1ll*fA[y]*fB[x]*fC[z]+s)%P,
    s=(1ll*fA[z]*fB[y]*fC[x]+s)%P;
    return s;
}
int main() {
    n=200000,calc_mu();
    scanf("%d",&T);
    while (T--) {
        scanf("%d%d%d",&A,&B,&C);
        n=max(A,max(B,C)),m=0,ans=0;
        clr(v2),clr(d),clr(fA),clr(fB),clr(fC),clr(a); // 多组数据,所有东西都要初始化一遍
        calc_f(),build(0,1,1,1);
        for (int i=1; i<=m; ++i) {
            int x=E[i].x,y=E[i].y,v=E[i].v;
            if (d[x]<d[y]||(d[x]==d[y]&&x<y)) a[x].pb({y,v});
            else a[y].pb({x,v});
            // 下面是计数(两个顶点重合的情况)
            t=1ll*fA[x]*fB[v]*fC[v]%P;
            t=(1ll*fA[v]*fB[x]*fC[v]+t)%P;
            t=(1ll*fA[v]*fB[v]*fC[x]+t)%P;
            add(mu[y],t);
            t=1ll*fA[y]*fB[v]*fC[v]%P;
            t=(1ll*fA[v]*fB[y]*fC[v]+t)%P;
            t=(1ll*fA[v]*fB[v]*fC[y]+t)%P;
            add(mu[x],t);
        }
        for (int i=1; i<=n; ++i) {
            if (!mu[i]) continue;
            add(mu[i],1ll*fA[i]*fB[i]*fC[i]%P); // 计数(三个顶点重合的情况)
            for (auto j:a[i]) v2[j.y]={i,j.v};
            for (auto j:a[i]) if (mu[j.y])
                for (auto k:a[j.y]) if (mu[k.y])
                    if (v2[k.y].y==i)
                        t=mu[i]*mu[j.y]*mu[k.y],
                        add(t,count(j.v,k.v,v2[k.y].v)); // 计数(三元环)
        }
        printf("%d\n",ans);
    }
    return 0;
}

推荐阅读