【学习笔记】图上的环计数问题

· · 算法·理论

毕竟这也算是一个冷门的知识点了吧,所以没找到几道题 qwq

无向图三元环计数

题目:给定一张 n 个点 m 条边的无向图 G,问存在多少个三元组 (a,b,c) 满足:

前置定义:\text{deg}_u 表示 u 在图中的度数。

一般法

对于 n 个点的无向图 G 中全部的边 u\leftrightarrow v,考虑将其定向。

记新图 G'i 点的出度\text{deg}'_i

此时有一个十分优美的结论:

然后考虑开始计数。考虑对于原无向图 G 中一个三元环 (i,j,k),钦定 \text{deg}_i\le\text{deg}_j\le\text{deg}_k,其在图 G' 中必然有三条有向边 i\to j,i\to k,j\to k

根据上面的结论一,对于每一个点,其出度不会超过 \bf\sqrt{2m}。因此枚举三元环中任意一条边 i\to j,枚举 j 的出边所连接的点 k,只需要判断图中是否存在一条 i\to k 的边。可以把点对压缩为 int 类型然后哈希表判断,做到理论 O(m\sqrt m) 的时间复杂度。

做到严格 O(m\sqrt m):考虑对于枚举到的每一条有向边 u\to v,对于其中的结点 u,将其可以恰好一步走到的所有结点 x 做标记。对于此时枚举到的所有 v 可以一步走到的结点 w,若 w 已被标记,则 (u,v,w) 在图 G 中就是一个三元环。

上述的标记结点部分也可以通过设置时间戳来做。

坑点:

代码:

int deg[N], vis[N];
vector<int> z[N];
pair<int, int> edge[N];
void run() {
    int n = read(), m = read();
    for (int i = 1; i <= m; ++i) {
        int x = read(), y = read();
        ++deg[x], ++deg[y];
        edge[i] = {x, y};
    }
    for (int i = 1; i <= m; ++i) {
        int x = edge[i].first, y = edge[i].second;
        if (deg[x] < deg[y] || deg[x] == deg[y] && x < y) z[x].eb(y);
        else z[y].eb(x);
    }
    int cnt = 0;
    for (int i = 1; i <= n; ++i) {
        for (auto &j : z[i]) vis[j] = 1;
        for (auto &j : z[i]) for (auto &k : z[j])
            if (vis[k]) ++cnt; // (i, j, k) 是一个三元环 
        for (auto &j : z[i]) vis[j] = 0;
    }
    cout << cnt << '\n';
}

bitset 法

考虑对于每一个点 i,维护一个 bitset 表示 i 点一步能走到的点组成的集合。于是三元环计数可以枚举两个点然后用 bitset 求两个集合交中 1 的个数,做到 O(\frac{n^3}{\omega}),其中 \omega32

稠密图的情况下很好用,且比普通法好写。

代码(AT_abc258_g):

bitset<3010> f[3010];
void run() {
    int n = read();
    for (int i = 1; i <= n; ++i) {
        string s; cin >> s;
        s = ' ' + s;
        for (int j = i + 1; j <= n; ++j)
            if (s[j] == '1') f[i].set(j);
    }
    int cnt = 0;
    for (int i = 1; i <= n; ++i)
        for (int j = i + 1; j <= n; ++j)
            if (f[i][j])
                cnt += (f[i] & f[j]).count();
    cout << cnt << '\n';
}

upd 2024/11/26

这里曾经有一张图片,但是由于某些原因丢失了。

无向图四元环计数

题目:给定一张 n 个点 m 条边的无向图 G,问存在多少个四元组 (a,b,c,d) 满足:

特殊的,认为 (a,b,c,d)(a,c,b,d) 是两个不同的四元环。

还是沿用上面的做法,把所有点按照度数从小到大排序,然后按照上面的做法建立新图,并建立其的反图。 考虑枚举到图上的一个点 d,然后在其反图中找到一个其可以一次到达的结点 a。现在只需要找到有多少组合法的四元组 (a,\times,\times,d) 满足其形成了一个四元环。

考虑枚举 a 新图中所有可以一步到达的结点 b,若 b 可以在新图中一次走到 d,则将其放入点集 V 中。最终上述合法的四元环中两个 \times 只需要在点集 V 中任取两个点 u,v 即可。

和上面的分析相似,时间复杂度还是 O(m\sqrt m)

int deg[N], buc[N];
vector<int> z[N], zz[N];
// 其中 z 为新图, zz 为新图的反图
pair<int, int> edge[N];
void run() {
    int n = read(), m = read();
    for (int i = 1; i <= m; ++i) {
        int x = read(), y = read();
        edge[i] = {x, y};
        ++deg[x], ++deg[y];
        z[x].eb(y), z[y].eb(x);
    }
    for (int i = 1; i <= m; ++i) {
        int x = edge[i].first, y = edge[i].second;
        if (deg[x] < deg[y] || deg[x] == deg[y] && x < y) zz[y].eb(x);
        else zz[x].eb(y);
    }
    int cnt = 0;
    for (int i = 1; i <= n; ++i) { // i 是环的最后一个点
        for (auto &j : zz[i]) // j 是环的第一个点
            for (auto &k : z[j])
                if (deg[i] > deg[k] || deg[i] == deg[k] && i > k)
                    cnt += buc[k]++;
        for (auto &j : zz[i])
            for (auto &k : z[j])
                buc[k] = 0;
    }
    cout << cnt << '\n';
}

例题

\boldsymbol{1.\ CF\ Gym\ 102028L\ Connected\ Subgraphs}

(不是为什么我刚学图上计数就开这么折磨的题啊)

题意:给定一张 n 个点 m 条边的无向图,问有多少种不同的选取四条边的方法,使得这四条边生成的导出子图连通。

数据范围:4\le n\le 10^5,4\le m\le 2\times 10^5

题解:

考虑分类讨论。发现上述选取的方案若合法则必然生成

容易发现五种情况的并集恰好为全集,且两两的交为空。因此只需要分别计算五种情况的答案并求和即可。

(1) 四元环

直接套四元环计数的板子即可。

(2) 三元环外加一个点

考虑继续沿用计算三元环的做法。因为三元环的每一个结点均与另两个结点有连边,且做三元环计数的时候枚举了每一个合法的三元环,因此考虑对于每一个三元环 (a,b,c),若选择其中一个点并外接一条边,不同的方案数为 \text{deg}_a+\text{deg}_b+\text{deg}_c-6。将所有这样的贡献相加即可得到这部分的答案。

(3) 菊花(十字)

考虑枚举菊花的重心结点 u,则该点对答案的贡献为 \binom{\text{deg}_u}{4}。将所有这样的贡献求和即可得到这部分的答案。

(4) 四个点的链外加一个点

