add context group to tc::Path and tc::Cosets

.gitignore

@@ -2,5 +2,3 @@
 
 *.[oa]
 cmake-build*
-
-Testing/Temporary

.gitmodules

@@ -1,3 +0,0 @@
-[submodule "ext/eigen"]
-	path = ext/eigen
-	url = https://gitlab.com/libeigen/eigen.git

CMakeLists.txt

@@ -1,22 +1,15 @@
 cmake_minimum_required(VERSION 3.10)
 project(toddcox-faster)
 
-set(CMAKE_CXX_STANDARD 17)
-
 option(TC_BUILD_EXAMPLE "Build example executables" OFF)
 
-add_subdirectory(ext)
+add_library(tc STATIC
+        src/groups.cpp
+        src/solve.cpp
+        src/core.cpp)
 
-add_library(tc INTERFACE)
+target_include_directories(tc PUBLIC include)
-target_link_libraries(tc INTERFACE eigen)
 
-target_include_directories(tc INTERFACE include)
+if (TC_BUILD_EXAMPLE)
-
+    add_subdirectory(example)
-if (TC_BUILD_EXAMPLES)
-    add_subdirectory(examples)
 endif ()
-
-if (TC_BUILD_TESTING)
-    enable_testing()
-    add_subdirectory(Testing)
-endif()

Testing/CMakeLists.txt

@@ -1,22 +0,0 @@
-cmake_policy(SET CMP0110 NEW)
-
-add_executable(schlafli schlafli.cpp)
-target_link_libraries(schlafli PRIVATE tc)
-
-add_executable(special special.cpp)
-target_link_libraries(special PRIVATE tc)
-
-add_test(NAME "schalfi [6]" COMMAND schlafli "6" "" "12")
-add_test(NAME "schalfi [6] [0]" COMMAND schlafli "6" "0" "6")
-add_test(NAME "schalfi [6] [1]" COMMAND schlafli "6" "1" "6")
-
-add_test(NAME "schalfi [5 3 3] []" COMMAND schlafli "5 3 3" "" "14400")
-add_test(NAME "schalfi [5 3 3] [0 1]" COMMAND schlafli "5 3 3" "0 1" "1440")
-add_test(NAME "schalfi [5 3 3] [1 2 3]" COMMAND schlafli "5 3 3" "1 2 3" "600")
-
-add_test(NAME "special E6" COMMAND special "E6" "" "51840")
-add_test(NAME "special E7" COMMAND special "E7" "" "2903040")
-#add_test(NAME "special E8" COMMAND special "E8" "" "696729600")  # too big.
-add_test(NAME "special B6" COMMAND special "B6" "" "46080")
-add_test(NAME "special B7" COMMAND special "B7" "" "645120")
-add_test(NAME "special B8" COMMAND special "B8" "" "10321920")

Testing/common.hpp

@@ -1,30 +0,0 @@
-#pragma once
-
-#include <string>
-#include <vector>
-#include <sstream>
-
-#include <tc/group.hpp>
-#include <tc/groups.hpp>
-#include <tc/solver.hpp>
-
-tc::Symbol parse_vec(const std::string &part) {
-    std::istringstream iss(part);
-
-    std::vector<unsigned int> vec;
-
-    std::string token;
-    while (std::getline(iss, token, ' ')) {
-        vec.push_back(std::stoul(token));
-    }
-
-    return Eigen::Map<tc::Symbol>(vec.data(), vec.size());
-}
-
-size_t compute(
-    const tc::Group &group,
-    const tc::Symbol &gens
-) {
-    auto table = tc::solve(group, gens);
-    return table.order();
-}

Testing/schlafli.cpp

@@ -1,17 +0,0 @@
-#include "common.hpp"
-
-#include <iostream>
-
-int main(int argc, char **argv) {
-    auto vsymbol = parse_vec(argv[1]);
-    auto vgens = parse_vec(argv[2]);
-    auto target = std::stoul(argv[3]);
-
-    tc::Symbol symbol(vsymbol.size());
-    symbol << Eigen::Map<tc::Symbol>(vsymbol.data(), vsymbol.size());
-    tc::Group group = tc::schlafli(symbol);
-    auto order = compute(group, vgens);
-
-    std::cout << "Order: " << order << ":" << target << std::endl;
-    return order != target;
-}

Testing/special.cpp

@@ -1,37 +0,0 @@
-#include "common.hpp"
-
-#include <iostream>
-
-int main(int argc, char **argv) {
-    auto name = std::string(argv[1]);
-    auto vgens = parse_vec(argv[2]);
-    auto target = std::stoul(argv[3]);
-
-    unsigned int order;
-
-    tc::Group group(0);
-
-    if (name == "E6") {
-        group = tc::group::E(6);
-    }
-    if (name == "E7") {
-        group = tc::group::E(7);
-    }
-    if (name == "E8") {
-        group = tc::group::E(8);
-    }
-    if (name == "B6") {
-        group = tc::group::B(6);
-    }
-    if (name == "B7") {
-        group = tc::group::B(7);
-    }
-    if (name == "B8") {
-        group = tc::group::B(8);
-    }
-
-    order = compute(group, vgens);
-
-    std::cout << "Order: " << order << ":" << target << std::endl;
-    return order != target;
-}

example/CMakeLists.txt

@@ -0,0 +1,5 @@
+add_executable(bench bench.cpp)
+target_link_libraries(bench PRIVATE tc)
+
+add_executable(path path.cpp)
+target_link_libraries(path PRIVATE tc)

example/bench.cpp

@@ -0,0 +1,35 @@
+#include "tc/core.hpp"
+#include "tc/groups.hpp"
+
+#include <ctime>
+#include <vector>
+#include <iostream>
+
+int main() {
+    std::vector<tc::Group> groups = {
+        tc::group::H(2),
+        tc::group::H(3),
+        tc::group::H(4),
+        tc::group::T(100),
+        tc::group::T(500),
+        tc::group::T(1000),
+        tc::group::E(6),
+        tc::group::E(7),
+        tc::group::B(6),
+        tc::group::B(7),
+        tc::group::B(8)
+    };
+
+    for (const auto &group : groups) {
+        auto s = std::clock(); // to measure CPU time
+        auto cosets = group.solve();
+        auto e = std::clock();
+
+        double diff = (double) (e - s) / CLOCKS_PER_SEC;
+        int order = cosets.size();
+
+        std::cout << group.name << "," << order << "," << diff << std::endl;
+    }
+
+    return 0;
+}

example/path.cpp

@@ -0,0 +1,21 @@
+#include "tc/groups.hpp"
+
+#include <ctime>
+#include <iostream>
+
+int main() {
+    auto cube = tc::group::B(3);
+    auto vars = cube.solve();
+
+    auto words = vars.path.walk<std::string, std::string>(
+        "",
+        {"a", "b", "c"},
+        [](auto a, auto b) { return a + b; }
+    );
+
+    for (const auto &word : words) {
+        std::cout << word << std::endl;
+    }
+
+    return 0;
+}

