【学习笔记】图上的环计数问题
Priestess_SLG · · 算法·理论
毕竟这也算是一个冷门的知识点了吧,所以没找到几道题 qwq
无向图三元环计数
题目:给定一张
前置定义:
一般法
对于
- 若
\text{deg}_u\le\text{deg}_v ,则令新图中边为u\to v 。 - 若
\text{deg}_u>\text{deg}_v ,则令新图中边为v\to u 。
记新图
此时有一个十分优美的结论:
-
结论
\bf1 :对于图中每一个点i(1\le i\le n) ,\text{deg}_i 都不超过\sqrt{2m} 。 -
证明: 考虑到对边定向并不会让每一个点的度数增加,因此
- 若
\text{deg}_i\le \sqrt {2m} ,则有\text{deg}'_i\le\text{deg}_i\le \sqrt {2m} 。 - 若
\text{deg}_i>\sqrt {2m} ,则根据上述的连边方式,i\to j 在新图中有边,当且仅当\text{deg}_j\ge\text{deg}_i>\sqrt {2m} 。因为\sum\text{deg}_i=2m ,所以满足上述条件的j 必不超过\sqrt{2m} 个。于是\text{deg}'_i\le\sqrt{2m} 。
- 若
然后考虑开始计数。考虑对于原无向图
根据上面的结论一,对于每一个点,其出度不会超过
做到严格
上述的标记结点部分也可以通过设置时间戳来做。
坑点:
-
m>10^5 - 排序的时候若度数相同,则需要根据点的编号大小排序
代码:
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 法
考虑对于每一个点
稠密图的情况下很好用,且比普通法好写。
代码(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
这里曾经有一张图片,但是由于某些原因丢失了。
无向图四元环计数
题目:给定一张
特殊的,认为
还是沿用上面的做法,把所有点按照度数从小到大排序,然后按照上面的做法建立新图,并建立其的反图。 考虑枚举到图上的一个点
考虑枚举
和上面的分析相似,时间复杂度还是
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}
(不是为什么我刚学图上计数就开这么折磨的题啊)
题意:给定一张
数据范围:
题解:
考虑分类讨论。发现上述选取的方案若合法则必然生成
- 四元环
- 三元环外加一个点
- 菊花(十字)
- 五个点的链
- 四个点的链外加一个点
容易发现五种情况的并集恰好为全集,且两两的交为空。因此只需要分别计算五种情况的答案并求和即可。
(1) 四元环
直接套四元环计数的板子即可。
(2) 三元环外加一个点
考虑继续沿用计算三元环的做法。因为三元环的每一个结点均与另两个结点有连边,且做三元环计数的时候枚举了每一个合法的三元环,因此考虑对于每一个三元环
(3) 菊花(十字)
考虑枚举菊花的重心结点
(4) 四个点的链外加一个点
考虑先枚举该形态的两个顶点
但是这显然是错的,因为
(5) 五个点的链
考虑枚举链的重心结点
令
于是讨论完了这题的所有情况。
加强版: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}
题意:给定一张
数据范围:
题解:
基础题。
考虑到计数三元环的时候是枚举了每一个三元环的,因此考虑先枚举每一个三元环,然后对每一条边记录有多少个三元环使用了该边。设有
代码:
不知道为什么 MLE 了,我寻思我只用了 7MB 的内存啊/yiw
\boldsymbol{3.\ QOJ 6534,\ Um\_nik\ mod\ 998\ 244\ 353\ Contest\ K\ 4}
题意:给定一张
数据范围:
题解:考虑枚举
然后把图上所有的边都定向,使用 bitset 对
时间复杂度为什么是对的?证明:考虑令
为了让
总时间复杂度即为
代码:
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}
题意:给定一张
数据范围:
题解:
容易发现为了让路径长度最短,有且仅可能有下面的三种情况:
- 全走长度为
a 的边(1 ) - 在原最短路上将相邻两条边替换为长度为
b 的边,剩下的走长度为a 的边(2 ) - 全走长度为
b 的边(3 )
首先
考虑枚举任意两条边
考虑这样的一个结论:
结论:对于三个点
x,y,z ,定义\text{dis}_i 表示从起点到i 点的最短路的长度,则若x\to y 和y\to z 两条边边权的和大于或等于x\to z 这条边边权,且x,y,z 三个点不两两之间有边,则可以直接删去y\to z 这条边。
感性理解容易发现结论正确。然后考虑到
证明:考虑三元环将边定向后得到的新有向图
于是直接在松弛的时候删边,对三元环的情况特判。还剩下
代码:
vector 存图想要删除某一个元素,若具体顺序不重要,则可以先把这个元素和末尾元素交换,然后再删除末尾元素,做到
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}
题意:给定
数据范围:
题解:
貌似是另一种套路?
考虑建图。先对序列离散化,然后将每一个序列中的元素向对应的序列连边。此时若找到任意一组四元环
然后问题就解决了。时间复杂度为
代码:
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}
题意:
有一张
其中 如:[T_TLucas_Yin 是可爱班花] 因括号内表达式为真所以答案为 (不是)
数据范围:
题解:
首先如果要对两两之间有边的情况计数则相对简单,直接计数三元环然后对于每一个三元环计算其贡献。但是问题要对两两之间无边计数。建立反图可以解决这个问题但是反图边数是
于是套路的考虑容斥答案。这样需要计数:
- 图中所有三元组对答案的贡献。容斥系数为
1 。 - 图中所有三元环对答案的贡献。容斥系数为
-1 。 - 图中所有三元组满足其中有两对点对有边另一对点对无边对答案的贡献。容斥系数为
1 。 - 图中所有三元组满足其中有一对点对有边另一对点对无边对答案的贡献。容斥系数为
-1 。
将上述部分的答案和容斥系数相乘后求和变为答案。问题是计算上述四个问题的答案。
(1) 所有三元组对答案的贡献
此时和图的形态无关,式子为:
考虑到最后的柿子
1.
考虑抽离出和
发现
可以直接计数。
2.
考虑抽离出和
发现该表达式和
此时式子只和
也可以直接计数。
3.
和
还是可以直接计数。
(2) 三元环对答案的贡献
考虑到做三元环计数的时候其实是枚举到了每一个三元环的,因此直接对于每一个三元环计算其对答案的贡献即可。
(3) 所有三元组满足其中有一对点对有边另两对点对无边对答案的贡献
考虑枚举该三元组中的一条边
只需要将上述七种情况的答案求和即可。
(4) 所有三元组满足其中有两对点对有边另一对点对无边对答案的贡献
最困难的一集
考虑枚举在这个大小为
一、
钦定三元组为
将上述三种情况的贡献求和即可得到
二、
考虑对
题意:给定一个
数据范围:
题解:
仍然考虑求三元环时的定向方法,让度数小的点连向度数大的点,形成一张新的有向图。此时每一个点的出度不超过
然后考虑枚举点
时间复杂度证明:显然为了让时间复杂度最大,让新图中点的出度的最大值为
代码:
\boldsymbol{8.\ CF11D\ A\ Simple\ Task}
题意:给定一张
数据范围:
题解:
这个好像没有什么特别优秀的解法……考虑到
于是设
写代码的话需要注意一下一个环会被正着计数一遍然后再被反着计数一遍,答案需要除以
代码:
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}
题意:给定一张
数据范围:
题解:
板子题。考虑到做三元环计数的时候具体枚举到了每一个三元环,所以直接对每一个三元环处理其贡献即可。时间复杂度为
代码:
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}
题意:给定一张有向稠密图
数据范围:
题解:有向图无法使用给边定向,但是发现
记录第一个 bitset:
记录第二个 bitset:
然后考虑枚举三元环中的两个结点
注意点:
bitset有向图会把三元环重复计算3 遍,而无向图中是6 遍。- 需要开文件读写! 输入文件名为 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\ 实力}
题意:给定整数
数据范围:
题解:简单构造。考虑下列性质:一张
证明:考虑到任选三个点
所以说考虑构造若干个完全图然后拼凑出
代码:
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]\ 旧试题}
题意:求
对
数据范围:
题解:
神奇反演+三元环计数题,去年noip刚学三元环计数的时候就在题单里了。首先有经典公式
后半部分的
总时间复杂度为
// #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;
}
有向图三元环计数的另解
考虑对有向图
容易证明若
因此只需要在新图