考虑先枚举该形态的两个顶点 x,y 满足 \text{deg}_x=2,\text{deg}_y=3。然后考虑找到 x 一步可以走到的结点和 y 一步可以走到的结点。此时对答案的贡献为 (\text{deg}_x-1)\times\binom{\text{deg}_y-1}{2}

但是这显然是错的,因为 xy 可能会有共同的可以走到的结点。于是考虑对答案容斥。考虑到 y 枚举的一步可以走到的两个结点互不相同,因此只有可能是 x 一步走到的结点和 y 一步走到的结点之间出现了重复,共有 2 种情况。然后可以发现若重复则变为“三元环外加一个点”的情况,因此减去两倍第二类情况的答案即可。

(5) 五个点的链

考虑枚举链的重心结点 x,于是只需要枚举和 x 有边相连的两个结点 y,z 满足 y\neq z。考虑到 y,z 两个点相邻的点已经被 x 所占用了一个,因此对答案的贡献就是 (\text{deg}_y-1)\times(\text{deg}_z-1)。但是同样的,可能会出现选择的点重复的情况,则同样考虑容斥。

y' 表示 y 向外连的点,z' 表示 z 向外连的点。

于是讨论完了这题的所有情况。

加强版:FJWC2019 子图,但是我没找到能提交这个题的地方 qwq

代码:

vector<int> z[N], zz[N];
int deg[N], n, m, vis[N];
pair<int, int> edge[N];
int visit[10];
int cycle_3() {
    if (visit[0] != -1) return visit[0];
    fill(deg + 1, deg + n + 1, 0ll);
    for (int i = 1; i <= m; ++i) {
        int x = edge[i].first, y = edge[i].second;
        ++deg[x], ++deg[y];
        edge[i] = {x, y};
    }
    for (int i = 1; i <= n; ++i) z[i].clear();
    for (int i = 1; i <= m; ++i) {
        int x = edge[i].first, y = edge[i].second;
        if (deg[x] < deg[y] || deg[x] == deg[y] && x < y) z[x].eb(y);
        else z[y].eb(x);
    }
    int cnt = 0;
    for (int i = 1; i <= n; ++i) {
        for (auto &j : z[i]) vis[j] = 1;
        for (auto &j : z[i]) for (auto &k : z[j])
            if (vis[k]) ++cnt;
        for (auto &j : z[i]) vis[j] = 0;
    }
    return visit[0] = cnt % mod;
}
int buc[N];
int cycle_4() {
    if (visit[1] != -1) return visit[1];
    fill(deg + 1, deg + n + 1, 0ll);
    for (int i = 1; i <= n; ++i) z[i].clear(), zz[i].clear();
    for (int i = 1; i <= m; ++i) {
        int x = edge[i].first, y = edge[i].second;
        ++deg[x], ++deg[y];
        z[x].eb(y), z[y].eb(x);
    }
    for (int i = 1; i <= m; ++i) {
        int x = edge[i].first, y = edge[i].second;
        if (deg[x] < deg[y] || deg[x] == deg[y] && x < y) zz[y].eb(x);
        else zz[x].eb(y);
    }
    int cnt = 0;
    for (int i = 1; i <= n; ++i) {
        for (auto &j : zz[i])
            for (auto &k : z[j])
                if (deg[i] > deg[k] || deg[i] == deg[k] && i > k)
                    cnt += buc[k]++, cnt %= mod;
        for (auto &j : zz[i])
            for (auto &k : z[j])
                buc[k] = 0;
    }
    return visit[1] = cnt % mod;
}
inline int binom4(int x) {
    return x * (x - 1) % mod * (x - 2) % mod * (x - 3) % mod * Inv(24, mod) % mod;
}
// 菊花(十字)
int count_flowers() {
    if (visit[2] != -1) return visit[2];
    fill(deg + 1, deg + n + 1, 0ll);
    for (int i = 1; i <= m; ++i) {
        int x = edge[i].first, y = edge[i].second;
        ++deg[x], ++deg[y];
    }
    int cnt = 0;
    for (int i = 1; i <= n; ++i)
        cnt += binom4(deg[i]);
    return visit[2] = cnt;
}
// 三元环外加一个点
int cycle_3_with_1() {
    if (visit[3] != -1) return visit[3];
    fill(deg + 1, deg + n + 1, 0ll);
    for (int i = 1; i <= m; ++i) {
        int x = edge[i].first, y = edge[i].second;
        ++deg[x], ++deg[y];
    }
    for (int i = 1; i <= n; ++i) z[i].clear();
    for (int i = 1; i <= m; ++i) {
        int x = edge[i].first, y = edge[i].second;
        if (deg[x] < deg[y] || deg[x] == deg[y] && x < y) z[x].eb(y);
        else z[y].eb(x);
    }
    int cnt = 0;
    for (int i = 1; i <= n; ++i) {
        for (auto &j : z[i]) vis[j] = 1;
        for (auto &j : z[i]) for (auto &k : z[j])
            if (vis[k]) cnt += (deg[i] + deg[j] + deg[k] - 6), cnt %= mod;
        for (auto &j : z[i]) vis[j] = 0;
    }
    return visit[3] = cnt;
}
int binom2(int x) {
    return x * (x - 1) % mod * Inv(2, mod) % mod;
}
// 四个点的链外加一个点
int link_4_with_1() {
    if (visit[4] != -1) return visit[4];
    fill(deg + 1, deg + n + 1, 0ll);
    for (int i = 1; i <= n; ++i) z[i].clear();
    for (int i = 1; i <= m; ++i) {
        int x = edge[i].first, y = edge[i].second;
        ++deg[x], ++deg[y];
        z[x].eb(y), z[y].eb(x);
    }
    int cnt = 0;
    for (int x = 1; x <= n; ++x)
        if (deg[x] >= 2)
            for (auto &y : z[x])
                if (x != y && deg[y] >= 3)
                    cnt += (deg[x] - 1) * binom2(deg[y] - 1), cnt %= mod;
    return visit[4] = (cnt + mod + mod - cycle_3_with_1() * 2) % mod;
}
// 五个点的链
int link_5() {
    if (visit[5] != -1) return visit[5];
    fill(deg + 1, deg + n + 1, 0ll);
    for (int i = 1; i <= m; ++i) {
        int x = edge[i].first, y = edge[i].second;
        ++deg[x], ++deg[y];
    }
    int cnt = 0;
    for (int i = 1; i <= n; ++i)
        if (deg[i] >= 2) {
            int pref_sum = 0;
            for (int i0 = 0; i0 < z[i].size(); ++i0) {
                int x = z[i][i0];
                if (deg[x] >= 2)
                    cnt = (cnt + (deg[x] - 1) * pref_sum % mod) % mod, pref_sum = (pref_sum + deg[x] - 1 + mod) % mod;
            }
        }
    return visit[5] = (cnt - 3 * cycle_3() - 4 * cycle_4() - 2 * cycle_3_with_1() + mod * 9) % mod;
}
void run() {
    // freopen("1.in", "r", stdin);
    int T = read(), ca = 1, tc = T, n1;
    while (T--) {
        cerr << (ca++) << '\n';
        memset(visit, -1, sizeof visit);
        n = read(), m = read();
        if (ca == 1) n1 = n;
        for (int i = 1; i <= m; ++i) {
            int x = read(), y = read();
            edge[i] = {x, y};
        }
        int x1 = cycle_4();
        int x2 = count_flowers();
        int x3 = link_5();
        int x4 = link_4_with_1();
        int x5 = cycle_3_with_1();
        cout << (x1 + x2 + x3 + x4 + x5) % mod << '\n';
    }
    cerr << clock() << '\n';
}