examples/CMakeLists.txt

@@ -1,11 +0,0 @@
-add_executable(bench bench.cpp)
-target_link_libraries(bench PRIVATE tc)
-
-add_executable(path path.cpp)
-target_link_libraries(path PRIVATE tc)
-
-add_executable(group group.cpp)
-target_link_libraries(group PRIVATE tc)
-
-add_executable(complex complex.cpp)
-target_link_libraries(complex PRIVATE tc)

examples/bench.cpp

@@ -1,41 +0,0 @@
-#include <ctime>
-#include <iostream>
-#include <iomanip>
-
-#include <tc/solver.hpp>
-#include <tc/groups.hpp>
-
-template<class G>
-void test(const G &group) {
-    tc::Symbol gens(1);
-    gens << 0;
-
-    auto s = std::clock();
-    auto cosets = tc::solve(group, gens);
-    auto e = std::clock();
-
-    double diff = (double) (e - s) / CLOCKS_PER_SEC;
-    int order = cosets.order();
-
-    std::cout
-        << std::setw(7) << group.name << ", "
-        << std::setw(7) << order << ", "
-        << std::fixed << std::setprecision(6) << diff << "s"
-        << std::endl;
-}
-
-int main() {
-    test(tc::group::H(2));
-    test(tc::group::H(3));
-    test(tc::group::H(4));
-    test(tc::group::T(100));
-    test(tc::group::T(500));
-    test(tc::group::T(1000));
-    test(tc::group::E(6));
-    test(tc::group::E(7));
-    test(tc::group::B(6));
-    test(tc::group::B(7));
-    test(tc::group::B(8));
-
-    return 0;
-}

examples/complex.cpp

@@ -1,39 +0,0 @@
-#include <ctime>
-#include <iostream>
-#include <iomanip>
-
-#include <tc/complex.hpp>
-#include <tc/groups.hpp>
-
-template<class G>
-void test(const G &group, unsigned int N) {
-    auto s = std::clock();
-    auto combos = tc::combinations(group.gens, N - 1);
-    auto data = tc::merge(tc::hull(group, combos));
-    auto e = std::clock();
-
-    double diff = (double) (e - s) / CLOCKS_PER_SEC;
-    int count = data.cols();
-
-    std::cout
-        << std::setw(2) << N << ", "
-        << std::setw(7) << group.name << ", "
-        << std::setw(7) << count << ", "
-        << std::fixed << std::setprecision(6) << diff << "s"
-        << std::endl;
-}
-
-int main() {
-    test(tc::group::H(4), 4);
-    test(tc::group::B(4), 4);
-    test(tc::group::B(5), 4);
-    test(tc::group::B(6), 4);
-    test(tc::group::E(6), 4);
-    test(tc::group::H(3), 3);
-    test(tc::group::B(4), 3);
-    test(tc::group::B(5), 3);
-    test(tc::group::B(6), 3);
-    test(tc::group::E(6), 3);
-
-    return 0;
-}

examples/group.cpp

@@ -1,35 +0,0 @@
-#include <iostream>
-
-#include <tc/solver.hpp>
-#include <tc/groups.hpp>
-
-int main() {
-    tc::Symbol symbol(4);
-    symbol << 5, 3, 2, 3;
-    tc::Group group = tc::schlafli(symbol);
-//    tc::Group group = tc::group::E(7);
-
-    size_t srank = 3;
-
-    tc::SubGroups subs = tc::subgroups(group, srank);
-
-    std::cout << "Group " << group.name << " (" << subs.size() << " subgroups)" << std::endl;
-    std::cout << group << std::endl;
-
-    for (const auto &sub: subs) {
-        for (int i = 0; i < srank; ++i) {
-            for (int j = 0; j < srank; ++j) {
-                auto sub_mult = sub(i, j);
-                auto src_mult = group(sub.gens(i), sub.gens(j));
-
-                if (sub_mult != src_mult) {
-                    std::cout << "Incorrect subgroup " << sub.name << std::endl;
-                    std::cout << sub << std::endl;
-                    return 1;
-                }
-            }
-        }
-    }
-
-    return 0;
-}

examples/path.cpp

@@ -1,22 +0,0 @@
-#include <vector>
-#include <string>
-#include <iostream>
-
-#include <tc/solver.hpp>
-#include <tc/groups.hpp>
-
-int main() {
-    tc::Symbol gens(0);
-    auto cube = tc::group::B(3);
-    auto vars = tc::solve(cube, gens);
-
-    std::string start;
-    std::vector<std::string> names = {"a", "b", "c"};
-    auto words = vars.path().walk(start, names, std::plus<>());
-
-    for (const auto &word: words) {
-        std::cout << (word.empty() ? "-" : word) << std::endl;
-    }
-
-    return 0;
-}

ext/CMakeLists.txt

@@ -1,2 +0,0 @@
-add_library(eigen INTERFACE)
-target_include_directories(eigen INTERFACE eigen)

include/tc/complex.hpp

