WIP: add todos for next tasks

Templating the groups doesn't improve performance that much, but it severely hinders usability. Switch things to be dynamically sized so that different ranked groups are the same type.

Could also now separate headers and implementation... but will probably leave it for now.

.gitignore

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

.gitmodules

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

CMakeLists.txt

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

Testing/CMakeLists.txt

@@ -0,0 +1,22 @@
+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

@@ -0,0 +1,30 @@
+#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

@@ -0,0 +1,17 @@
+#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

@@ -0,0 +1,37 @@
+#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

@@ -1,5 +0,0 @@
-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

@@ -1,35 +0,0 @@
-#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

@@ -1,21 +0,0 @@
-#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

@@ -0,0 +1,11 @@
+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

@@ -0,0 +1,41 @@
+#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

@@ -0,0 +1,39 @@
+#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

@@ -0,0 +1,35 @@
+#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

@@ -0,0 +1,22 @@
+#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

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

include/tc/complex.hpp

@@ -0,0 +1,170 @@
+#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

@@ -1,165 +0,0 @@
-#pragma once
-
-#include <array>
-#include <functional>
-#include <vector>
-#include <string>
-
-namespace tc {
-    struct Action {
-        int from_idx = -1;
-        int gen = -1;
-
-        Action() = default;
-
-        Action(const Action &) = default;
-
-        Action(int from_idx, int gen);
-    };
-
-    struct Path {
-        std::vector<Action> path;
-
-        Path() = default;
-
-        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 C, class T, class E>
-        void walk(
-            C& res,
-            T start,
-            std::vector<E> gens,
-            std::function<T(const T &, const E &)> op
-        ) const {
-            size_t s = size();
-            res.reserve(s);
-            res.push_back(start);
-
-            for (int i = 1; i < s; ++i) {
-                auto &action = path[i];
-                auto &from = res.get(action.from_idx);
-                auto &val = gens[action.gen];
-                res.push_back(op(from,val));
-            }
-        }
-
-        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;
-            res.reserve(size());
-            res.push_back(start);
-
-            for (int i = 1; i < size(); ++i) {
-                auto &action = path[i];
-                auto &from = res[action.from_idx];
-                auto &val = gens[action.gen];
-                res.push_back(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;
-            res.reserve(size());
-            res.push_back(start);
-
-            for (int i = 1; i < 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 {
-        int ngens;
-        std::vector<int> data;
-        Path path;
-
-        Cosets(const Cosets &) = default;
-
-        explicit Cosets(int ngens);
-
-        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 SubGroup;
-
-    struct Group {
-        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

@@ -0,0 +1,103 @@
+#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

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

@@ -0,0 +1,240 @@
+#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

@@ -1,161 +0,0 @@
-#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) {
-    }
-
-    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(int ngens)
-        : ngens(ngens) {
-    }
-
-    void Cosets::add_row() {
-        data.resize(data.size() + ngens, -1);
-        path.add_row();
-    }
-
-    void Cosets::put(int coset, int gen, int target) {
-        data[coset * ngens + gen] = target;
-        data[target * 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 / ngens;
-        int gen = idx % ngens;
-        data[idx] = target;
-        data[target * 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 * 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

@@ -1,120 +0,0 @@
-#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

@@ -1,189 +0,0 @@
-#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(ngens);
-        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;
-    }
-}