(甚至是最劣解)

\boldsymbol{2.\ HDU\ 6184\ Counting\ Stars}

题意:给定一张 n 个点 m 条边的无向图 G,问可以在其中找到多少个下面的图形:

数据范围:1\le n\le 3\times 10^5,\sum n\le 6\times 10^5

题解:

基础题。

考虑到计数三元环的时候是枚举了每一个三元环的,因此考虑先枚举每一个三元环,然后对每一条边记录有多少个三元环使用了该边。设有 a_i 个三元环使用了 i 边,则对答案的贡献是 \binom{a_i}{2}。将所有的贡献相加即可得到答案。

代码:

不知道为什么 MLE 了,我寻思我只用了 7MB 的内存啊/yiw

\boldsymbol{3.\ QOJ 6534,\ Um\_nik\ mod\ 998\ 244\ 353\ Contest\ K\ 4}

题意:给定一张 n 个点 m 条边的无向图 G,问其有多少个子图是 K_4

数据范围:4\le n\le 10^50\le m\le 10^5

题解:考虑枚举 K_4 中任意一个结点 i,问题变为计数有多少个三元环都和 i 结点有边相连。

然后把图上所有的边都定向,使用 bitset 对 i 的所有出度可以通往的结点做三元环计数即可。

时间复杂度为什么是对的?证明:考虑令 i 在定向后图的出度为 \text{deg}_i\le\sqrt{2m},视作 O(\sqrt m) 级别。然后对这 \text{deg}_i 个点做 bitset 三元环计数,时间复杂度为 O(\frac{{\text{deg}_i}^3}{\omega}),其中 \omega32

为了让 \text{deg}_i 尽量大,所以钦定 \text{deg}_i\approx O(\sqrt m)。又因为 \sum\limits_{i=1}^n\text{deg}_i=2m\approx O(m),因此此时不同的 i 的数量为 O(\sqrt m) 级别。

总时间复杂度即为 O(\sqrt m\times \frac{m^\frac{3}{2}}{\omega})=O(\frac{m^2}{\omega}),可以通过本题。

代码:

vector<int> z[N];
pair<int, int> edge[N];
bitset<1010> bit[100100];
int to[N];
int deg[N];
void run() {
    int n = read(), m = read();
    for (int i = 1; i <= m; ++i) {
        edge[i].first = read(), edge[i].second = read();
        ++deg[edge[i].first], ++deg[edge[i].second];
    }
    for (int i = 1; i <= m; ++i) {
        int x = edge[i].first, y = edge[i].second;
        if (deg[x] < deg[y] || deg[x] == deg[y] && x < y) z[x].eb(y);
        else z[y].eb(x);
    }
    memset(to, -1, sizeof to);
    int cnt = 0;
    for (int i = 1; i <= n; ++i) {
        for (int j = 0; j < z[i].size(); ++j)
            to[z[i][j]] = j;
        for (auto &j : z[i]) {
            bit[j].reset();
            for (auto &k : z[j])
                if (~to[k]) bit[j].set(to[k]);
        }
        for (auto &j : z[i])
            for (auto &k : z[j])
                if (~to[j] && ~to[k])
                    cnt += (bit[j] & bit[k]).count();
        for (int j = 0; j < z[i].size(); ++j)
            to[z[i][j]] = -1;
    }
    cout << cnt << '\n';
}

\boldsymbol{4.\ P3547\ [POI2013]\ CEN-Price\ List}

题意:给定一张 n 个点 m 条边的无向图 G,初始每一条边边权都为 a。然后对于任意两个点 x,y,若 x\to y 的最短路长度恰好为 2a,则再加一条 x\to y,边权为 b 的无向边。问在这张新的图中,k 点到其他每一个点的最短路长度是多少。

数据范围:n,m\le10^51\le a,b\le 1000。保证最后得到的新图连通。

题解:

容易发现为了让路径长度最短,有且仅可能有下面的三种情况:

首先 (1),(2) 两种情况可以一遍 dijkstra 或者直接广搜跑出来,问题在于 (3) 情况。

考虑枚举任意两条边 x\to yy\to z,可以在 x\to z 之间添加一条长度为 b 的边。但是这样是 O(m^2) 的,十分的不优秀。

考虑这样的一个结论:

结论:对于三个点 x,y,z,定义 \text{dis}_i 表示从起点到 i 点的最短路的长度,则若 x\to yy\to z 两条边边权的和大于或等于 x\to z 这条边边权,x,y,z 三个点不两两之间有边,则可以直接删去 y\to z 这条边。

感性理解容易发现结论正确。然后考虑到 x,y,z 三个点两两有边则 (x,y,z) 为图的一个三元环。考虑下一个性质:

证明:考虑三元环将边定向后得到的新有向图 G',每一个点的出度都不超过 \sqrt{2m}\approx O(\sqrt m)。考虑枚举三元环中的一条边,然后枚举该边中某一个顶点可以一步走到的顶点,这样的顶点数量级是不会超过 O(\sqrt m) 的。又因为新图也有 m 条有向边,所以合法的三元环数量为 O(m\sqrt m) 级别。

于是直接在松弛的时候删边,对三元环的情况特判。还剩下 O(m\sqrt m) 条边。对于这张图再跑最短路,因为边权相同所以可以直接广搜,时间复杂度为 O(m\sqrt m) 可以通过。

代码:

vector 存图想要删除某一个元素,若具体顺序不重要,则可以先把这个元素和末尾元素交换,然后再删除末尾元素,做到 O(1) 删除任意元素。