@@ -1,170 +0,0 @@
-#pragma once
-
-#include <tc/group.hpp>
-#include <tc/solver.hpp>
-#include <cmath>
-#include <optional>
-#include <numeric>
-#include <iostream>
-#include <utility>
-
-namespace tc {
-    std::vector<Symbol> combinations(const Symbol &symbol, size_t srank) {
-        size_t rank = symbol.size();
-
-        std::vector<bool> mask(rank, false);
-        std::fill_n(mask.begin(), srank, true);
-
-        std::vector<Symbol> combos;
-        combos.reserve(choose(rank, srank));
-
-        Symbol row(srank);
-        do {
-            for (int j = 0, k = 0; j < rank; ++j) {
-                if (mask[j]) {
-                    row(k++) = symbol(j);
-                }
-            }
-            combos.emplace_back(row);
-        } while (std::prev_permutation(mask.begin(), mask.end()));
-
-        return combos;
-    }
-
-    // todo remove
-    ArrayXui fan(const ArrayXui &prims, int root) {
-        ArrayXui res(prims.rows() + 1, prims.cols());
-        res.topRows(prims.rows()) << prims;
-        res.bottomRows(1).fill(root);
-        return res;
-    }
-
-    // todo visitor
-    void apply(const tc::Cosets &table, unsigned int gen, ArrayXui &prims) {
-        for (Eigen::Index i = 0; i < prims.size(); ++i) {
-            prims(i) = table.get(prims(i), gen);
-        }
-    }
-
-/**
- * Convert the indexes of this mesh to those of a different context, using g_gens to build the parent context and sg_gens to build this context.
- */
-    void recontext(
-        ArrayXui &prims,
-        const tc::Group &context,
-        const Symbol &g_gens,
-        const Symbol &sg_gens
-    ) {
-        const auto proper_sg_gens = recontext_gens(context.rank(), g_gens, sg_gens);
-
-        const auto table = solve(context, g_gens, Symbol(0));
-        const auto path = solve(context, sg_gens, Symbol(0)).path();
-
-        auto map = path.walk(0U, proper_sg_gens, [&table](auto coset, auto gen) {
-            return table.get(coset, gen);
-        });
-
-        // todo visitor
-        for (Eigen::Index i = 0; i < prims.size(); ++i) {
-            prims(i) = map[prims(i)];
-        }
-    }
-
-    // todo remove
-    ArrayXui merge(const std::vector<ArrayXui> &meshes) {
-        Eigen::Index cols = 0;
-        for (const auto &mesh: meshes) {
-            cols += mesh.cols();
-        }
-
-        ArrayXui res(meshes[0].rows(), cols);
-        Eigen::Index offset = 0;
-        for (const auto &mesh: meshes) {
-            res.middleCols(offset, mesh.cols()) << mesh;
-            offset += mesh.cols();
-        }
-
-        return res;
-    }
-
-    std::vector<ArrayXui> each_tile(
-        ArrayXui base,
-        const tc::Group &context,
-        const Symbol &g_gens,
-        const Symbol &sg_gens
-    ) {
-        recontext(base, context, g_gens, sg_gens);
-
-        const auto table = solve(context, g_gens, Symbol(0));
-        const auto path = solve(context, g_gens, sg_gens).path();
-
-        auto _gens = context.gens;
-
-        auto res = path.walk(base, _gens, [&table](auto from, auto &gen) {
-            apply(table, gen, from);
-            return from;
-        });
-
-        return res;
-    }
-
-    /**
- * Produce a mesh of primitives that fill out the volume of the subgroup generated by generators g_gens within the group context
- */
-    ArrayXui triangulate(
-        const tc::Group &context,
-        const Symbol &g_gens
-    ) {
-        if (g_gens.size() == 0) {
-            return ArrayXui::Zero(1, 1);
-        }
-
-        const auto &combos = combinations(g_gens, g_gens.size() - 1);
-
-        std::vector<ArrayXui> meshes;
-
-        // todo inline logic
-        //  erase, merge, and fan can be inlined
-        //  for 1..#tiles
-        //   cols += cols
-        //   parts.append(cols)
-        //  result(rows, cols)
-        //  for 1..#parts
-        //   result.middlecols << part
-        //  result.bottomrow.fill(0)
-
-        // todo subgroup/coset metadata
-        //  would be good to also output which coset of which subgroup each primitive is a part of.
-        //  this could be used in shaders etc to make rendering things easier
-
-        for (const auto &sg_gens: combos) {
-            auto base = triangulate(context, sg_gens);
-            auto tiles = each_tile(base, context, g_gens, sg_gens);
-            tiles.erase(tiles.begin(), tiles.begin() + 1);
-            auto raised = merge(tiles);
-            auto fanned = fan(raised, 0);
-            meshes.push_back(fanned);
-        }
-
-        const ArrayXui &result = merge(meshes);
-        return result;
-    }
-
-    template<class T>
-    auto hull(const tc::Group &group, T all_sg_gens) {
-        std::vector<ArrayXui> parts;
-        auto g_gens = group.gens;
-
-        // todo inline logic
-        //  should be able to inline in a similar way as is possible in triangulate
-
-        for (const Symbol &sg_gens: all_sg_gens) {
-            const auto &base = triangulate(group, sg_gens);
-            const auto &tiles = each_tile(base, group, g_gens, sg_gens);
-            for (const auto &tile: tiles) {
-                parts.push_back(tile);
-            }
-        }
-        return parts;
-    }
-}

include/tc/core.hpp

@@ -0,0 +1,151 @@
+#pragma once
+
+#include <array>
+#include <functional>
+#include <vector>
+#include <string>
+
+namespace tc {
+    struct Action;
+    struct Path;
+    struct Cosets;
+    struct Rel;
+    struct Group;
+    struct SubGroup;
+
+    struct Action {
+        int from_idx = -1;
+        int gen = -1;
+
+        Action() = default;
+
+        Action(const Action &) = default;
+
+        Action(int from_idx, int gen);
+    };
+
+    struct Path {
+        const tc::Group &context;
+        std::vector<Action> path;
+
+        explicit Path(const tc::Group &context);
+
+        Path(const Path &) = default;
+
+        void add_row();
+
+        [[nodiscard]] Action get(int to_idx) const;
+
+        void put(int from_idx, int gen, int to_idx);
+
+        template<class T, class E>
+        [[nodiscard]] std::vector<T> walk(
+            T start,
+            std::vector<E> gens,
+            std::function<T(const T &, const E &)> op
+        ) const {
+            std::vector<T> res(size());
+            res[0] = start;
+
+            for (int i = 1; i < res.size(); ++i) {
+                auto &action = path[i];
+                auto &from = res[action.from_idx];
+                auto &val = gens[action.gen];
+                res[i] = op(from, val);
+            }
+
+            return res;
+        }
+
+        template<class T>
+        [[nodiscard]] std::vector<T> walk(
+            T start,
+            std::function<T(const T &, const int &)> op
+        ) const {
+            std::vector<T> res(size());
+            res[0] = start;
+
+            for (int i = 1; i < res.size(); ++i) {
+                auto &action = path[i];
+                auto &from = res[action.from_idx];
+                auto &val = action.gen;
+                res[i] = op(from, val);
+            }
+
+            return res;
+        }
+
+        [[nodiscard]] size_t size() const;
+    };
+
+    struct Cosets {
+        const tc::Group &context;
+
+        std::vector<int> data;
+        Path path;
+
+        Cosets(const Cosets &) = default;
+
+        explicit Cosets(const tc::Group &context);
+
+        void add_row();
+
+        void put(int coset, int gen, int target);
+
+        void put(int idx, int target);
+
+        [[nodiscard]] int get(int coset, int gen) const;
+
+        [[nodiscard]] int get(int idx) const;
+
+        [[nodiscard]] size_t size() const;
+    };
+
+    struct Rel {
+        std::array<int, 2> gens;
+        int mult;
+
+        Rel() = default;
+
+        Rel(const Rel &) = default;
+
+        Rel(int a, int b, int m);
+
+        [[nodiscard]] Rel shift(int off) const;
+    };
+
+    struct Group {
+        const int ngens;
+        std::vector<std::vector<int>> _mults;
+        std::string name;
+
+        Group(const Group &) = default;
+
+        explicit Group(int ngens, const std::vector<Rel> &rels = {}, std::string name = "G");
+
+        void set(const Rel &r);
+
+        [[nodiscard]] int get(int a, int b) const;
+
+        [[nodiscard]] std::vector<Rel> rels() const;
+
+        [[nodiscard]] SubGroup subgroup(const std::vector<int> &gens) const;
+
+        [[nodiscard]] Group product(const Group &other) const;
+
+        [[nodiscard]] Group power(int p) const;
+
+        [[nodiscard]] Cosets solve(const std::vector<int> &sub_gens = {}) const;
+    };
+
+    struct SubGroup : public Group {
+        std::vector<int> gen_map;
+        const Group &parent;
+
+        SubGroup(const Group &parent, std::vector<int> gen_map);
+    };
+
+    Group operator*(const Group &g, const Group &h);
+
+    Group operator^(const Group &g, int p);
+}

