Betti Numbers

1
#include <bits/stdc++.h>
2
using namespace std;
3
4
struct DSU {
5
    int n; vector<int> p, r;
6
    DSU(int n=0): n(n), p(n), r(n,0) { iota(p.begin(), p.end(), 0); }
7
    int find(int x){ return p[x]==x? x: p[x]=find(p[x]); }
8
    bool unite(int a, int b){ a=find(a); b=find(b); if(a==b) return false; if(r[a]<r[b]) swap(a,b); p[b]=a; if(r[a]==r[b]) r[a]++; return true; }
9
};
10
11
// Example graph: two components, one contains a triangle (one independent cycle)
12
// V = {0,1,2,3,4}, E = {(0,1),(1,2),(2,0),(3,4)}
13
int main(){
14
    int V = 5;
15
    vector<pair<int,int>> edges = {{0,1},{1,2},{2,0},{3,4}};
16
17
    DSU dsu(V);
18
    for (auto [u,v] : edges) dsu.unite(u,v);
19
20
    // Count components (β0)
21
    unordered_set<int> reps;
22
    for (int v=0; v<V; ++v) reps.insert(dsu.find(v));
23
    int beta0 = (int)reps.size();
24
25
    // β1 for a graph: |E| - |V| + β0
26
    int beta1 = (int)edges.size() - V + beta0;
27
28
    cout << "Graph Betti numbers:\n";
29
    cout << "beta0 (components) = " << beta0 << "\n";
30
    cout << "beta1 (independent cycles) = " << beta1 << "\n";
31
32
    return 0;
33
}
34

1
#include <bits/stdc++.h>
2
using namespace std;
3
4
// Bit-packed matrix over GF(2), stored as rows of uint64_t words
5
struct BitMatrix {
6
    int rows, cols, W; // W = words per row
7
    vector<vector<uint64_t>> a; // a[rows][W]
8
    BitMatrix(int r=0, int c=0): rows(r), cols(c) {
9
        W = (cols + 63) / 64;
10
        a.assign(rows, vector<uint64_t>(W, 0ULL));
11
    }
12
    inline void set1(int r, int c){
13
        int w = c >> 6, b = c & 63;
14
        a[r][w] |= (1ULL << b);
15
    }
16
    inline bool get(int r, int c) const {
17
        int w = c >> 6, b = c & 63;
18
        return (a[r][w] >> b) & 1ULL;
19
    }
20
    inline void xor_row(int i, const vector<uint64_t>& pivot){
21
        for(int w=0; w<W; ++w) a[i][w] ^= pivot[w];
22
    }
23
};
24
25
int rank_mod2(BitMatrix M){
26
    int r = 0; // current pivot row
27
    for(int col=0; col<M.cols && r < M.rows; ++col){
28
        int sel = -1;
29
        for(int i=r; i<M.rows; ++i){ if (M.get(i,col)) { sel = i; break; } }
30
        if (sel == -1) continue; // no pivot in this column
31
        swap(M.a[sel], M.a[r]);
32
        // Eliminate this column in all other rows
33
        for(int i=0; i<M.rows; ++i){ if (i!=r && M.get(i,col)) M.xor_row(i, M.a[r]); }
34
        ++r;
35
    }
36
    return r; // number of pivots = rank
37
}
38
39
struct PairHash { size_t operator()(const pair<int,int>& p) const noexcept { return (uint64_t)p.first<<32 ^ (uint64_t)p.second; } };
40
41
// Build edges and boundary matrices for a 2D simplicial complex
42
// Input: V vertices labeled 0..V-1, list of triangles (each as 3 vertex indices)
43
// Output: beta0, beta1, beta2 over Z2
44
struct Betti2D {
45
    int V, E, T;
46
    int beta0, beta1, beta2;
47
};
48
49
Betti2D compute_betti_2D(int V, const vector<array<int,3>>& triangles){
50
    // Deduplicate edges from triangles
51
    unordered_map<pair<int,int>, int, PairHash> edge_id;
52
    vector<pair<int,int>> edges;
53
    auto add_edge = [&](int u, int v){ if(u>v) swap(u,v); auto key = make_pair(u,v); if(!edge_id.count(key)){ int id = (int)edges.size(); edge_id[key]=id; edges.push_back(key); } };
54
55
    for (auto tri : triangles){
56
        array<int,3> t = tri;
57
        sort(t.begin(), t.end()); // canonicalize triangle vertex order (not strictly needed over Z2)
58
        add_edge(t[0], t[1]); add_edge(t[0], t[2]); add_edge(t[1], t[2]);
59
    }
60
    int E = (int)edges.size();
61
    int T = (int)triangles.size();
62
63
    // Build boundary matrices over Z2: \partial_1: C1->C0 is V x E; \partial_2: C2->C1 is E x T
64
    BitMatrix d1(V, E);
65
    for (int e=0; e<E; ++e){
66
        auto [u,v] = edges[e];
67
        d1.set1(u, e);
68
        d1.set1(v, e);
69
    }
70
71
    BitMatrix d2(E, T);
72
    for (int t=0; t<T; ++t){
73
        int a = triangles[t][0], b = triangles[t][1], c = triangles[t][2];
74
        // canonicalize edges (u<v) to match edge_id keys
75
        auto key1 = make_pair(min(a,b), max(a,b));
76
        auto key2 = make_pair(min(a,c), max(a,c));
77
        auto key3 = make_pair(min(b,c), max(b,c));
78
        int e1 = edge_id[key1];
79
        int e2 = edge_id[key2];
80
        int e3 = edge_id[key3];
81
        d2.set1(e1, t); d2.set1(e2, t); d2.set1(e3, t);
82
    }
83
84
    int rank_d1 = rank_mod2(d1);
85
    int rank_d2 = rank_mod2(d2);
86
87
    // Betti numbers over Z2 for a 2D complex (no 3-simplices):
88
    // beta0 = m0 - rank(d1), beta1 = m1 - rank(d1) - rank(d2), beta2 = m2 - rank(d2)
89
    Betti2D out; out.V=V; out.E=E; out.T=T;
90
    out.beta0 = V - rank_d1;
91
    out.beta1 = E - rank_d1 - rank_d2;
92
    out.beta2 = T - rank_d2; // since \partial_3 = 0 (no tetrahedra)
93
    return out;
94
}
95
96
int main(){
97
    // Example: boundary of a tetrahedron (a sphere S^2)
98
    // Vertices: 0,1,2,3; Triangles: all 4 faces
99
    int V = 4;
100
    vector<array<int,3>> tris = {
101
        {0,1,2}, {0,1,3}, {0,2,3}, {1,2,3}
102
    };
103
    auto betti = compute_betti_2D(V, tris);
104
105
    cout << "2D complex stats: V=" << betti.V << ", E=" << betti.E << ", T=" << betti.T << "\n";
106
    cout << "Betti numbers over Z2:" << "\n";
107
    cout << "beta0 (components) = " << betti.beta0 << "\n";
108
    cout << "beta1 (tunnels)     = " << betti.beta1 << "\n";
109
    cout << "beta2 (voids)       = " << betti.beta2 << "\n";
110
111
    // Expected for sphere: beta0=1, beta1=0, beta2=1
112
    return 0;
113
}
114