int dis[N], vis[N], dis2[N];
vector<int> z[N], zz[N];
void run() {
    int n = read(), m = read(), k = read(), a = read(), b = read();
    for (int i = 0; i < m; ++i) {
        int x = read(), y = read();
        z[x].eb(y), z[y].eb(x);
        zz[x].eb(y), zz[y].eb(x);
    }
    queue<int> q;
    memset(dis, 0x3f, sizeof dis);
    dis[k] = 0;
    q.push(k);
    while (q.size()) {
        int t = q.front();
        q.pop();
        for (auto &g : z[t])
            if (dis[g] > dis[t] + 1) {
                dis[g] = dis[t] + 1;
                q.push(g);
            }
    }
    memset(dis2, 0x3f, sizeof dis2);
    dis2[k] = 0;
    q.push(k);
    while (q.size()) {
        int t = q.front();
        q.pop();
        for (auto &g : z[t]) vis[g] = 1;
        for (auto &g : z[t])
            for (int i = 0; i < zz[g].size(); ++i) {
                int gg = zz[g][i];
                if (!vis[gg] && dis2[gg] > dis2[t] + b) {
                    dis2[gg] = dis2[t] + b;
                    q.push(gg);
                    swap(zz[g].back(), zz[g][i]);
                    zz[g].pop_back();
                    --i;
                }
            }
        for (auto &g : z[t]) vis[g] = 0;
    }
    for (int i = 1; i <= n; ++i)
        // cout << dis2[i] << '\n';
        cout << min({dis[i] * a, dis[i] / 2 * b + (dis[i] & 1) * a, dis2[i]}) << '\n';
} }

\boldsymbol{5.\ CF1468M\ Similar\ Sets}

题意:给定 n 个序列 a_i,同一序列中元素两两不同,问是否可以在其中找出两个序列满足它们至少有两个相同的元素。不需要计数

数据范围:1\le n\le 10^5,1\le a_{i,j}\le 2\times 10^9,\sum|a_i|\le 2\times 10^5

题解:

貌似是另一种套路?

考虑建图。先对序列离散化,然后将每一个序列中的元素向对应的序列连边。此时若找到任意一组四元环 (i,j,k,l),一定可以被循环移位直到该四元环为“元素-序列-元素-序列”的形式。

然后问题就解决了。时间复杂度为 O(m\sqrt m)。这里 m 表示建出来的图的边数。

代码:

int deg[N], id[N];
pair<int, int> edge[N];
vector<int> z[N], e[N], zz[N];
void run() {
    int T = read();
    while (T--) {
        int n = read();
        for (int i = 1; i <= n; ++i) {
            int oo = read();
            e[i].resize(oo);
            for (auto &j : e[i]) j = read();
        }
        int idx = n;
        set<int> se;
        map<int, int> mp;
        for (int i = 1; i <= n; ++i)
            for (auto &j : e[i]) se.insert(j);
        for (auto &x : se) mp[x] = ++idx;
        for (int i = 1; i <= n; ++i)
            for (auto &j : e[i]) j = mp[j];
        int idxx = 0;
        fill(deg + 1, deg + idx + 1, 0ll);
        for (int i = 1; i <= n; ++i)
            for (auto &j : e[i]) edge[++idxx] = {i, j}, ++deg[i], ++deg[j];
        for (int i = 1; i <= idx; ++i) z[i].clear(), zz[i].clear();
        for (int i = 1; i <= idxx; ++i) {
            auto [x, y] = edge[i];
            if (deg[x] < deg[y] || deg[x] == deg[y] && x < y) zz[y].eb(x);
            else zz[x].eb(y);
            z[x].eb(y), z[y].eb(x);
        }
        for (int i = 1; i <= idx; ++i) {
            for (auto &j : zz[i])
                for (auto &k : z[j])
                    if (deg[i] > deg[k] || deg[i] == deg[k] && i > k) {
                        if (!id[k]) id[k] = j + 99;
                        else {
                            // i j id[k] k
                            int t[4] = {i, j, id[k] - 99, k}; 
                            int lng = 0;
                            for (int d = 0; d < 4; ++d)
                                if (t[d] <= n) ++lng;
                            if (lng != 2) continue;
                            for (int d = 0; d < 4; ++d)
                                if (t[d] <= n) cout << t[d] << ' ';
                            cout << '\n';
                            for (auto &j : zz[i])
                                for (auto &k : z[j])
                                    id[k] = 0;
                            goto ee;
                        }
                    }
            for (auto &j : zz[i])
                for (auto &k : z[j])
                    id[k] = 0;
        }
        cout << "-1\n";
        ee:;
    }
} }

\boldsymbol{6.\ CF985G\ Team\ Players}

题意:

有一张 n 个点 m 条边的图,点的编号从 0 开始。求下列表达式的值:

\sum\limits_{i=0}^{n-1}\sum\limits_{j=i+1}^{n-1}\sum\limits_{k=j+1}^{n-1} [(i,j,k)三个点之间互相没有边相连](Ai+Bj+Ck)

其中 [expr] 为艾弗森括号,若 expr 这个布尔表达式为真则值为 1,否则为 0如:[T_TLucas_Yin 是可爱班花] 因括号内表达式为真所以答案为 1,其反命题的值则为 0(不是)

数据范围:3\le n\le 200000,0\le m\le 200000

题解:

首先如果要对两两之间有边的情况计数则相对简单,直接计数三元环然后对于每一个三元环计算其贡献。但是问题要对两两之间无边计数。建立反图可以解决这个问题但是反图边数是 O(n^2) 级别的,显然不行。

于是套路的考虑容斥答案。这样需要计数:

将上述部分的答案和容斥系数相乘后求和变为答案。问题是计算上述四个问题的答案。

(1) 所有三元组对答案的贡献

此时和图的形态无关,式子为:

\sum\limits_{i=0}^{n-1} \sum\limits_{j=i+1}^{n-1} \sum\limits_{k=j+1}^{n-1} (Ai+Bj+Ck)

考虑到最后的柿子 Ai+Bj+CkA,B,C 三个元素两两独立,因此考虑分开计算 A,B,C 对答案的贡献。

1. A 部分对答案的贡献

考虑抽离出和 A 有关的部分:

\sum\limits_{i=0}^{n-1} \sum\limits_{j=i+1}^{n-1} \sum\limits_{k=j+1}^{n-1} Ai

发现 j,k 的值都和最后答案的值没有关系,只需要考虑对于每一个 i,有多少种不同的合法的二元组 (j,k) 对答案存在贡献。容易发现 i 固定下上述合法的二元组数量是 \frac{(n-i+1)\times(n-i)}{2},因此实际的贡献为:

A\times \sum\limits_{i=0}^{n-1}(i\times\frac{(n-i+1)\times(n-i)}{2})

可以直接计数。

2. B 部分对答案的贡献

考虑抽离出和 B 有关的部分:

\sum\limits_{i=0}^{n-1}\sum\limits_{j=i+1}^{n-1}\sum\limits_{k=j+1}^{n-1}Bj