include/tc/cosets.hpp

@@ -1,103 +0,0 @@
-#pragma once
-
-#include <vector>
-
-namespace tc {
-    class Path;
-
-    class Cosets {
-    private:
-        std::vector<int> data;
-        size_t _rank;
-
-    public:
-        Cosets(const Cosets &) = default;
-
-        explicit Cosets(size_t rank) : _rank(rank) {}
-
-        void add_row() {
-            data.resize(data.size() + rank(), -1);
-        }
-
-        void put(int coset, int gen, int target) {
-            data[coset * rank() + gen] = target;
-            data[target * rank() + gen] = coset;
-        }
-
-        [[nodiscard]] int get(int coset, int gen) const {
-            return data[coset * rank() + gen];
-        }
-
-        [[nodiscard]] size_t rank() const {
-            return _rank;
-        }
-
-        [[nodiscard]] size_t order() const {
-            if (!_rank) return 0;
-            return data.size() / _rank;
-        }
-
-        Path path() const;
-    };
-
-    class Path {
-    private:
-        friend class Cosets;
-
-        std::vector<unsigned int> source;
-        std::vector<unsigned int> gen;
-        size_t _rank;
-        size_t _order;
-
-        explicit Path(size_t rank, size_t order) : _rank(rank), _order(order), source(order), gen(order) {}
-
-    public:
-        [[nodiscard]] size_t rank() const {
-            return _rank;
-        }
-
-        [[nodiscard]] size_t order() const {
-            return _order;
-        }
-
-        template<class T, class F>
-        std::vector<T> walk(const T &start, const F &op) const {
-            std::vector<T> res;
-            res.reserve(order());
-            res.push_back(start);
-
-            for (size_t i = 1; i < order(); ++i) {
-                auto val = op(res[source[i]], gen[i]);
-                res.push_back(val);
-            }
-
-            return res;
-        }
-
-        template<class T, class E, class F>
-        std::vector<T> walk(const T &start, const E &gens, const F &op) const {
-            return walk(start, [&](const T &s, const int g) {
-                return op(s, gens[g]);
-            });
-        }
-    };
-
-    Path Cosets::path() const {
-        Path res(rank(), order());
-        std::vector<bool> set(order());
-
-        for (int coset = 0; coset < order(); ++coset) {
-            for (int gen = 0; gen < rank(); ++gen) {
-                int target = get(coset, gen);
-
-                if (!set[target]) {
-                    res.source[target] = coset;
-                    res.gen[target] = gen;
-                    set[target] = true;
-                }
-            }
-        }
-
-        return res;
-    }
-}

include/tc/group.hpp

@@ -1,174 +0,0 @@
-#pragma once
-
-#include <string>
-#include <sstream>
-#include <numeric>
-
-#include <Eigen/Eigen>
-
-namespace tc {
-    template<class T>
-    std::string stringify(const T &vec) {
-        std::stringstream ss;
-        ss << "[" << vec.transpose() << "]";
-        return ss.str();
-    }
-}
-
-namespace tc {
-    using Symbol = Eigen::Vector<unsigned int, Eigen::Dynamic>;
-
-    using MatrixXui = Eigen::Matrix<unsigned int, Eigen::Dynamic, Eigen::Dynamic>;
-    using ArrayXui = Eigen::Array<unsigned int, Eigen::Dynamic, Eigen::Dynamic>;
-
-    /// A Coxeter Matrix
-    class Group : public MatrixXui {
-    public:
-        using Base = MatrixXui;
-
-        std::string name = "G";
-        Symbol gens;
-
-        explicit Group(size_t rank) : Base(rank, rank), gens(rank) {
-            for (Eigen::Index i = 0; i < rank; ++i) {
-                gens(i) = i;
-            }
-        }
-
-        [[nodiscard]] size_t rank() const {
-            return rows();
-        }
-    };
-
-    Group subgroup(const Group &group, const Symbol &gens) {
-        size_t rank = group.size();
-        size_t srank = gens.size();
-
-        Group res(srank);
-        res.name = group.name + ":" + stringify(gens);
-        res.gens = gens;
-
-        for (Eigen::Index i = 0; i < srank; ++i) {
-            for (Eigen::Index j = 0; j < srank; ++j) {
-                res(i, j) = group(gens[i], gens[j]);
-            }
-        }
-
-        return res;
-    }
-
-    unsigned int factorial(unsigned int n) {
-        unsigned int res = 1;
-        for (int i = 1; i <= n; ++i) {
-            res *= i;
-        }
-        return res;
-    }
-
-    unsigned int choose(unsigned int n, unsigned int k) {
-        return factorial(n) / factorial(k) / factorial(n - k);
-    }
-
-    using SubGroups = std::vector<Group>;
-
-    SubGroups subgroups(const Group &group, size_t srank) {
-        size_t rank = group.rank();
-
-        std::vector<bool> mask(rank, false);
-        std::fill_n(mask.begin(), srank, true);
-
-        SubGroups res;
-        res.reserve(choose(rank, srank));
-
-        Symbol row(srank);
-        do {
-            for (int j = 0, k = 0; j < rank; ++j) {
-                if (mask[j])
-                    row(k++) = j;
-            }
-            res.push_back(subgroup(group, row));
-        } while (std::prev_permutation(mask.begin(), mask.end()));
-
-        return res;
-    }
-
-    /**
-     * Determine which of g_gens are the correct names for sg_gens within the current context
-     */
-    Symbol recontext_gens(
-        size_t rank,
-        Symbol g_gens,
-        Symbol sg_gens
-    ) {
-        std::sort(g_gens.begin(),  g_gens.end());
-        std::sort(sg_gens.begin(),  sg_gens.end());
-
-        int inv_gen_map[rank];
-        for (int i = 0; i < g_gens.size(); ++i) {
-            inv_gen_map[g_gens[i]] = i;
-        }
-
-        Symbol s_sg_gens(sg_gens.size());
-        for (int i = 0; i < sg_gens.size(); ++i) {
-            s_sg_gens[i] = inv_gen_map[sg_gens[i]];
-        }
-        std::sort(s_sg_gens.begin(),  s_sg_gens.end());
-
-        return s_sg_gens;
-    }
-
-    /**
-     * Create a named coxeter matrix from a simplified schlafli symbol
-     */
-    Group schlafli(const Symbol &mults, const std::string &name) {
-        size_t rank = mults.size() + 1;
-
-        Group res(rank);
-        res.name = name;
-
-        res.fill(2);
-        res.diagonal().fill(1);
-        res.topRightCorner(rank - 1, rank - 1).diagonal() << mults;
-        res.bottomLeftCorner(rank - 1, rank - 1).diagonal() << mults;
-
-        return res;
-    }
-
-    /**
-     * Create a coxeter matrix from a simplified schlafli symbol.
-     */
-    Group schlafli(const Symbol &mults) {
-        return schlafli(mults, stringify(mults));
-    }
-
-    Group product(const Group &g, const Group &h) {
-        Group res(g.rank() + h.rank());
-        res.name = g.name + "*" + h.name;
-
-        res.fill(2);
-
-        Eigen::Index off = 0;
-        res.block(off, off, g.rank(), g.rank()) << g.array() + off;
-        off += (Eigen::Index) g.rank();
-
-        res.block(off, off, h.rank(), h.rank()) << h.array() + off;
-        off += (Eigen::Index) h.rank();
-
-        return res;
-    }
-
-    Group power(const Group &g, size_t p) {
-        Group res(g.rank() * p);
-
-        res.name = g.name + "^" + std::to_string(p);
-
-        res.fill(2);
-
-        for (Eigen::Index k = 0; k < p; ++k) {
-            auto off = (Eigen::Index) g.rank() * k;
-            res.block(off, off, g.rank(), g.rank()) << g.array() + off;
-        }
-
-        return res;
-    }
-}

