First implementation of the Todd-Coxeter algorithm for Coxeter groups. Based on github.com/allemangD/toddcox-fast

solver.cpp

@@ -0,0 +1,231 @@
+#include <vector>
+#include <utility>
+#include <algorithm>
+#include <functional>
+#include <iostream>
+
+struct Cosets {
+    int ngens;
+    std::vector<int> data;
+    int len;
+
+    Cosets(int ngens, std::vector<int> data): ngens(ngens), data(data) {
+        len = data.size() / ngens;
+    }
+
+    void add_row() {
+        len++;
+        data.resize(data.size()+ngens, -1);
+    }
+
+    void put(int coset, int gen, int target) {
+        data[coset * ngens + gen] = target;
+        data[target * ngens + gen] = coset;
+    }
+
+    void put(int idx, int target) {
+        int coset = idx / ngens;
+        int gen = idx % ngens;
+        data[idx] = target;
+        data[target * ngens + gen] = coset;
+    }
+
+    int get(int coset, int gen) {
+        return data[coset * ngens + gen];
+    }
+
+    int get(int idx) {
+        return data[idx];
+    }
+};
+
+struct Mult {
+    int gen0, gen1, mult;
+};
+
+struct RelTable {
+    int gens[2];
+    int mult;
+    std::vector<int> fam, gen, lst;
+    
+    RelTable(Mult m): mult(m.mult) {
+        gens[0] = m.gen0;
+        gens[1] = m.gen1;
+    }
+
+    int add_row() {
+        int idx = fam.size();
+        fam.push_back(-1);
+        gen.push_back(-1);
+        lst.push_back(-1);
+        return idx;
+    }
+};
+
+struct {
+    bool operator()(std::pair<int,int> p1, std::pair<int,int> p2) {
+        return p1.first == p2.first ? p1.second > p2.second : p1.first > p2.first;
+    }
+} intpair_greater;
+
+struct Group {
+    int ngens;
+    std::vector<std::vector<int>> _mults;
+
+    Group(int ngens, std::vector<Mult> rels = {}): ngens(ngens) {
+        _mults.resize(ngens);
+        for (int i = 0; i < ngens; i++) {
+            _mults[i].resize(ngens, 2);
+        }
+
+        for (Mult m : rels) {
+            if (m.gen0 < m.gen1)
+                _mults[m.gen0][m.gen1] = m.mult;
+            else
+                _mults[m.gen1][m.gen0] = m.mult;
+        }
+    }
+
+    void setmult(Mult m) {
+        if (m.gen0 < m.gen1)
+            _mults[m.gen0][m.gen1] = m.mult;
+        else
+            _mults[m.gen1][m.gen0] = m.mult;
+    }
+
+    std::vector<Mult> get_mults() const {
+        std::vector<Mult> mults; 
+        for (int i = 0; i < ngens - 1; i++) {
+            for (int j = i+1; j < ngens; j++) {
+                mults.push_back({i,j,_mults[i][j]});
+            }
+        }
+        return mults;
+    }
+
+    Group operator*(const Group &other) {
+        int off = ngens;
+
+        Group g(ngens + other.ngens, get_mults());
+        
+        for (Mult m : other.get_mults()) {
+            g.setmult({off + m.gen0, off + m.gen1, m.mult});
+        }
+
+        return g;
+    }
+
+    Group operator^(int p) {
+        Group g(ngens * p);
+
+        for (Mult m : get_mults()) {
+            for (int off = 0; off < g.ngens; off += ngens) {
+                g.setmult({off + m.gen0, off + m.gen1, m.mult});
+            }
+        }
+
+        return g;
+    }
+
+    static Group schlafli(std::vector<int> mults) {
+        int ngens = mults.size() + 1;
+        Group g(ngens);
+        for (int i = 0; i < mults.size(); i++) {
+            g.setmult({i, i+1, mults[i]});
+        }
+        return g;
+    }
+
+    Cosets solve(std::vector<int> sub_gens = {}) {
+        std::vector<int> init_row(ngens, -1);
+        for (int i : sub_gens) {
+            init_row[i] = 0;
+        }
+
+        Cosets cosets(ngens, init_row);
+        std::vector<RelTable> rel_tables;
+        for (Mult m : get_mults()) {
+            rel_tables.emplace_back(m);
+        }
+
+        for (RelTable &rel : rel_tables) {
+            int idx = rel.add_row();
+            
+            rel.fam[idx] = 0;
+            rel.gen[idx] = 0;
+            rel.lst[idx] = 0;
+
+            if ( (cosets.get(rel.gens[0]) == 0) xor (cosets.get(rel.gens[1]) == 0) ) {
+                rel.gen[idx] = -1;
+            }
+        }
+
+        int idx = 0;
+        while (true) {
+            while (idx < cosets.data.size() and cosets.get(idx) >= 0)
+                idx++;
+
+            if (idx == cosets.data.size())
+                break;
+
+            int coset = idx / ngens;
+            int gen = idx % ngens;
+            int target = cosets.len;
+
+            cosets.add_row();
+            
+            for (RelTable &rel : rel_tables) {
+                rel.add_row();
+            }
+
+            std::vector<std::pair<int,int>> facts;
+            facts.push_back(std::make_pair(coset, gen));
+
+            while (!facts.empty()) {
+                auto top = facts.back();
+                coset = top.first;
+                gen = top.second;
+                facts.pop_back();
+
+                cosets.put(coset, gen, target);
+
+                for (RelTable &rel : rel_tables) {
+                    if ( (gen == rel.gens[0] or gen == rel.gens[1]) and (rel.fam[target] == -1) ) {
+                        rel.fam[target] = rel.fam[coset];
+                        rel.gen[target] = rel.gen[coset] + 1;
+                        
+                        if (rel.gen[coset] < 0)
+                            rel.gen[target] -= 2;
+
+                        if (rel.gen[target] == rel.mult) {
+                            int lst = rel.lst[rel.fam[target]];
+                            int gen_ = rel.gens[(int)(rel.gens[0] == gen)];
+                            facts.push_back(std::make_pair(lst, gen_));
+                        }
+                        else if (rel.gen[target] == -rel.mult) {
+                            int gen_ = rel.gens[rel.gens[0] == gen];
+                            facts.push_back(std::make_pair(target, gen_));
+                        }
+                        else if (rel.gen[target] == rel.mult - 1) {
+                            rel.lst[rel.fam[target]] = target;
+                        }
+                    }
+                }
+
+                std::sort(facts.begin(), facts.end(), intpair_greater);
+            }
+
+            for (RelTable &rel : rel_tables) {
+                if (rel.fam[target] == -1) {
+                    rel.fam[target] = target;
+                    rel.gen[target] = 0;
+
+                    if ( (cosets.get(target, rel.gens[0]) == target) xor (cosets.get(target, rel.gens[1]) == target) )
+                        rel.gen[target] = -1;
+                }
+            }
+        }
+        
+        return cosets;
+    }
+};