发现该表达式和 i,j 的值均有关系。因此考虑经典套路交换求和顺序,得到

\sum\limits_{j=0}^{n-1}\sum\limits_{i=0}^{j-1}\sum\limits_{k=j+1}^{n-1} Bj

此时式子只和 j 的值强相关,因此考虑固定 j 计数合法的二元组 (i,k) 的数量。容易发现答案为

B\times \sum\limits_{j=0}^{n-1}(j\times j\times (n-j-1))=B\times\sum\limits_{j=0}^{n-1}(j\times (n-j-1))

也可以直接计数。

3. C 部分对答案的贡献

B 部分相似的,交换求和顺序然后计数 (i,j) 二元组对答案的贡献。答案为

C\times \sum\limits_{k=0}^{n-1}\frac{k\times(k-1)}{2}

还是可以直接计数。

(2) 三元环对答案的贡献

考虑到做三元环计数的时候其实是枚举到了每一个三元环的,因此直接对于每一个三元环计算其对答案的贡献即可。

(3) 所有三元组满足其中有一对点对有边另两对点对无边对答案的贡献

考虑枚举该三元组中的一条边 u\to v,这里钦定 u<v。(显然 u\neq v。)以及三元组为 (x,y,z)1\le x<y<z\le n)。考虑分讨环上另一个元素 w 的下标对答案的影响:

只需要将上述七种情况的答案求和即可。

(4) 所有三元组满足其中有两对点对有边另一对点对无边对答案的贡献

最困难的一集

考虑枚举在这个大小为 3 的生成子图中度数为 2 的点 u,以及其可以一步走到的点 v,另一个没有提到的结点为 w。这里令 u 可以一步到达的结点的数目为 cv 点在其中排名为 pu\to v 这条边对答案的贡献可以分为 u 对答案的贡献和 v 对答案的贡献。

一、u 对答案的贡献

钦定三元组为 (x,y,z) 且满足 1\le x<y<z\le n,则:

将上述三种情况的贡献求和即可得到 u 对答案的贡献。

二、v 对答案的贡献

考虑对 u,v 的大小关系进行分类讨论。容易发现 u\neq v

+ $w<v$,则 $v$ 对答案的贡献为 $Cv(p-2)$。 + $w>v$,则 $v$ 对答案的贡献为 $Bv(c-p)$。 + $w=v$。容易发现此情况不存在。 否则即 $u>v$。同样需要继续分类讨论 $w$ 和 $v$ 的大小关系: + $w<v$,则 $v$ 对答案的贡献为 $Bv(p-1)$。 + $w>v$,则 $v$ 对答案的贡献为 $Av(c-i-1)$。 + $w=v$。容易发现此情况不存在。 于是将上述几种情况相加就得到了该类情况对答案的贡献。 于是这个 duliu 题就被做完了。 代码: 题解区里代码对于最后一种情况的处理十分简洁,于是学习了一下。 [![pA44foR.png](https://s21.ax1x.com/2024/11/27/pA44foR.png)](https://imgse.com/i/pA44foR) 上为使用数组,下为是使用 vector。警示后人,不要乱用 vector。 upd:调过了,更新一下代码。 ```cpp pair<int, int> edge[N]; vector<int> z[N], zz[N]; int a[N], n, m, A, B, C, deg[N]; using ull = unsigned long long; ull count_all() { ull s = 0; for (int i = 1; i <= n; ++i) s += A * (n - i) * (n - i - 1) / 2 * (i - 1) + B * (i - 1) * (i - 1) * (n - i) + C * (i - 1) * (i - 2) / 2 * (i - 1); return s; } int vis[N]; ull cycle_3() { ull s = 0; for (int i = 1; i <= n; ++i) { for (auto &j : zz[i]) vis[j] = 1; for (auto &j : zz[i]) for (auto &k : zz[j]) if (vis[k]) { int vec[3] = {i, j, k}; sort(vec, vec + 3); // vector<int> vec = {i, j, k}; sort(vec.begin(), vec.end()); s += A * vec[0] + B * vec[1] + C * vec[2] - A - B - C; } for (auto &j : zz[i]) vis[j] = 0; } return s; } ull cycle_1() { ull s = 0; for (int i = 1; i <= n; ++i) for (auto &j : z[i]) { int x = i, y = j; if (x > y) continue; s += A * ((x - 1) * (n - x - 1) + (x - 1) * (x - 2) / 2) + B * ((x - 1) * (x - 1) + (y - 1) * (n - y) + (y + x - 2) * (y - x - 1) / 2) + C * ((y - 1) * (y - 2) + (n + y - 1) * (n - y) / 2); } return s; } ull cycle_2() { ull s = 0; for (int i = 1; i <= n; ++i) { int allsize = z[i].size(); z[i].eb(i); sort(z[i].begin(), z[i].end()); for (int x = 0; x <= allsize; ++x) { int j = z[i][x]; if (i > j) s += A * (j - 1) * (allsize - x - 1) + B * (j - 1) * x; else if (i < j) s += B * (j - 1) * (allsize - x) + C * (j - 1) * (x - 1); else s += A * (allsize - x) * (allsize - x - 1) / 2 * (i - 1) + B * x * (allsize - x) * (i - 1) + C * x * (x - 1) / 2 * (i - 1); } } return s; } void run() { n = read(), m = read(); A = read(), B = read(), C = read(); for (int i = 1; i <= m; ++i) { edge[i].first = read() + 1, edge[i].second = read() + 1; z[edge[i].first].eb(edge[i].second), z[edge[i].second].eb(edge[i].first); ++deg[edge[i].first], ++deg[edge[i].second]; } for (int i = 1; i <= m; ++i) { auto [x, y] = edge[i]; if (deg[x] < deg[y] || deg[x] == deg[y] && x < y) zz[x].eb(y); else zz[y].eb(x); } ull all = count_all(), c1 = cycle_1(), c2 = cycle_2(), c3 = cycle_3(); // cout << "qwq " << all << ' ' << c1 << ' ' << c2 << ' ' << c3 << '\n'; cout << (all - c1 + c2 - c3) << '\n'; } ``` ### $\boldsymbol{7.\ CF\ Gym\ 104651C\ Clique\ Challenge}

题意:给定一个 n 个点 m 条边的无向图,问该图团的数量。

数据范围:1\le n\le 10001\le m\le 1000

题解:

仍然考虑求三元环时的定向方法,让度数小的点连向度数大的点,形成一张新的有向图。此时每一个点的出度不超过 \sqrt{2m}

然后考虑枚举点 u,并钦定点 u 是团中的点,且在新图中不存在任何点能够到达 u。根据团的性质此时其余的点都必然可以从 u 开始在新图中通过恰好一条边到达。又因为新图的性质,满足这样条件的点不会超过 \sqrt{2m}<45 个。考虑对这些点进行折半求出可以组成的团的数目,于是这个题就做完了。

时间复杂度证明:显然为了让时间复杂度最大,让新图中点的出度的最大值为 \sqrt{2m} 最优。此时最多可以得到 O(\sqrt m) 个这样的点。又因为折半搜索时间复杂度为 O(2^{\frac{\sqrt{2m}}{2}}),因此总的时间复杂度为 O(\sqrt m\times 2^\frac{\sqrt{2m}}{2}) 可以通过。

代码:

\boldsymbol{8.\ CF11D\ A\ Simple\ Task}

题意:给定一张 n 个点 m 条边的无向图,求该图中简单环的数量。

数据范围:1\le n\le 19m\ge 0,图无重边自环。

题解:

这个好像没有什么特别优秀的解法……考虑到 n\le 192^{19}<10^6,看上去十分的状压。

于是设 f_{i,j} 表示当前选择了 i 集合内的点,当前位于 j 点的方案数。这里钦定该环的起点为 \text{lowbit}(i)。转移是容易的。于是这个问题就做完了。

写代码的话需要注意一下一个环会被正着计数一遍然后再被反着计数一遍,答案需要除以 2

代码:

int n, m;
int mp[50][50], f[1ll << 20][20];
void add(int &x, int y) {
    x += y;
}
void run() {
    n = read(), m = read();
    while (m--) {
        int x = read() - 1, y = read() - 1;
        mp[x][y] = mp[y][x] = 1;
    }
    for (int i = 0; i < n; ++i) f[1ll << i][i] = 1;
    for (int i = 1; i < (1ll << n) - 1; ++i)
        for (int j = 0; j < n; ++j) if (i >> j & 1)
            for (int k = __lg(i &- i) + 1; k < n; ++k) if (~i >> k & 1)
                if (mp[j][k]) add(f[i | (1ll << k)][k], f[i][j]);
    int s = 0;
    for (int i = 1; i < (1ll << n); ++i)
        if (__builtin_popcountll(i) >= 3)
            for (int j = 0; j < n; ++j)
                if (mp[__lg(i &- i)][j])
                    s += f[i][j];
    cout << s / 2 << '\n';
} }

\boldsymbol{9.\ [PA2009]\ Cakes}

题意:给定一张 n 个点 m 条边的无向图,点有点权。定义三元组 (i,j,k) 的贡献是 \max(a_i,a_j,a_k)。问所有三元环 (i,j,k) 满足 i<j<k 对答案的贡献的和是多少。

数据范围:n<10^5,m<2.5\times 10^5

题解:

板子题。考虑到做三元环计数的时候具体枚举到了每一个三元环,所以直接对每一个三元环处理其贡献即可。时间复杂度为 O(m\sqrt m)

代码:

vector<int> z[N];
pair<int, int> edge[N];
int n, m, a[N], deg[N], vis[N];
void run() {
    int n = read(), m = read();
    for (int i = 1; i <= n; ++i) a[i] = read();
    for (int i = 0; i < m; ++i) {
        int x = read(), y = read();
        edge[i] = {x, y};
        ++deg[x], ++deg[y];
    }
    for (int i = 0; i < m; ++i) {
        auto [x, y] = edge[i];
        if (deg[x] < deg[y] || deg[x] == deg[y] && x < y) z[x].eb(y);
        else z[y].eb(x);
    }
    int score = 0;
    for (int i = 1; i <= n; ++i) {
        for (auto &j : z[i]) vis[j] = 1;
        for (auto &j : z[i]) for (auto &k : z[j])
            if (vis[k]) score += max({a[i], a[j], a[k]});
        for (auto &j : z[i]) vis[j] = 0;
    }
    cout << score << '\n';
} }