include/tc/groups.hpp

@@ -1,148 +1,70 @@
 #pragma once
 
-#include "group.hpp"
+#include "core.hpp"
 
-namespace tc::group {
+
+namespace tc {
     /**
-     * Universal Coxeter Group
+     * Construct a group from a (simplified) Schlafli Symbol of the form [a, b, ..., c]
+     * @param mults: The sequence of multiplicites between adjacent generators.
      */
-    Group U(size_t rank) {
+    Group schlafli(const std::vector<int> &mults, const std::string &name);
-        std::string name = "U(" + std::to_string(rank) + ")";
-
-        Group res(rank);
-        res.name = name;
-        res.fill(2);
-        return res;
-    }
 
     /**
-     * Simplex
+     * Construct a group from a (simplified) Schlafli Symbol of the form [a, b, ..., c]
+     * @param mults: The sequence of multiplicites between adjacent generators.
      */
-    Group A(size_t rank) {
+    Group schlafli(const std::vector<int> &mults);
-        std::string name = "A(" + std::to_string(rank) + ")";
 
-        if (rank == 0) {
+    namespace group {
-            Group res(rank);
+        /**
-            res.name = name;
+         * Simplex
-            return res;
+         */
-        }
+        Group A(int dim);
 
-        tc::Symbol symbol(rank - 1);
+        /**
-        symbol.fill(3);
+         * Cube, Orthoplex
+         */
+        Group B(int dim);
 
-        return schlafli(symbol, name);
+        /**
-    }
+         * Demicube, Orthoplex
+         */
+        Group D(int dim);
 
-    /**
+        /**
-     * Cube, Orthoplex
+         * E groups
-     */
+         */
-    Group B(size_t rank) {
+        Group E(int dim);
-        std::string name = "B(" + std::to_string(rank) + ")";
 
-        tc::Symbol symbol(rank - 1);
+        /**
-        symbol.fill(3);
+         * 24 Cell
-        symbol(0) = 4;
+         */
+        Group F4();
 
-        return schlafli(symbol, name);
+        /**
-    }
+         * Hexagon
+         */
+        Group G2();
 
-    /**
+        /**
-     * Demicube, Orthoplex
+         * Icosahedron
-     */
+         */
-    Group D(size_t rank) {
+        Group H(int dim);
-        std::string name = "D(" + std::to_string(rank) + ")";
 
-        tc::Symbol symbol(rank - 1);
+        /**
-        symbol.fill(3);
+         * Polygonal
-        symbol((Eigen::Index) rank - 2) = 2;
+         */
+        Group I2(int n);
 
-        Group g = schlafli(symbol, name);
+        /**
-        g(1, (Eigen::Index) rank - 1) = 3;
+         * Toroidal. I2(n) * I2(m)
-        g((Eigen::Index) rank - 1, 1) = 3;
+         */
+        Group T(int n, int m);
 
-        return g;
+        /**
-    }
+         * Toroidal. T(n, n)
-
+         */
-    /**
+        Group T(int n);
-     * E groups
-     */
-    Group E(size_t rank) {
-        std::string name = "E(" + std::to_string(rank) + ")";
-
-        tc::Symbol symbol(rank - 1);
-        symbol.fill(3);
-        symbol((Eigen::Index) rank - 2) = 2;
-
-        Group g = schlafli(symbol, name);
-        g(2, (Eigen::Index) rank - 1) = 3;
-        g((Eigen::Index) rank - 1, 2) = 3;
-
-        return g;
-    }
-
-    /**
-     * 24 Cell
-     */
-    Group F4() {
-        tc::Symbol symbol(3);
-        symbol << 3, 4, 3;
-        return schlafli(symbol, "F4");
-    }
-
-    /**
-     * Hexagon
-     */
-    Group G2() {
-        tc::Symbol symbol(1);
-        symbol << 6;
-        return schlafli(symbol, "G2");
-    }
-
-    /**
-     * Icosahedron
-     */
-    Group H(size_t rank) {
-        std::string name = "H(" + std::to_string(rank) + ")";
-
-        tc::Symbol symbol(rank - 1);
-        symbol.fill(3);
-        symbol(0) = 5;
-
-        return schlafli(symbol, name);
-    }
-
-    /**
-     * Polygonal
-     */
-    Group I2(unsigned int n) {
-        std::string name = "I2(" + std::to_string(n) + ")";
-
-        tc::Symbol symbol(1);
-        symbol << n;
-        return schlafli(symbol, name);
-    }
-
-    /**
-     * Toroidal. I2(n) * I2(m)
-     */
-    Group T(unsigned int n, unsigned int m) {
-        std::string name = "T(" + std::to_string(n) + "," + std::to_string(m) + ")";
-
-        tc::Symbol symbol(3);
-        symbol << n, 2, m;
-
-        return schlafli(symbol, name);
-    }
-
-    /**
-     * Toroidal. T(n, n)
-     */
-    Group T(unsigned int n) {
-        std::string name = "T(" + std::to_string(n) + ")";
-
-        tc::Symbol symbol(3);
-        symbol << n, 2, n;
-
-        return schlafli(symbol, name);
     }
 }

include/tc/solver.hpp

@@ -1,240 +0,0 @@
-#pragma once
-
-#include <algorithm>
-#include <array>
-#include <memory>
-#include <vector>
-#include <queue>
-
-#include "group.hpp"
-#include "cosets.hpp"
-
-namespace {
-    struct Row {
-        int gnr;
-        int *lst;
-    };
-
-    struct Table {
-    private:
-    public:
-        int i, j, mult;
-
-        std::vector<Row> rows;
-
-    public:
-        explicit Table(int i, int j, int mult) :
-            i(i), j(j), mult(mult) {
-        }
-    };
-
-    template<class T, size_t BlockSize = 4096>
-    class BlockAllocator {
-        /// 4096 seems to be the best (on my machine anway) from profiling.
-    private:
-        int block = 0;
-        int next = 0;
-        std::vector<T *> data = {build()};
-
-        T *build() {
-            T *blk = new T[BlockSize];
-            std::fill_n(blk, BlockSize, 0);
-            return blk;
-        }
-
-    public:
-        T *operator()() {
-            if (next >= BlockSize) {
-                data.push_back(build());
-                block++;
-                next = 0;
-            }
-
-            return &data[block][next++];
-        }
-
-        ~BlockAllocator() {
-            for (auto &blk: data) {
-                delete[] blk;
-            }
-        }
-    };
-
-    class Tables {
-    private:
-        int *null_lst_ptr = new int;
-        BlockAllocator<int> alloc;
-
-        std::vector<std::shared_ptr<Table>> tables;
-        std::vector<std::vector<std::shared_ptr<Table>>> deps;
-
-        size_t _rank;
-        size_t _rels;
-
-    public:
-        explicit Tables(const tc::Group &group) : _rank(group.rank()), _rels(rank() * (rank() + 1) / 2 - rank()) {
-            deps.resize(rank());
-            for (int i = 0; i < rank() - 1; ++i) {
-                for (int j = i + 1; j < rank(); ++j) {
-                    auto table = std::make_shared<Table>(i, j, group(i, j));
-                    tables.push_back(table);
-                    deps[i].push_back(table);
-                    deps[j].push_back(table);
-                }
-            }
-        }
-
-        [[nodiscard]] size_t rank() const {
-            return _rank;
-        }
-
-        [[nodiscard]] size_t rels() const {
-            return _rels;
-        }
-
-        void add_row() {
-            // std::vector already does block allocation.
-            for (const auto &table: tables) {
-                table->rows.emplace_back();
-            }
-        }
-
-        void initialize(int target, const tc::Cosets &cosets) {
-            for (auto &table: tables) {
-                Row &row = table->rows[target];
-
-                if (row.lst == nullptr) {
-                    if (cosets.get(target, table->i) != target and
-                        cosets.get(target, table->j) != target) {
-                        row.lst = alloc();
-                        row.gnr = 0;
-                    } else {
-                        row.lst = null_lst_ptr;
-                        row.gnr = -1;
-                    }
-                }
-            }
-        }
-
-        ~Tables() {
-            delete null_lst_ptr;
-        }
-
-        void learn(int coset, int gen, int target, const tc::Cosets &cosets, std::priority_queue<size_t> &facts) {
-            if (target == coset) {
-                for (auto &table: deps[gen]) {
-                    Row &target_row = table->rows[target];
-
-                    if (target_row.lst == nullptr) {
-                        target_row.gnr = -1;
-                    }
-                }
-            }
-
-            for (auto &table: deps[gen]) {
-                Row &target_row = table->rows[target];
-                Row &coset_row = table->rows[coset];
-
-                if (target_row.lst == nullptr) {
-                    target_row.lst = coset_row.lst;
-                    target_row.gnr = coset_row.gnr + 1;
-
-                    if (coset_row.gnr < 0) {
-                        target_row.gnr -= 2;
-                    }
-
-                    if (target_row.gnr == table->mult) {
-                        // forward learn
-                        int lst = *target_row.lst;
-                        int gen_ = (table->i == gen) ? table->j : table->i;
-                        facts.push(lst * rank() + gen_);
-                    } else if (target_row.gnr == -table->mult) {
-                        // stationary learn
-                        int gen_ = (table->i == gen) ? table->j : table->i;
-                        facts.push(target * rank() + gen_);
-                    } else if (target_row.gnr == table->mult - 1) {
-                        // determined family
-                        *target_row.lst = target;
-                    }
-                }
-            }
-        }
-    };
-}
-
-namespace tc {
-    /**
-     * Assumes that g is a coxeter group - that is, self-adjoint and the diagonal is 2.
-     */
-    tc::Cosets solve(const Group &group, const Symbol &s_gens) {
-        size_t rank = group.rank();
-
-        tc::Cosets cosets(rank);
-        cosets.add_row();
-
-        if (rank == 0) {
-            return cosets;
-        }
-
-        for (unsigned int gen: s_gens) {
-            if (gen < rank)
-                cosets.put(0, gen, 0);
-        }
-
-        Tables tables(group);
-        tables.add_row();
-        tables.initialize(0, cosets);
-
-        std::priority_queue<size_t> facts;
-
-        for (int coset = 0; coset < cosets.order(); coset++) {
-            for (int gen = 0; gen < rank; ++gen) {
-                if (cosets.get(coset, gen) >= 0) continue; // todo vector<bool> set
-
-                int target = cosets.order();
-                cosets.add_row();
-                tables.add_row();
-
-                facts.push(coset * rank + gen);
-
-                // todo nothing before the current coset will be used.
-                //  delete all table rows using old cosets to free memory early.
-                //  probably some unrolled linked list would be good; just drop
-                //  old blocks.
-
-                while (!facts.empty()) {
-                    int fact_idx = facts.top();
-                    facts.pop();
-
-                    int coset_ = fact_idx / rank;
-                    int gen_ = fact_idx % rank;
-
-                    if (cosets.get(coset_, gen_) != -1)
-                        continue;
-
-                    cosets.put(coset_, gen_, target);
-
-                    tables.learn(coset_, gen_, target, cosets, facts);
-                }
-
-                tables.initialize(target, cosets);
-            }
-        }
-
-        return cosets;
-    }
-
-/**
- * Solve the cosets generated by sg_gens within the subgroup generated by g_gens of the group context
- */
-    Cosets solve(
-        const Group &context,
-        const Symbol &g_gens,
-        const Symbol &sg_gens
-    ) {
-        const Symbol &proper_sg_gens = recontext_gens(context.rank(), g_gens, sg_gens);
-        const Group &group = subgroup(context, g_gens);
-
-        return solve(group, proper_sg_gens);
-    }
-}

src/core.cpp

@@ -0,0 +1,165 @@
+#include "tc/core.hpp"
+
+#include <utility>
+#include <sstream>
+#include <algorithm>
+
+namespace tc {
+    Action::Action(int from_idx, int gen)
+        : from_idx(from_idx), gen(gen) {
+    }
+
+    Path::Path(const tc::Group &context)
+        : context(context), path() {
+    }
+
+    void Path::add_row() {
+        path.resize(path.size() + 1);
+    }
+
+    Action Path::get(int to_idx) const {
+        return path[to_idx];
+    }
+
+    void Path::put(int from_idx, int gen, int to_idx) {
+        path[to_idx] = Action(from_idx, gen);
+    }
+
+    size_t Path::size() const {
+        return path.size();
+    }
+
+    Cosets::Cosets(const tc::Group &context)
+        : path(context), context(context), data() {
+    }
+
+    void Cosets::add_row() {
+        data.resize(data.size() + context.ngens, -1);
+        path.add_row();
+    }
+
+    void Cosets::put(int coset, int gen, int target) {
+        data[coset * context.ngens + gen] = target;
+        data[target * context.ngens + gen] = coset;
+
+        if (path.get(target).from_idx == -1) {
+            path.put(coset, gen, target);
+        }
+    }
+
+    void Cosets::put(int idx, int target) {
+        int coset = idx / context.ngens;
+        int gen = idx % context.ngens;
+        data[idx] = target;
+        data[target * context.ngens + gen] = coset;
+
+        if (path.get(target).from_idx == -1) {
+            path.put(coset, gen, target);
+        }
+    }
+
+    int Cosets::get(int coset, int gen) const {
+        return data[coset * context.ngens + gen];
+    }
+
+    int Cosets::get(int idx) const {
+        return data[idx];
+    }
+
+    size_t Cosets::size() const {
+        return path.size();
+    }
+
+    Rel::Rel(int a, int b, int m)
+        : gens({a, b}), mult(m) {
+    }
+
+    Rel Rel::shift(int off) const {
+        return Rel(gens[0] + off, gens[1] + off, mult);
+    }
+
+    Group::Group(int ngens, const std::vector<Rel> &rels, std::string name)
+        : ngens(ngens), name(std::move(name)) {
+        _mults.resize(ngens);
+
+        for (auto &mult : _mults) {
+            mult.resize(ngens, 2);
+        }
+
+        for (const auto &rel : rels) {
+            set(rel);
+        }
+    }
+
+    void Group::set(const Rel &r) {
+        _mults[r.gens[0]][r.gens[1]] = r.mult;
+        _mults[r.gens[1]][r.gens[0]] = r.mult;
+    }
+
+    int Group::get(int a, int b) const {
+        return _mults[a][b];
+    }
+
+    std::vector<Rel> Group::rels() const {
+        std::vector<Rel> res;
+        for (int i = 0; i < ngens - 1; ++i) {
+            for (int j = i + 1; j < ngens; ++j) {
+                res.emplace_back(i, j, get(i, j));
+            }
+        }
+        return res;
+    }
+
+    SubGroup Group::subgroup(const std::vector<int> &gens) const {
+        return SubGroup(*this, gens);
+    }
+
+    Group Group::product(const Group &other) const {
+        std::stringstream ss;
+        ss << name << "*" << other.name;
+
+        Group g(ngens + other.ngens, rels(), ss.str());
+
+        for (const auto &rel : other.rels()) {
+            g.set(rel.shift(ngens));
+        }
+
+        return g;
+    }
+
+    Group Group::power(int p) const {
+        std::stringstream ss;
+        ss << name << "^" << p;
+
+        Group g(ngens * p, {}, ss.str());
+        for (const auto &rel : rels()) {
+            for (int off = 0; off < g.ngens; off += ngens) {
+                g.set(rel.shift(off));
+            }
+        }
+
+        return g;
+    }
+
+    SubGroup::SubGroup(const Group &parent, std::vector<int> gen_map)
+        : Group(gen_map.size()), parent(parent) {
+
+        std::sort(gen_map.begin(), gen_map.end());
+        this->gen_map = gen_map;
+
+        for (size_t i = 0; i < gen_map.size(); ++i) {
+            for (size_t j = 0; j < gen_map.size(); ++j) {
+                int mult = parent.get(gen_map[i], gen_map[j]);
+                set(Rel(i, j, mult));
+            }
+        }
+    }
+
+    Group operator*(const Group &g, const Group &h) {
+        return g.product(h);
+    }
+
+    Group operator^(const Group &g, int p) {
+        return g.power(p);
+    }
+}

src/groups.cpp

@@ -0,0 +1,120 @@
+#include "tc/groups.hpp"
+
+#include <sstream>
+
+namespace tc {
+    Group schlafli(const std::vector<int> &mults, const std::string &name) {
+        int ngens = (int) mults.size() + 1;
+
+        Group g(ngens, {}, name);
+
+        for (int i = 0; i < (int) mults.size(); i++) {
+            g.set(Rel(i, i + 1, mults[i]));
+        }
+
+        return g;
+    }
+
+    Group schlafli(const std::vector<int> &mults) {
+        std::stringstream ss;
+        ss << "[";
+        if (!mults.empty()) {
+            for (size_t i = 0; i < mults.size() - 1; ++i) {
+                ss << mults[i] << ",";
+            }
+            ss << mults.back();
+        }
+        ss << "]";
+
+        return schlafli(mults, ss.str());
+    }
+
+    namespace group {
+        Group A(const int dim) {
+            std::stringstream ss;
+            ss << "A(" << dim << ")";
+
+            if (dim == 0)
+                return Group(0, {}, ss.str());
+
+            const std::vector<int> &mults = std::vector<int>(dim - 1, 3);
+
+            return schlafli(mults, ss.str());
+        }
+
+        Group B(const int dim) {
+            std::stringstream ss;
+            ss << "B(" << dim << ")";
+
+            std::vector<int> mults(dim - 1, 3);
+            mults[0] = 4;
+
+            return schlafli(mults, ss.str());
+        }
+
+        Group D(const int dim) {
+            std::stringstream ss;
+            ss << "D(" << dim << ")";
+
+            std::vector<int> mults(dim - 1, 3);
+            mults[dim - 2] = 2;
+
+            Group g = schlafli(mults, ss.str());
+            g.set(Rel(1, dim - 1, 3));
+
+            return g;
+        }
+
+        Group E(const int dim) {
+            std::stringstream ss;
+            ss << "E(" << dim << ")";
+
+            std::vector<int> mults(dim - 1, 3);
+            mults[dim - 2] = 2;
+
+            Group g = schlafli(mults, ss.str());
+            g.set(Rel(2, dim - 1, 3));
+
+            return g;
+        }
+
+        Group F4() {
+            return schlafli({3, 4, 3}, "F4");
+        }
+
+        Group G2() {
+            return schlafli({6}, "G2");
+        }
+
+        Group H(const int dim) {
+            std::stringstream ss;
+            ss << "H(" << dim << ")";
+
+            std::vector<int> mults(dim - 1, 3);
+            mults[0] = 5;
+
+            return schlafli(mults, ss.str());
+        }
+
+        Group I2(const int n) {
+            std::stringstream ss;
+            ss << "I2(" << n << ")";
+
+            return schlafli({n}, ss.str());
+        }
+
+        Group T(const int n, const int m) {
+            std::stringstream ss;
+            ss << "T(" << n << "," << m << ")";
+
+            return schlafli({n, 2, m}, ss.str());
+        }
+
+        Group T(const int n) {
+            std::stringstream ss;
+            ss << "T(" << n << ")";
+
+            return schlafli({n, 2, n}, ss.str());
+        }
+    }
+}

src/solve.cpp

@@ -0,0 +1,189 @@
+#include "tc/core.hpp"
+
+#include <algorithm>
+
+namespace tc {
+    struct RelTablesRow {
+        int *gnrs;
+        int **lst_ptrs;
+
+        RelTablesRow(int N, int *gnrs, int **lst_ptrs) : gnrs(gnrs), lst_ptrs(lst_ptrs) {
+            for (int i = 0; i < N; i++) {
+                lst_ptrs[i] = nullptr;
+            }
+        }
+    };
+
+    struct RelTables {
+        static const int ROW_BLOCK_SIZE = 64;
+        std::vector<Rel> rels;
+        std::vector<RelTablesRow *> rows;
+        int start = 0;
+        int num_tables;
+        int buffer_rows = 0;
+
+        explicit RelTables(const std::vector<Rel> &rels)
+            : num_tables(rels.size()), rels(rels) {
+        }
+
+        void add_row() {
+            if (buffer_rows == 0) {
+                int *gnrs_alloc = new int[num_tables * RelTables::ROW_BLOCK_SIZE];
+                int **lst_ptrs_alloc = new int *[num_tables * RelTables::ROW_BLOCK_SIZE];
+                for (int i = 0; i < RelTables::ROW_BLOCK_SIZE; i++) {
+                    rows.push_back(
+                        new RelTablesRow(num_tables, &gnrs_alloc[i * num_tables], &lst_ptrs_alloc[i * num_tables]));
+                }
+                buffer_rows = RelTables::ROW_BLOCK_SIZE;
+            }
+
+            buffer_rows--;
+        }
+
+        void del_rows_to(int idx) {
+            const int del_to = (idx / RelTables::ROW_BLOCK_SIZE) * RelTables::ROW_BLOCK_SIZE;
+            for (int i = start; i < del_to; i += RelTables::ROW_BLOCK_SIZE) {
+                delete[] rows[i]->gnrs;
+                delete[] rows[i]->lst_ptrs;
+                for (int j = 0; j < RelTables::ROW_BLOCK_SIZE; j++) {
+                    delete rows[i + j];
+                }
+                start += RelTables::ROW_BLOCK_SIZE;
+            }
+        }
+
+        ~RelTables() {
+            while (start < rows.size()) {
+                delete[] rows[start]->gnrs;
+                delete[] rows[start]->lst_ptrs;
+                for (int j = 0; j < RelTables::ROW_BLOCK_SIZE; j++) {
+                    delete rows[start + j];
+                }
+                start += RelTables::ROW_BLOCK_SIZE;
+            }
+        }
+    };
+
+    Cosets Group::solve(const std::vector<int> &sub_gens) const {
+        Cosets cosets(*this);
+        cosets.add_row();
+
+        if (ngens == 0) {
+            return cosets;
+        }
+
+        for (int g : sub_gens) {
+            if (g < ngens)
+                cosets.put(0, g, 0);
+        }
+
+        RelTables rel_tables(rels());
+        std::vector<std::vector<int>> gen_map(ngens);
+        int rel_idx = 0;
+        for (Rel m : rels()) {
+            gen_map[m.gens[0]].push_back(rel_idx);
+            gen_map[m.gens[1]].push_back(rel_idx);
+            rel_idx++;
+        }
+
+        int null_lst_ptr;
+        rel_tables.add_row();
+        RelTablesRow &row = *(rel_tables.rows[0]);
+        for (int table_idx = 0; table_idx < rel_tables.num_tables; table_idx++) {
+            Rel &ti = rel_tables.rels[table_idx];
+
+            if (cosets.get(ti.gens[0]) + cosets.get(ti.gens[1]) == -2) {
+                row.lst_ptrs[table_idx] = new int;
+                row.gnrs[table_idx] = 0;
+            } else {
+                row.lst_ptrs[table_idx] = &null_lst_ptr;
+                row.gnrs[table_idx] = -1;
+            }
+        }
+
+        int idx = 0;
+        int coset, gen, target, fact_idx, lst, gen_;
+        while (true) {
+            while (idx < cosets.data.size() and cosets.get(idx) >= 0)
+                idx++;
+
+            if (idx == cosets.data.size()) {
+                rel_tables.del_rows_to(idx / ngens);
+                break;
+            }
+
+            target = cosets.size();
+            cosets.add_row();
+            rel_tables.add_row();
+
+            std::vector<int> facts;
+            facts.push_back(idx);
+
+            coset = idx / ngens;
+            gen = idx % ngens;
+
+            rel_tables.del_rows_to(coset);
+
+            RelTablesRow &target_row = *(rel_tables.rows[target]);
+            while (!facts.empty()) {
+                fact_idx = facts.back();
+                facts.pop_back();
+
+                if (cosets.get(fact_idx) != -1)
+                    continue;
+
+                cosets.put(fact_idx, target);
+
+                coset = fact_idx / ngens;
+                gen = fact_idx % ngens;
+
+                if (target == coset)
+                    for (int table_idx : gen_map[gen])
+                        if (target_row.lst_ptrs[table_idx] == nullptr)
+                            target_row.gnrs[table_idx] = -1;
+
+                RelTablesRow &coset_row = *(rel_tables.rows[coset]);
+                for (int table_idx : gen_map[gen]) {
+                    if (target_row.lst_ptrs[table_idx] == nullptr) {
+                        Rel &ti = rel_tables.rels[table_idx];
+                        target_row.lst_ptrs[table_idx] = coset_row.lst_ptrs[table_idx];
+                        target_row.gnrs[table_idx] = coset_row.gnrs[table_idx] + 1;
+
+                        if (coset_row.gnrs[table_idx] < 0)
+                            target_row.gnrs[table_idx] -= 2;
+
+                        if (target_row.gnrs[table_idx] == ti.mult) {
+                            lst = *(target_row.lst_ptrs[table_idx]);
+                            delete target_row.lst_ptrs[table_idx];
+                            gen_ = ti.gens[(int) (ti.gens[0] == gen)];
+                            facts.push_back(lst * ngens + gen_);
+                        } else if (target_row.gnrs[table_idx] == -ti.mult) {
+                            gen_ = ti.gens[ti.gens[0] == gen];
+                            facts.push_back(target * ngens + gen_);
+                        } else if (target_row.gnrs[table_idx] == ti.mult - 1) {
+                            *(target_row.lst_ptrs[table_idx]) = target;
+                        }
+                    }
+                }
+
+                std::sort(facts.begin(), facts.end(), std::greater<>());
+            }
+
+            for (int table_idx = 0; table_idx < rel_tables.num_tables; table_idx++) {
+                Rel &ti = rel_tables.rels[table_idx];
+                if (target_row.lst_ptrs[table_idx] == nullptr) {
+                    if ((cosets.get(target, ti.gens[0]) != target) and
+                        (cosets.get(target, ti.gens[1]) != target)) {
+                        target_row.lst_ptrs[table_idx] = new int;
+                        target_row.gnrs[table_idx] = 0;
+                    } else {
+                        target_row.lst_ptrs[table_idx] = &null_lst_ptr;
+                        target_row.gnrs[table_idx] = -1;
+                    }
+                }
+            }
+        }
+
+        return cosets;
+    }
+}