\boldsymbol{10.\ CF\ Gym\ 100342J\ Triatrip}

题意:给定一张有向稠密图 Gn 个点 m 条边,问图中有多少个三元环。

数据范围:3\le n\le 1500m\le n^2-n

题解:有向图无法使用给边定向,但是发现 n 十分小,因此考虑 bitset 做法。

记录第一个 bitset:F_i 表示 i 可以一步到达的结点组成的集合。

记录第二个 bitset:G_i 表示可以一步到达 i 的结点组成的集合。

然后考虑枚举三元环中的两个结点 i,j,钦定边的方向为 i\to j。只需要统计 k 的数量满足 i 可以被 k 一步到达,且 j 可以一步到达 k。答案即为 G_iF_j 两个 bitset 的交集的大小。直接这样做时间复杂度为 O(\frac{n^3}{\omega}),可以通过该题。

注意点:

  1. bitset 有向图会把三元环重复计算 3 遍,而无向图中是 6 遍。
  2. 需要开文件读写! 输入文件名为 triatrip.in,输出文件名为 triatrip.out。

代码:

char s[2010];
bitset<1510> f[1510], g[1510];
void run() {
    freopen("triatrip.in", "r", stdin);
    freopen("triatrip.out", "w", stdout);
    int n = read();
    for (int i = 1; i <= n; ++i) {
        scanf("%s", s + 1);
        for (int j = 1; j <= n; ++j)
            if (s[j] == '+') f[i].set(j), g[j].set(i);
    }
    int cnt = 0;
    for (int i = 1; i <= n; ++i)
        for (int j = 1; j <= n; ++j)
            if (i != j && f[i][j])
                cnt += (g[i] & f[j]).count();
    cout << cnt / 3 << '\n';
} }

\boldsymbol{11.\ 某场\ NOIP\ 模拟赛的\ T1\ 实力}

题意:给定整数 n,要求构造一张无向简单图,满足结点数 n 不超过 500,且恰有 m 个本质不同的三元环。

数据范围:1\le m\le 2\times 10^6

题解:简单构造。考虑下列性质:一张 n 个点的完全图恰有 \binom{n}{3} 个本质不同的三元环。

证明:考虑到任选三个点 i,j,ki<j<k)它们之间都两两有边相连,且任意两组点对都本质不同。显然这样的数对 (i,j,k) 的数量为 \binom{n}{3}。因此得证。

所以说考虑构造若干个完全图然后拼凑出 m。这个可以贪心的让完全图的大小在保证三元环数目之和不超过 m 的前提下尽量的大。直接做打表可以得到三元环数目远远不到 500 个,可以通过。

代码:

int n;
int b3(int x) {
    return x * (x - 1) * (x - 2) / 6;
}
vector<pair<int, int>> edge;
int mp[4145][4145];
void run() {
    freopen("capability.in", "r", stdin);
    freopen("capability.out", "w", stdout);
    n = read();
    int i = 3, idx = 1;
    while (b3(i) <= n) {
        ++i;
    }
    --i;
    while (i >= 3) { 
        int k = n % b3(i);
        int kk = n / b3(i);
        for (int ij = 0; ij < kk; ++ij) {
            int x = idx, y = idx + i - 1;
            for (int i = x; i <= y; ++i)
                for (int j = i + 1; j <= y; ++j)
                    edge.eb(i, j);
            idx += i;
        }
        n = k;
        --i;
    }
    cout << idx - 1 << '\n';
    for (auto &x : edge) {
        int t = x.first, w = x.second;
        mp[t][w - t] = 1;
    }
    for (int i = 1; i < idx - 1; ++i) {
        for (int j = 1; j < idx - i; ++j)
            cout << mp[i][j] << ' ';
        cout << '\n';
    }
}

\boldsymbol{12.\ P4619\ [SDOI2018]\ 旧试题}

题意:求

\sum\limits_{i=1}^A\sum\limits_{j=1}^B\sum\limits_{k=1}^C d(ijk)

10^9+7 取模后的结果。

数据范围:10 组数据,1\le A.B.C\le 10^5\bf{5s}

题解:

神奇反演+三元环计数题,去年noip刚学三元环计数的时候就在题单里了。首先有经典公式 d(ijk)=\sum\limits_{x\mid i}\sum\limits_{y\mid j}\sum\limits_{z\mid k}\epsilon(\gcd(x,y))\epsilon(\gcd(x,z))\epsilon(\gcd(y,z)),因此可以开始推柿子:

\newcommand\sm{\sum\limits_} \begin{aligned} &\sm{i=1}^A\sm{j=1}^B\sm{k=1}^Cd(ijk)\\ =&\sm{i=1}^A\sm{j=1}^B\sm{k=1}^C\sm{x\mid i}\sm{y\mid j}\sm{z\mid k}\epsilon(\gcd(x,y))\epsilon(\gcd(x,z))\epsilon(\gcd(y,z))\\ =&\sm{x=1}^A\sm{i=1}^{\lfloor\frac A{x}\rfloor}\sm{y=1}^B\sm{j=1}^{\lfloor\frac B{y}\rfloor}\sm{z=1}^C\sm{j=1}^{\lfloor\frac C{z}\rfloor}\epsilon(\gcd(x,y))\epsilon(\gcd(x,z))\epsilon(\gcd(y,z))\\ =&\sm{x=1}^A\sm{y=1}^B\sm{z=1}^C\epsilon(\gcd(x,y))\epsilon(\gcd(x,z))\epsilon(\gcd(y,z))\lfloor\frac Ax\rfloor\lfloor\frac By\rfloor\lfloor\frac Cz\rfloor\\ =&\sm{x=1}^A\sm{y=1}^B\sm{z=1}^C\lfloor\frac Ax\rfloor\lfloor\frac By\rfloor\lfloor\frac Cz\rfloor\sm{a\mid x\land a\mid y}\mu(a)\sm{b\mid x\land b\mid y}\mu(b)\sm{c\mid x\land c\mid y}\mu(c)\\ =&\sm{a=1}^{\min(A,B)}\sm{b=1}^{\min(A,C)}\sm{c=1}^{\min(B,C)}\mu(a)\mu(b)\mu(c)\sm{\operatorname{lcm}(a,b)\mid x}\lfloor\frac Ax\rfloor\sm{\operatorname{lcm}(a,c)\mid y}\lfloor\frac By\rfloor\sm{\operatorname{lcm}(b,c)\mid z}\lfloor\frac Cz\rfloor \end{aligned}

后半部分的 \newcommand\sm{\sum\limits_}\sm{\operatorname{lcm}(a,b)\mid x}\lfloor\frac Ax\rfloor\sm{\operatorname{lcm}(a,c)\mid y}\lfloor\frac By\rfloor\sm{\operatorname{lcm}(b,c)\mid z}\lfloor\frac Cz\rfloor 可以分别预处理后 O(1) 查询,因此此时直接做时间复杂度为 O(n^3)。考虑优化。容易发现两点 i,j 之间可能发生转移当且仅当 \mu(i)\neq 0\land \mu(j)\neq 0\land \text{lcm}(i,j)\le\max(A,B,C)。为了不重不漏的计数考虑枚举 i,j 的公约数 g,同时满足 \gcd(i,j)=1,i\times j\times g\le\max(a,b,c) 然后给 i,j 建一条双向边。打个表发现在 a=b=c=10^5 时边数也只有 760741 条。然后考虑用这个图来优化上面的莫反柿子,可以想到上面中每一个对答案有贡献的三元组 (a,b,c) 在图上就是一个三元环,直接搞三元环计数即可。但是由于这里没法处理 a=b,a=c,b=ca=b=c 的情况,所以暴力算出上面两种情况再累加到答案中即可。

总时间复杂度为 O(Tm\sqrt m),其中 m 为边数 \le 760741,卡常后可以通过。

// #pragma GCC optimize(3, "Ofast", "inline", "unroll-loops")
#include <bits/stdc++.h>
#define int long long
using namespace std;

const int N = 100010;
const int mod = 1e9 + 7;

int mu[N], isp[N], pr[N], idx, dfn[N], f[N];
void sieve()
{
    isp[1] = mu[1] = 1;
    for (int i = 2; i < N; ++i)
    {
        if (!isp[i])
        {
            pr[++idx] = i;
            mu[i] = -1;
        }
        for (int j = 1; j <= idx && i * pr[j] < N; ++j)
        {
            isp[i * pr[j]] = 1;
            if (i % pr[j] == 0)
                break;
            mu[i * pr[j]] = -mu[i];
        }
    }
    for (int i = 1; i < N; ++i)
        for (int j = i; j < N; j += i)
            ++f[j];
    for (int i = 1; i < N; ++i)
        f[i] = (f[i] + f[i - 1]) % mod;
}
int fx[N], fy[N], fz[N];
inline int calc(int x, int y)
{
    int cnt = 0;
    for (int i = x; i <= y; i += x)
        cnt += y / i;
    return cnt % mod;
}
inline int calcx(int x, int y, int *fx)
{
    if (x <= y)
        return fx[x];
    return 0;
}
inline int gcd(int a, int b)
{
    int _a = __builtin_ctz(a), _b = __builtin_ctz(b), c = min(_a, _b);
    b >>= _b;
    while (a)
    {
        a >>= _a;
        int diff = a - b;
        _a = __builtin_ctz(diff);
        b = min(a, b), a = abs(diff);
    }
    return b << c;
}
inline int lcm(int a, int b)
{
    return a / gcd(a, b) * b;
}

vector<pair<int, int>> adj[N];
vector<pair<int, int>> adj2[N];
int deg[N];

signed main()
{
    cin.tie(0)->sync_with_stdio(false);
    int T;
    cin >> T;
    sieve();
    // cout << "mu: ";
    // for (int i = 1; i <= 20; ++i)
    //     cout << mu[i] << ' ';
    // cout << '\n';
    while (T--)
    {
        int a, b, c;
        cin >> a >> b >> c;
        int o[3] = {a, b, c};
        sort(o, o + 3);
        tie(a, b, c) = tie(o[0], o[1], o[2]);
        for (int i = 0; i < N; ++i)
        {
            deg[i] = dfn[i] = 0;
            adj[i].clear();
            adj2[i].clear();
            fx[i] = fy[i] = fz[i] = 0;
        }
        for (int i = 1; i <= a; ++i)
            fx[i] = calc(i, a);
        for (int i = 1; i <= b; ++i)
            fy[i] = calc(i, b);
        for (int i = 1; i <= c; ++i)
            fz[i] = calc(i, c);
        // cout << "fx: ";
        // for (int i = 1; i <= a; ++i)
        //     cout << fx[i] << ' ';
        // cout << "\nfy: ";
        // for (int i = 1; i <= b; ++i)
        //     cout << fy[i] << ' ';
        // cout << "\nfz: ";
        // for (int i = 1; i <= c; ++i)
        //     cout << fz[i] << ' ';
        // cout << '\n';
        // for (int g = 1; g <= c; ++g)
        //     if (!isp[g] || g == 1)
        //     {
        //         vector<int> use;
        //         for (int i = g; i <= c; i += g)
        //             if (mu[i])
        //                 use.emplace_back(i);
        //         for (int i = 0; i < use.size(); ++i)
        //             for (int j = i + 1; j < use.size(); ++j)
        //                 if (lcm(use[i], use[j]) <= c)
        //                 {
        //                     adj[use[i]].emplace_back(use[j]);
        //                     adj[use[j]].emplace_back(use[i]);
        //                 }
        //     }
        for (int g = 1; g <= c; ++g)
            for (int i = 1; i * g <= c; ++i)
                if (mu[i * g])
                    for (int j = i + 1; i * j * g <= c; ++j)
                        if (mu[j * g] && gcd(i, j) == 1)
                        {
                            adj[i * g].emplace_back(j * g, i * j * g);
                            adj[j * g].emplace_back(i * g, i * j * g);
                            // cout << "ae " << i * g << ' ' << j * g << '\n';
                            ++deg[i * g], ++deg[j * g];
                        }
        for (int i = 1; i <= c; ++i)
            for (auto &[j, _] : adj[i])
                if (i < j)
                {
                    int u = i, v = j;
                    if (deg[u] > deg[v] || deg[u] == deg[v] && u > v)
                        u ^= v ^= u ^= v;
                    adj2[u].emplace_back(v, _);
                    // cout << "be " << u << ' ' << v << '\n';
                }
        // cout << "deg: ";
        // for (int i = 1; i <= c; ++i) cout << deg[i] << ' '; cout << '\n';
        int s = 0;
        for (int x = 1; x <= c; ++x)
            if (mu[x])
            {
                for (auto &[i, _] : adj2[x])
                    dfn[i] = _;
                for (auto &[y, w1] : adj2[x])
                    if (mu[y])
                        for (auto &[z, w2] : adj2[y])
                            if (dfn[z] && mu[z])
                            {
                                int w3 = dfn[z], w = mu[x] * mu[y] * mu[z];
                                // cout << "ae " << x << ' ' << y << ' ' << z << ' ' << w1 << ' ' << w2 << ' ' << w3 << ' ' << dfn[z] << ' ' << s << ' ' << fx[w1] << ' ' << fy[w2] << ' ' << fz[w3] << "  ";
                                s += w * (fx[w1] * fy[w2] % mod * fz[w3] % mod), s %= mod; //cout << s << ' ';
                                s += w * (fx[w1] * fy[w3] % mod * fz[w2] % mod), s %= mod; //cout << s << ' ';
                                s += w * (fx[w2] * fy[w1] % mod * fz[w3] % mod), s %= mod; //cout << s << ' ';
                                s += w * (fx[w2] * fy[w3] % mod * fz[w1] % mod), s %= mod; //cout << s << ' ';
                                s += w * (fx[w3] * fy[w1] % mod * fz[w2] % mod), s %= mod; //cout << s << ' ';
                                s += w * (fx[w3] * fy[w2] % mod * fz[w1] % mod), s %= mod; //cout << s << ' ';
                                // cout << '\n';
                            }
                for (auto &[i, _] : adj2[x])
                    dfn[i] = 0;
            }
        // cout << s << '\n';
        int last = s;
        for (int g = 1; g <= c; ++g)
            for (int i = 1; i * g <= c; ++i)
                if (mu[i * g])
                    for (int j = i + 1; i * j * g <= c; ++j)
                        if (mu[j * g] && __gcd(i, j) == 1)
                        {
                            s = (s + mu[i * g] * mu[i * g] * mu[j * g] * fx[i * g] * fy[i * j * g] % mod * fz[i * j * g] % mod) % mod;
                            s = (s + mu[i * g] * mu[i * g] * mu[j * g] * fx[i * j * g] * fy[i * g] % mod * fz[i * j * g] % mod) % mod;
                            s = (s + mu[i * g] * mu[i * g] * mu[j * g] * fx[i * j * g] * fy[i * j * g] % mod * fz[i * g] % mod) % mod;
                            s = (s + mu[i * g] * mu[j * g] * mu[j * g] * fx[j * g] * fy[i * j * g] % mod * fz[i * j * g] % mod) % mod;
                            s = (s + mu[i * g] * mu[j * g] * mu[j * g] * fx[i * j * g] * fy[j * g] % mod * fz[i * j * g] % mod) % mod;
                            s = (s + mu[i * g] * mu[j * g] * mu[j * g] * fx[i * j * g] * fy[i * j * g] % mod * fz[j * g] % mod) % mod;
                            // cout << "wa " << i << ' ' << j << ' ' << g << ' ' << s << '\n';
                        }
        // cout << s << '\n';
        // cout << a << ' ' << b << ' ' << c << '\n';
        last = s;
        for (int i = 1; i <= a; ++i)
            s = (s + mu[i] * mu[i] * mu[i] * fx[i] * fy[i] % mod * fz[i] % mod) % mod;
        // cout << s - last << '\n';
        cout << (s % mod + mod) % mod << '\n';
    }
    return 0;
}

有向图三元环计数的另解

考虑对有向图 G,将其中所有的有向边均建立反边,得到无向图 G'

容易证明若 (x,y,z)G 中是合法的三元环,则在 G' 中也为合法的三元环。

因此只需要在新图 G' 中跑三元环计数,然后对于得到的每一个三元环,暴力判断其在有向图 G 中是否也是三元环即可。时间复杂度同样也是 O(m\sqrt m)