COMP: Replace toddcox-faster by reimplemented toddcox

CMakeLists.txt

@@ -1,16 +1,19 @@
 cmake_minimum_required(VERSION 3.10)
 project(toddcox-faster)
 
-set(CMAKE_CXX_STANDARD 17)
+set(CMAKE_CXX_STANDARD 20)
 
 include(FetchContent)
 
-include(vendor/toddcox.cmake)
+include(vendor/fmt.cmake)
-include(vendor/yaml-cpp.cmake)
+include(vendor/peglib.cmake)
 include(vendor/glad.cmake)
 include(vendor/glfw.cmake)
 include(vendor/glm.cmake)
+include(vendor/gtest.cmake)
+include(vendor/yaml-cpp.cmake)
 
+add_subdirectory(tc)
 add_subdirectory(vis)
 
 add_subdirectory(examples)

tc/CMakeLists.txt

@@ -0,0 +1,17 @@
+add_library(tc
+    include/tc/core.hpp
+    include/tc/groups.hpp
+
+    src/cosets.cpp
+    src/group.cpp
+    src/groups.cpp
+    src/lang.cpp
+    src/solve.cpp
+    )
+target_link_libraries(tc peglib::peglib fmt::fmt)
+target_include_directories(tc PUBLIC include)
+
+add_library(tc::tc ALIAS tc)
+
+add_subdirectory(test)
+add_subdirectory(bench)

tc/bench/CMakeLists.txt

@@ -0,0 +1,5 @@
+add_executable(benchmark benchmark.cpp)
+target_link_libraries(benchmark PUBLIC tc fmt::fmt)
+
+add_executable(named named.cpp)
+target_link_libraries(named PUBLIC tc fmt::fmt)

tc/bench/benchmark.cpp

@@ -0,0 +1,138 @@
+#include <ctime>
+#include <string>
+#include <vector>
+
+#include <fmt/core.h>
+#include <fmt/ranges.h>
+
+#include <tc/core.hpp>
+#include <tc/groups.hpp>
+
+void bench(
+    const std::string &group_expr,
+    const std::string &symbol,
+    const std::vector<size_t> &gens,
+    const size_t bound = SIZE_MAX
+) {
+    tc::Group<> group = tc::coxeter(symbol);
+    
+    std::clock_t s = std::clock();
+    tc::Cosets<> cosets = group.solve(gens, bound);
+    std::clock_t e = std::clock();
+
+    auto time = (double) (e - s) / CLOCKS_PER_SEC;
+    size_t order = cosets.order();
+    auto cos_s = (size_t) (order / time);
+
+    bool complete = cosets.complete();
+
+    std::string name = fmt::format("{}/{}", group_expr, gens);
+    std::string row = fmt::format(
+        "{:>24},{:>10},{:>6},{:>8.3f}s,{:>10L}",
+        name, order, complete, time, cos_s
+    );
+    fmt::print("{}\n", row);
+}
+
+int main(int argc, char *argv[]) {
+    std::vector<std::string> args(argv + 1, argv + argc);
+
+    fmt::print("{:>24},{:>10},{:>6},{:>9},{:>10}\n", "NAME", "ORDER", "COMPL", "TIME", "COS/S");
+
+    // Finite Groups
+    
+    // A_n: 3 * `n-1`           ; n >= 1
+    bench("A_5", "3 * 4", {});
+    bench("A_6", "3 * 5", {});
+    bench("A_7", "3 * 6", {});
+    bench("A_8", "3 * 7", {});
+    // B_n: 4 3 * `n-2`         ; n >= 2
+    bench("B_5", "4 3 * 3", {});
+    bench("B_6", "4 3 * 4", {});
+    bench("B_7", "4 3 * 5", {});
+    bench("B_8", "4 3 * 6", {});
+    // D_n: 3 * [1 1 `n-3`]     ; n >= 4
+    bench("D_5", "3 * [1 1 2]", {});
+    bench("D_6", "3 * [1 1 3]", {});
+    bench("D_7", "3 * [1 1 4]", {});
+    bench("D_8", "3 * [1 1 5]", {});
+    // E_n: 3 * [1 2 `n-4`]     ; n >= 6
+    bench("E_6", "3 * [1 2 2]", {});
+    bench("E_7", "3 * [1 2 3]", {});
+//    bench("E_8", "3 * [1 2 4]", {}); // too big
+    // H_n: 5 3 * `n-2`         ; n >= 2
+    bench("H_3", "5 3 * 0", {});
+    bench("H_4", "5 3 * 1", {});
+    bench("H_5", "5 3 * 2", {});
+    // grid: `p` `q`            ; 2(p+q) > pq
+    // triangle: `p` `q` `r`    ; 1/p + 1/q + 1/r > 1
+
+    // Special Finite Groups
+    bench("F_4", "3 4 3", {});
+    bench("G_2", "6", {});
+    // I_2(p): `p`              ; p >= 2
+    bench("I_2(100)", "100", {});
+    bench("I_2(1000)", "1000", {});
+    // "Torus": `p` 2 `q`       ; p, q >= 2
+    bench("T(100)", "100 2 100", {});
+    bench("T(1000)", "1000 2 1000", {});
+
+    // Affine Groups
+    
+    // ~A_n: {3 * `n+1`}
+    bench("~A_5", "{3 * 6}", {}, 10'000'000);
+    bench("~A_6", "{3 * 7}", {}, 10'000'000);
+    bench("~A_7", "{3 * 8}", {}, 10'000'000);
+    bench("~A_8", "{3 * 9}", {}, 10'000'000);
+    // ~B_n: 4 3 * `n-3` 3 * [1 1]
+    bench("~B_5", "4 3 * 2 3 * [1 1]", {}, 10'000'000);
+    bench("~B_6", "4 3 * 3 3 * [1 1]", {}, 10'000'000);
+    bench("~B_7", "4 3 * 4 3 * [1 1]", {}, 10'000'000);
+    bench("~B_8", "4 3 * 5 3 * [1 1]", {}, 10'000'000);
+    // ~B_n: 4 3 * `n-2` 4
+    bench("~C_5", "4 3 * 3 4", {}, 10'000'000);
+    bench("~C_6", "4 3 * 4 4", {}, 10'000'000);
+    bench("~C_7", "4 3 * 5 4", {}, 10'000'000);
+    bench("~C_8", "4 3 * 6 4", {}, 10'000'000);
+    // ~D_n: 3 * [1 1] 3 * `n-4` 3 * [1 1]
+    bench("~D_5", "3 * [1 1] 3 * 1 3 * [1 1]", {}, 10'000'000);
+    bench("~D_6", "3 * [1 1] 3 * 2 3 * [1 1]", {}, 10'000'000);
+    bench("~D_7", "3 * [1 1] 3 * 3 3 * [1 1]", {}, 10'000'000);
+    bench("~D_8", "3 * [1 1] 3 * 4 3 * [1 1]", {}, 10'000'000);
+    // grid: `p` `q`            ; 2(p+q) = pq
+    // triangle: `p` `q` `r`    ; 1/p + 1/q + 1/r = 1
+    
+    // Special Affine Groups
+    bench("~I_1", "-", {}, 10'000'000);
+    bench("~E_6", "3 * [2 2 2]", {}, 10'000'000);
+    bench("~E_7", "3 * [1 3 3]", {}, 10'000'000);
+    bench("~E_8", "3 * [1 2 5]", {}, 10'000'000);
+//    bench("E_9",  "3 * [1 2 5]", {}, 10'000'000);  // ~E_8 == E_9
+    bench("~F_4", "3 4 3 3", {}, 10'000'000);
+    bench("~G_2", "6 3", {}, 10'000'000);
+
+    // Hyperbolic Groups
+    // grid: `p` `q`            ; 2(p+q) < pq
+    // triangle: `p` `q` `r`    ; 1/p + 1/q + 1/r < 1
+    
+    // Special Hyperbolic Groups
+    bench("-BH_3", "4 3 5", {}, 10'000'000);
+    bench("-K_3", "5 3 5", {}, 10'000'000);
+    bench("-J_3", "3 5 3", {}, 10'000'000);
+//    bench("~H_3", "3 5 3", {}, 10'000'000);  // -J_3 == ~H_3
+    bench("-DH_3", "5 3 * [1 1]", {}, 10'000'000);
+    bench("^AB_3", "{3 3 3 4}", {}, 10'000'000);
+    bench("^AH_3", "{3 3 3 5}", {}, 10'000'000);
+    bench("^BB_3", "{3 4 3 4}", {}, 10'000'000);
+    bench("^BH_3", "{3 4 3 5}", {}, 10'000'000);
+    bench("^HH_3", "{3 5 3 5}", {}, 10'000'000);
+    bench("-H_4", "5 3 3 3", {}, 10'000'000);
+//    bench("~H_4", "5 3 3 3", {}, 10'000'000);  // -H_4 == ~H_4 == H_5 
+//    bench("H_5", "5 3 3 3", {}, 10'000'000);
+    bench("-BH_4", "4 3 3 5", {}, 10'000'000);
+    bench("-K_4", "5 3 3 5", {}, 10'000'000);
+    bench("-DH_4", "5 3 3 * [1 1]", {}, 10'000'000);
+    bench("^AF_4", "{3 3 3 3 4}", {}, 10'000'000);
+    
+    return EXIT_SUCCESS;
+}

tc/bench/named.cpp

@@ -0,0 +1,149 @@
+#include <string>
+#include <vector>
+
+#include <fmt/core.h>
+#include <fmt/ranges.h>
+
+#include <tc/core.hpp>
+
+using namespace std::string_literals;
+
+template<typename Gen_>
+struct fmt::formatter<tc::Group<Gen_>> {
+    template<typename ParseContext>
+    constexpr auto parse(ParseContext &ctx) {
+        return ctx.begin();
+    }
+
+    template<typename FormatContext>
+    constexpr auto format(tc::Group<Gen_> const &g, FormatContext &ctx) {
+        auto gens = g.gens();
+
+        fmt::format_to(ctx.out(), "  | ");
+        for (const auto &gen: gens) {
+            fmt::format_to(ctx.out(), "{} ", gen);
+        }
+        fmt::format_to(ctx.out(), "\n");
+
+        for (int i = 0; i < gens.size(); ++i) {
+            auto u = gens[i];
+
+            fmt::format_to(ctx.out(), "{} | ", u);
+
+            for (int j = 0; j < gens.size(); ++j) {
+                auto v = gens[j];
+
+                if (i <= j) {
+                    fmt::format_to(ctx.out(), "{} ", g.get(u, v));
+                } else {
+                    fmt::format_to(ctx.out(), "  ");
+                }
+            }
+
+            fmt::format_to(ctx.out(), "\n");
+        }
+
+        return ctx.out();
+    }
+};
+
+template<typename Gen_>
+struct fmt::formatter<tc::Cosets<Gen_>> {
+    template<typename ParseContext>
+    constexpr auto parse(ParseContext &ctx) {
+        return ctx.begin();
+    }
+
+    template<typename FormatContext>
+    constexpr auto format(tc::Cosets<Gen_> const &c, FormatContext &ctx) {
+        auto gens = c.gens();
+
+        auto width = (size_t) (std::log10(c.order() + 1));
+
+        fmt::format_to(ctx.out(), "{:>{}} | ", "-", width);
+        for (const auto &gen: gens) {
+            fmt::format_to(ctx.out(), "{:>{}} ", gen, width);
+        }
+        fmt::format_to(ctx.out(), "\n");
+
+        for (size_t cos = 0; cos < c.order(); ++cos) {
+            fmt::format_to(ctx.out(), "{:>{}} | ", cos, width);
+
+            for (const auto &gen: gens) {
+                auto target = c.get(cos, gen);
+
+                fmt::format_to(ctx.out(), "{:>{}} ", target, width);
+            }
+
+            fmt::format_to(ctx.out(), "\n");
+        }
+
+        return ctx.out();
+    }
+};
+
+int main() {
+    tc::Group<char> group(4, {'r', 'g', 'b', 'y'});
+
+    group.set('r', 'g', 5);
+    group.set('g', 'b', 3);
+    group.set('b', 'y', 3);
+
+    fmt::print("{}\n", group);
+
+//    auto res = group.sub({'r', 'g', 'b', 'y'}).solve();
+    auto res = group.sub({'r', 'g'}).solve();
+    fmt::print("rank: {}, order: {}\n{}\n",
+               res.rank(), res.order(), res);
+
+//    tc::Path<> path(res);
+//    auto gens = group.gens();
+//    std::vector<std::string> words(path.order());
+//    path.walk(
+//        "-"s,
+//        [&](const std::string &cos, size_t gen) -> std::string {
+//            return cos + gens[gen];
+//        },
+//        words.begin());
+
+    tc::Path<char> typed(res);
+    std::vector<std::string> words(typed.order());
+    typed.walk("-"s, std::plus<>(), words.begin());
+
+    fmt::print("words: {}\n", words);
+    fmt::print("size: {}\n", words.size());
+
+//    path.walk(
+//        words.begin(),
+//        "-",
+//        [&gens](const std::string &cos, size_t gen) {
+//            return cos + gens[gen];
+//        }
+//    );
+//    for (const auto &word: words) {
+//        fmt::print("'{}'\n", word);
+//    }
+
+//    auto sub = group.sub({'r', 'g', 'y'});
+//    fmt::print("{}\n", sub);
+//
+//    auto res = sub.solve({});
+//    fmt::print("res order: {}\n", res.order());
+//
+//    auto cos = sub.solve({'r', 'y'});
+
+//    fmt::print("order: {}\n", cos.order());
+//
+//    fmt::print("{}\n", cos);
+
+//    tc::Group<u_int8_t> group(4, {0, 1, 2, 3});
+//
+//    group.set(0, 1, 5);
+//    group.set(1, 2, 4);
+//    group.set(2, 3, 3);
+//
+//    show(group);
+//
+//    auto sub = group.sub({3, 2, 0, 1});
+//    show(sub);
+}

tc/include/tc/core.hpp

@@ -0,0 +1,306 @@
+#pragma once
+
+#include <vector> //todo clean up includes. lots of duplicate cstdint, cassert.
+#include <cstdint>
+#include <cassert>
+#include <tuple>
+
+#include <limits>
+
+namespace tc {
+    using Mult = u_int16_t;
+    constexpr Mult FREE = 0;
+
+    /**
+     * @brief Mapping from "global" generator names or objects to indexes used for value lookup. 
+     * @tparam Gen_ 
+     */
+    template<typename Gen_=void>
+    struct Index;
+
+    /**
+     * @brief Complete representation of a quotient group. Describes the action of each generator on each coset.
+     * @tparam Gen_ 
+     */
+    template<typename Gen_=void>
+    struct Cosets;
+
+    /** 
+     * @brief Manage the presentation of a Coxeter group and enforce constraints
+     * on the multiplicities of its relations.
+     * <ul>
+     *   <li>
+     *     <code>m_ij = 1</code> iff <code>i != j</code>
+     *   </li>
+     *   <li>
+     *     <code>m_ij = m_ji</code>
+     *   </li>
+     *   <li>
+     *     If <code>m_ij == inf</code> (<code>tc::FREE</code>) then no relation is imposed.
+     *   </li>
+     * </ul>
+     * @see
+     * <a href="https://en.wikipedia.org/wiki/Coxeter_group#Definition">Coxeter Group (Wikipedia)</a>
+     */
+    template<typename Gen_=void>
+    struct Group;
+
+    /**
+     * @brief Support generating values given a Cosets and transformation callback.
+     * @tparam Gen_ 
+     */
+    template<typename Gen_=void>
+    struct Path;
+
+    template<>
+    struct Index<> {
+        size_t operator()(size_t const &idx) const {
+            return idx;
+        }
+    };
+
+    template<typename Gen_>
+    struct Index {
+        using Gen = Gen_;
+
+        std::vector<Gen> _gens{};
+
+        explicit Index(std::vector<Gen> gens) : _gens(gens) {}
+
+        template<typename It>
+        Index(It begin, It end) : _gens(begin, end) {}
+
+        size_t operator()(Gen const &gen) const {
+            auto it = std::find(_gens.begin(), _gens.end(), gen);
+            assert(it != _gens.end());
+            return it - _gens.begin();
+        }
+    };
+
+    template<>
+    struct Cosets<> {
+        static constexpr size_t UNSET = std::numeric_limits<size_t>::max();
+
+    private:
+        size_t _rank;
+        size_t _order;
+        bool _complete;
+        std::vector<size_t> _data;
+
+    public:
+        Cosets(Cosets const &) = default;
+
+        Cosets(Cosets &&) noexcept = default;
+
+        ~Cosets() = default;
+
+        void set(size_t coset, size_t gen, size_t target);
+
+        [[nodiscard]] size_t get(size_t coset, size_t gen) const;
+
+        [[nodiscard]] bool isset(size_t coset, size_t gen) const;
+
+        [[nodiscard]] size_t rank() const;
+
+        [[nodiscard]] size_t order() const;
+
+        [[nodiscard]] bool complete() const;
+
+        [[nodiscard]] size_t size() const;
+
+        friend Group<>;  // only constructible via Group<>::solve
+
+    private:
+        explicit Cosets(size_t rank);
+
+        void add_row();
+
+        void set(size_t idx, size_t target);
+
+        [[nodiscard]] size_t get(size_t idx) const;
+
+        [[nodiscard]] bool isset(size_t idx) const;
+    };
+
+    template<>
+    struct Group<> {
+        using Rel = std::tuple<size_t, size_t, Mult>;
+
+    private:
+        size_t _rank;
+        std::vector<size_t> _mults;
+
+    public:
+        Group(Group const &) = default;
+
+        Group(Group &&) noexcept = default;
+
+        ~Group() = default;
+
+        explicit Group(size_t rank);
+
+        void set(size_t, size_t, Mult);
+
+        [[nodiscard]] Mult get(size_t, size_t) const;
+
+        [[nodiscard]] size_t rank() const;
+
+        [[nodiscard]] Group sub(std::vector<size_t> const &idxs) const;
+
+        [[nodiscard]] Cosets<> solve(std::vector<size_t> const &idxs = {}, size_t bound = SIZE_MAX) const;
+    };
+
+    template<>
+    struct Path<> {
+        using Action = std::tuple<size_t, size_t>;  // coset, gen
+
+    private:
+        std::vector<Action> _data;
+
+    public:
+        // todo be smarter about move semantics
+
+        explicit Path(Cosets<> const &cosets) : _data() {
+            _data.resize(cosets.order());
+
+            std::vector<bool> complete(cosets.order(), false);
+
+            complete[0] = true;
+
+            for (size_t cos = 0; cos < cosets.order(); ++cos) {
+                for (size_t gen = 0; gen < cosets.rank(); ++gen) {
+                    size_t tgt = cosets.get(cos, gen);
+                    if (!complete[tgt]) {
+                        _data[tgt] = {cos, gen};
+                        complete[tgt] = true;
+                    }
+                }
+            }
+        }
+
+        template<typename Elem, typename BinaryOp, std::random_access_iterator It>
+        void walk(Elem const &start, BinaryOp op, It out) {
+            // todo how to work for std::back_insert_iterator?
+            // todo how to work for non-default-constructible values?
+
+            out[0] = start;
+
+            for (int tgt = 1; tgt < _data.size(); ++tgt) {
+                auto [cos, gen] = _data[tgt];
+                out[tgt] = op(out[cos], gen);
+            }
+        }
+
+        template<typename Elem, typename BinaryOp, std::random_access_iterator It, std::random_access_iterator Gens>
+        void walk(Elem const &start, BinaryOp op, It out, Gens const &gens) {
+            out[0] = start;
+
+            for (int tgt = 1; tgt < _data.size(); ++tgt) {
+                auto [cos, gen] = _data[tgt];
+                out[tgt] = op(out[cos], gens[gen]);
+            }
+        }
+
+        [[nodiscard]] size_t order() const {
+            return _data.size();
+        }
+    };
+
+    template<typename Gen_>
+    struct Cosets : public Cosets<> {
+        using Gen = Gen_;
+
+    private:
+        Index<Gen> _index;
+
+        friend Path<Gen>;
+
+    public:
+        Cosets(Cosets<> g, std::vector<Gen> gens) : Cosets<>(g), _index(gens) {}
+
+        void set(size_t coset, Gen const &gen, size_t target) {
+            Cosets<>::set(coset, _index(gen), target);
+        }
+
+        [[nodiscard]] size_t get(size_t coset, Gen const &gen) const {
+            return Cosets<>::get(coset, _index(gen));
+        }
+
+        [[nodiscard]] bool isset(size_t coset, Gen const &gen) const {
+            return Cosets<>::isset(coset, _index(gen));
+        }
+
+        [[nodiscard]] std::vector<Gen> gens() const {
+            return _index._gens;
+        }
+
+    private:
+        Cosets(size_t rank, std::vector<Gen> gens) : Cosets<>(rank), _index(gens) {}
+    };
+
+    template<typename Gen_>
+    struct Group : public Group<> {
+        using Gen = Gen_;
+        using Rel = std::tuple<Gen, Gen, Mult>;
+
+    private:
+        Index<Gen> _index;
+
+    public:
+        Group(Group const &) = default;
+
+        Group(Group &&) noexcept = default;
+
+        Group(Group<> g, std::vector<Gen> gens) : Group<>(g), _index(gens) {}
+
+        Group(size_t rank, std::vector<Gen> gens) : Group<>(rank), _index(gens) {}
+
+        ~Group() = default;
+
+        void set(Gen const &u, Gen const &v, Mult m) {
+            Group<>::set(_index(u), _index(v), m);
+        }
+
+        [[nodiscard]] Mult get(Gen const &u, Gen const &v) const {
+            return Group<>::get(_index(u), _index(v));
+        }
+
+        [[nodiscard]] std::vector<Gen> gens() const {
+            return _index._gens;
+        }
+
+        [[nodiscard]] Group sub(std::vector<Gen> const &gens) const {
+            std::vector<size_t> idxs(gens.size());
+            std::transform(gens.begin(), gens.end(), idxs.begin(), _index);
+            return Group(Group<>::sub(idxs), gens);
+        }
+
+        [[nodiscard]] Cosets<Gen> solve(std::vector<Gen> const &gens = {}, size_t bound = SIZE_MAX) const {
+            std::vector<size_t> idxs(gens.size());
+            std::transform(gens.begin(), gens.end(), idxs.begin(), _index);
+
+            return Cosets<Gen>(Group<>::solve(idxs, bound), _index._gens);
+        }
+    };
+
+    template<typename Gen_>
+    struct Path : public Path<> {
+        using Gen = Gen_;
+
+    private:
+        Index<Gen> _index;
+
+    public:
+        // todo be smarter about move semantics
+        explicit Path(Cosets<Gen> const &cosets) : Path<>(cosets), _index(cosets._index) {}
+
+        template<typename T>
+        Path(Cosets<> const &cosets, T const &gens)
+            : Path<>(cosets), _index(gens.begin(), gens.end()) {}
+
+        template<typename Elem, typename BinaryOp, std::random_access_iterator It>
+        void walk(Elem const &start, BinaryOp op, It out) {
+            Path<>::walk(start, op, out, _index._gens.begin());
+        }
+    };
+}

tc/include/tc/groups.hpp

@@ -0,0 +1,24 @@
+#pragma once
+
+#include <vector>
+#include <string>
+
+#include <tc/core.hpp>
+
+namespace tc {
+    /**
+     * Construct a group from a (simplified) Schlafli Symbol of the form [a, b, ..., c]
+     * @param mults: The sequence of multiplicites between adjacent generators.
+     */
+    Group<> schlafli(const std::vector<unsigned int> &mults);
+
+    Group<> coxeter(const std::string &symbol);
+
+    Group<> vcoxeter(const std::string &symbol, const std::vector<unsigned int> &values);
+
+    template<typename ...Args>
+    Group<> coxeter(const std::string &symbol, const Args &... args) {
+        std::vector<unsigned int> values = {{args...}};
+        return vcoxeter(symbol, values);
+    }
+}

tc/src/cosets.cpp

@@ -0,0 +1,55 @@
+#include <tc/core.hpp>
+
+namespace tc {
+    Cosets<>::Cosets(size_t rank)
+        : _rank(rank), _order(0), _complete(false), _data() {}
+
+    void Cosets<>::set(size_t coset, size_t gen, size_t target) {
+        set(coset * rank() + gen, target);
+    }
+
+    [[nodiscard]] size_t Cosets<>::get(size_t coset, size_t gen) const {
+        return get(coset * rank() + gen);
+    }
+
+    [[nodiscard]] bool Cosets<>::isset(size_t coset, size_t gen) const {
+        return isset(coset * rank() + gen);
+    }
+
+    [[nodiscard]] size_t Cosets<>::rank() const {
+        return _rank;
+    }
+
+    [[nodiscard]] size_t Cosets<>::order() const {
+        return _order;
+    }
+
+    [[nodiscard]] bool Cosets<>::complete() const {
+        return _complete;
+    }
+
+    [[nodiscard]] size_t Cosets<>::size() const {
+        return _data.size();
+    }
+
+    void Cosets<>::add_row() {
+        _data.resize(_data.size() + rank(), UNSET);
+        _order++;
+    }
+
+    void Cosets<>::set(size_t idx, size_t target) {
+        size_t coset = idx / rank();
+        size_t gen = idx % rank();
+        _data[idx] = target;
+        _data[target * rank() + gen] = coset;
+    }
+
+    [[nodiscard]] size_t Cosets<>::get(size_t idx) const {
+        return _data[idx];
+    }
+
+    [[nodiscard]] bool Cosets<>::isset(size_t idx) const {
+        return get(idx) != UNSET;
+    }
+
+}

tc/src/group.cpp

@@ -0,0 +1,42 @@
+#include <tc/core.hpp>
+
+#include <cassert>
+
+namespace tc {
+    Group<>::Group(size_t rank) : _rank(rank), _mults(_rank * _rank, 2) {
+        for (int idx = 0; idx < rank; ++idx) {
+            set(idx, idx, 1);
+        }
+    }
+
+    void Group<>::set(size_t u, size_t v, Mult m) {
+        assert(u < rank());
+        assert(v < rank());
+
+        _mults[u * rank() + v] = m;
+        _mults[v * rank() + u] = m;
+    }
+
+    [[nodiscard]] Mult Group<>::get(size_t u, size_t v) const {
+        assert(u < rank());
+        assert(v < rank());
+
+        return _mults[u * rank() + v];
+    }
+
+    [[nodiscard]] size_t Group<>::rank() const {
+        return _rank;
+    }
+
+    [[nodiscard]] Group<> Group<>::sub(std::vector<size_t> const &idxs) const {
+        Group<> res(idxs.size());
+
+        for (int i = 0; i < idxs.size(); ++i) {
+            for (int j = i; j < idxs.size(); ++j) {
+                res.set(i, j, get(idxs[i], idxs[j]));
+            }
+        }
+
+        return res;
+    }
+}

tc/src/groups.cpp

@@ -0,0 +1,25 @@
+#include <tc/groups.hpp>
+
+#include <fmt/args.h>
+#include <fmt/core.h>
+#include <numeric>
+
+namespace tc {
+    Group<> schlafli(const std::vector<unsigned int> &mults) {
+        Group<> res(mults.size() + 1);
+        for (size_t i = 0; i < mults.size(); ++i) {
+            res.set(i, i + 1, mults[i]);
+        }
+        return res;
+    }
+
+    Group<> vcoxeter(const std::string &symbol, const std::vector<unsigned int> &values) {
+        fmt::dynamic_format_arg_store<fmt::format_context> ds;
+
+        for (const auto &value: values) {
+            ds.push_back(value);
+        }
+
+        return coxeter(fmt::vformat(symbol, ds));
+    }
+}

tc/src/lang.cpp

@@ -0,0 +1,299 @@
+#include <vector>
+#include <stack>
+#include <string>
+#include <cassert>
+
+#include <tc/core.hpp>
+#include <tc/groups.hpp>
+
+#include <peglib.h>
+
+#include <fmt/core.h>
+#include <fmt/ranges.h>
+#include <numeric>
+
+struct Graph {
+    size_t rank{};
+    std::vector<tc::Group<>::Rel> edges{};
+};
+
+struct Op {
+    enum Code {
+        LINK,
+        PUSH,
+        POP,
+        LOOP,
+        FREE,
+    };
+
+    Code code: 4;
+    unsigned int value: 12;
+
+    explicit Op(Code code, unsigned int value = 0)
+        : code(code), value(value) {}
+};
+
+template<>
+struct fmt::formatter<Op> {
+    template<typename ParseContext>
+    constexpr auto parse(ParseContext &ctx) {
+        return ctx.begin();
+    }
+
+    template<typename FormatContext>
+    auto format(const Op &op, FormatContext &ctx) {
+        switch (op.code) {
+            case Op::LINK:
+                return fmt::format_to(ctx.out(), "link({})",
+                                      (unsigned int) op.value);
+            case Op::PUSH:
+                return fmt::format_to(ctx.out(), "push()");
+            case Op::POP:
+                return fmt::format_to(ctx.out(), "pop()");
+            case Op::LOOP:
+                return fmt::format_to(ctx.out(), "loop()");
+            case Op::FREE:
+                return fmt::format_to(ctx.out(), "free()");
+            default:
+                return fmt::format_to(ctx.out(), "[{}]({})",
+                                      (unsigned int) op.code,
+                                      (unsigned int) op.value);
+        }
+    }
+};
+
+struct Factor {
+    unsigned int mode;
+    std::vector<unsigned int> orders;
+};
+
+struct codegen {
+    std::vector<Op> ops;
+
+    void link(unsigned int order) {
+        ops.emplace_back(Op::LINK, order);
+    }
+
+    void push() {
+        ops.emplace_back(Op::PUSH);
+    }
+
+    void pop() {
+        ops.emplace_back(Op::POP);
+    }
+
+    void loop() {
+        ops.emplace_back(Op::LOOP);
+    }
+
+    void free() {
+        ops.emplace_back(Op::FREE);
+    }
+
+    template<typename It>
+    void insert(It begin, It end) {
+        std::copy(begin, end, std::back_inserter(ops));
+    }
+
+    void insert(const codegen &o) {
+        insert(o.ops.begin(), o.ops.end());
+    }
+};
+
+static const std::string GRAMMAR = R"(
+    root    <- term+
+    term    <- product / op
+    op      <- block / link
+    block   <- '(' root ')' / '{' root '}' / '[' root ']'
+    link    <- int / '-'
+    product <- op '*' factor
+    factor  <- int / '{' int+ '}' / '[' int+ ']'
+    int     <- < [0-9]+ >
+
+    %whitespace <- [ \t\n\r]*
+)";
+
+peg::parser build_parser() {
+    peg::parser parser;
+    parser.set_logger([](size_t line, size_t col, const std::string &msg, const std::string &rule) {
+        fmt::print(stderr, "{}:{} [{}] {}\n", line, col, rule, msg);
+    });
+    auto ok = parser.load_grammar(GRAMMAR);
+    assert(ok);
+
+    parser["int"] = [](const peg::SemanticValues &vs) -> std::any {
+        return vs.token_to_number<unsigned int>();
+    };
+
+    parser["link"] = [](const peg::SemanticValues &vs) -> std::any {
+        codegen cg;
+
+        if (vs.choice() == 0) {
+            auto order = std::any_cast<unsigned int>(vs[0]);
+            cg.link(order);
+        } else {
+            cg.free();
+        }
+
+        return cg;
+    };
+
+    parser["root"] = [](const peg::SemanticValues &vs) -> std::any {
+        codegen cg;
+
+        for (const auto &sub: vs) {
+            cg.insert(std::any_cast<codegen>(sub));
+        }
+
+        return cg;
+    };
+
+    parser["block"] = [](const peg::SemanticValues &vs) -> std::any {
+        if (vs.choice() == 0) return vs[0];
+
+        codegen cg;
+
+        cg.push();
+        cg.insert(std::any_cast<codegen>(vs[0]));
+        if (vs.choice() == 1) cg.loop();
+        cg.pop();
+
+        return cg;
+    };
+
+    parser["factor"] = [](const peg::SemanticValues &vs) -> std::any {
+        return Factor{
+            (unsigned int) vs.choice(),
+            vs.transform<unsigned int>(),
+        };
+    };
+
+    parser["product"] = [](const peg::SemanticValues &vs) -> std::any {
+        auto sub = std::any_cast<codegen>(vs[0]);
+        auto fac = std::any_cast<Factor>(vs[1]);
+
+        codegen cg;
+
+        for (const auto &order: fac.orders) {
+            if (fac.mode == 0) {
+                for (unsigned int i = 0; i < order; ++i) {
+                    cg.insert(sub);
+                }
+            } else {
+                cg.push();
+                for (unsigned int i = 0; i < order; ++i) {
+                    cg.insert(sub);
+                }
+                if (fac.mode == 1) cg.loop();
+                cg.pop();
+            }
+        }
+
+        return cg;
+    };
+
+    return parser;
+}
+
+#ifndef NDEBUG
+
+peg::parser build_ast_parser() {
+    peg::parser parser;
+    parser.set_logger([](size_t line, size_t col, const std::string &msg, const std::string &rule) {
+        fmt::print(stderr, "{}:{} [{}] {}\n", line, col, rule, msg);
+    });
+    auto ok = parser.load_grammar(GRAMMAR);
+    assert(ok);
+
+    parser.enable_ast();
+
+    return parser;
+}
+
+#endif
+
+std::vector<Op> compile(const std::string &source) {
+#ifndef NDEBUG
+    static peg::parser ast_parser = build_ast_parser();
+    std::shared_ptr<peg::Ast> ast;
+    bool ast_ok = ast_parser.parse(source, ast);
+    assert(ast_ok);
+//    std::cout << peg::ast_to_s(ast) << std::endl;
+#endif
+
+    static peg::parser parser = build_parser();
+    codegen cg;
+    bool ok = parser.parse(source, cg);
+    assert(ok);
+    return cg.ops;
+}
+
+Graph eval(const std::vector<Op> &ops) {
+    std::vector<std::stack<size_t>> stacks(1);
+
+    Graph g;
+    stacks.back().emplace(g.rank++);
+
+    for (const auto &op: ops) {
+        switch (op.code) {
+            case Op::FREE:
+            case Op::LINK: {
+                tc::Mult order = tc::FREE;
+
+                if (op.code == Op::LINK) {
+                    order = op.value;
+                }
+
+                auto top = stacks.back().top();
+                auto curr = g.rank++;
+
+                stacks.back().emplace(curr);
+                g.edges.emplace_back(top, curr, order);
+
+                break;
+            }
+            case Op::PUSH: {
+                stacks.emplace_back();
+
+                auto ptop = stacks[stacks.size() - 2].top();
+                stacks.back().emplace(ptop);
+
+                break;
+            }
+            case Op::POP: {
+                stacks.pop_back();
+
+                break;
+            }
+            case Op::LOOP: {
+                g.rank--;
+
+                auto ptop = stacks[stacks.size() - 2].top();
+                auto &[top, _, order] = g.edges.back();
+
+                stacks.back().pop();
+                g.edges.pop_back();
+
+                g.edges.emplace_back(top, ptop, order);
+
+                break;
+            }
+            default:
+                throw std::runtime_error("Invalid opcode");
+        }
+    }
+
+    return g;
+}
+
+namespace tc {
+    Group<> coxeter(const std::string &symbol) {
+        auto ops = compile(symbol);
+        auto diagram = eval(ops);
+        Group<> res(diagram.rank);
+        for (const auto &[i, j, m]: diagram.edges) {
+            res.set(i, j, m);
+        }
+        return res;
+    }
+}

tc/src/solve.cpp

@@ -0,0 +1,211 @@
+#include <algorithm>
+#include <queue>
+#include <utility>
+#include <vector>
+
+#include <tc/core.hpp>
+
+namespace tc {
+    /**
+     * Each coset is associated a row in each table.
+     * Rows document the "loops" formed by 
+     */
+    struct Row {
+        bool free: 1;
+        bool idem: 1;
+        unsigned int gnr: 14;  // progress through the loop
+        unsigned int lst_idx: 32;  // the coset that would complete the loop
+
+        Row() : free(true), idem(false), gnr(0), lst_idx(0) {}
+    };
+
+    struct Tables {
+        std::vector<Group<>::Rel> rels;
+        std::vector<std::vector<Row>> rows;
+
+        explicit Tables(std::vector<Group<>::Rel> rels)
+            : rels(std::move(rels)), rows() {
+        }
+
+        [[nodiscard]] size_t size() const {
+            return rels.size();
+        }
+
+        void add_row() {
+            rows.emplace_back(rels.size());
+        }
+    };
+
+    [[nodiscard]] Cosets<> Group<>::solve(std::vector<size_t> const &idxs, size_t bound) const {
+        // region Initialize Cosets Table
+        Cosets<> cosets(rank());
+        cosets.add_row();
+
+        if (rank() == 0) {
+            cosets._complete = true;
+            return cosets;
+        }
+
+        for (size_t g: idxs) {
+            if (g < rank())
+                cosets.set(0, g, 0);
+        }
+        // endregion
+
+        // region Initialize Relation Tables
+        std::vector<Group<>::Rel> rels;
+        for (int i = 0; i < rank(); ++i) {
+            for (int j = i + 1; j < rank(); ++j) {
+                // The algorithm only works for Coxeter groups; multiplicities m_ii=1 are assumed. Relation tables
+                // _may_ be added for them, but they are redundant and hurt performance so are skipped.
+                if (i == j) continue;
+
+                // Coxeter groups admit infinite multiplicities, represented by contexpr tc::FREE. Relation tables
+                // for these should be skipped.
+                auto m = get(i, j);
+
+                if (m == FREE) {
+                    continue;
+                }
+
+                rels.emplace_back(i, j, m);
+            }
+        }
+
+        Tables rel_tables(rels);
+        std::vector<std::vector<size_t>> tables_for(rank());
+        int rel_idx = 0;
+        for (const auto &[i, j, m]: rels) {
+            tables_for[i].push_back(rel_idx);
+            tables_for[j].push_back(rel_idx);
+            rel_idx++;
+        }
+
+        std::vector<size_t> lst_vals;
+        rel_tables.add_row();
+        for (int table_idx = 0; table_idx < rel_tables.size(); ++table_idx) {
+            const auto &[i, j, m] = rel_tables.rels[table_idx];
+            Row &row = rel_tables.rows[0][table_idx];
+
+            if (!cosets.isset(0, i) && !cosets.isset(0, j)) {
+                row.lst_idx = lst_vals.size();
+                lst_vals.push_back(0);
+                row.free = false;
+                row.gnr = 0;
+            } else {
+                row.free = false;
+                row.gnr = 1;
+                row.idem = true;
+            }
+        }
+        // endregion
+
+        size_t idx = 0;
+        size_t fact_idx;
+        size_t coset, gen, target, lst;
+
+        while (true) {
+            // find next unknown product
+            while (idx < cosets.size() and cosets.isset(idx))
+                idx++;
+
+            if (cosets.order() >= bound) {
+                return cosets;
+            }
+
+            // if there are none, then return
+            if (idx == cosets.size()) {
+                // todo unrolled linked list interval
+//                rel_tables.del_rows_to(idx / ngens);  
+                break;
+            }
+
+            // the unknown product must be a new coset, so add it
+            target = cosets.order();
+            cosets.add_row();
+            rel_tables.add_row();
+
+            // queue of products that equal target
+            std::queue<size_t> facts;
+            facts.push(idx);  // new product should be recorded and propagated
+
+            // todo unrolled linked list interval
+//            rel_tables.del_rows_to(coset);
+
+            // find all products which also lead to target
+            while (!facts.empty()) {
+                fact_idx = facts.front();
+                facts.pop();
+
+                // skip if this product was already learned
+                if (cosets.get(fact_idx) != -1) continue;
+
+                cosets.set(fact_idx, target);
+
+                coset = fact_idx / rank();
+                gen = fact_idx % rank();
+
+                // If the product stays within the coset todo
+                for (size_t table_idx: tables_for[gen]) {
+                    auto &[i, j, m] = rel_tables.rels[table_idx];
+                    auto &trow = rel_tables.rows[target][table_idx];
+                    auto &crow = rel_tables.rows[coset][table_idx];
+
+                    size_t other_gen = (i == gen) ? j : i;
+
+                    // Test if loop is closed
+                    if (trow.free) {
+                        trow = crow;
+                        trow.gnr++;
+
+                        if (target == coset) {
+                            trow.idem = true;
+                        }
+
+                        if (trow.idem) {
+                            if (trow.gnr == m) {
+                                // loop is closed, but idempotent, so the target links to itself via the other generator.
+                                // todo might be able to move this logic up into the (target == coset) block and avoid those computations.
+                                facts.push(target * rank() + other_gen);
+                            }
+                        } else {
+                            if (trow.gnr == m - 1) {
+                                // loop is almost closed. record that the target closes this loop.
+                                lst_vals[trow.lst_idx] = target;
+                            } else if (trow.gnr == m) {
+                                // loop is closed. We know the last element in the loop must link with this one. 
+                                lst = lst_vals[trow.lst_idx];
+//                            delete trow.lst_ptr;
+                                facts.push(lst * rank() + other_gen);
+                            }
+                        }
+                    }
+                }
+            }
+
+            // If any target row wasn't identified with a loop,
+            // then assign it a new loop.
+            for (size_t table_idx = 0; table_idx < rel_tables.size(); table_idx++) {
+                auto &[i, j, m] = rel_tables.rels[table_idx];
+                auto &trow = rel_tables.rows[target][table_idx];
+
+                if (trow.free) {
+                    if ((cosets.get(target, i) != target) and
+                        (cosets.get(target, j) != target)) {
+                        trow.lst_idx = lst_vals.size();
+                        trow.free = false;
+                        lst_vals.push_back(0);
+                        trow.gnr = 0;
+                    } else {
+                        trow.free = false;
+                        trow.gnr = 1;
+                        trow.idem = true;
+                    }
+                }
+            }
+        }
+
+        cosets._complete = true;
+        return cosets;
+    }
+}

tc/test/CMakeLists.txt

@@ -0,0 +1,18 @@
+include(GoogleTest)
+
+add_executable(test_solve test_solve.cpp)
+target_link_libraries(test_solve tc GTest::gtest_main)
+
+add_executable(test_lang test_lang.cpp)
+target_link_libraries(test_lang PUBLIC tc::tc GTest::gtest_main)
+
+set(MIN_DEBUG_CPS 200000)
+set(MIN_RELEASE_CPS 1000000)
+
+target_compile_definitions(
+    test_solve PUBLIC
+    MIN_COS_PER_SEC=$<IF:$<CONFIG:Debug>,${MIN_DEBUG_CPS},${MIN_RELEASE_CPS}>
+)
+
+gtest_discover_tests(test_solve)
+gtest_discover_tests(test_lang)

tc/test/test_lang.cpp

@@ -0,0 +1,54 @@
+#include <tc/core.hpp>
+#include <tc/groups.hpp>
+
+#include <gtest/gtest.h>
+
+// "5 3 3"
+// "5 (3 3)"
+// "[5 3 3]"
+// "[4 3 [3 5] 3]"
+// "{3 4 5 6 7 8 9}"
+// "3 {3 3 [4] 3} 5"
+// "5 * 3"
+// "5 * [3]"
+// "5 * {2 3}"
+// "5 * [3 2]"
+// "(5 2) * [3 2]"
+// "4 [3 * [2 3]] 5"
+// "{3 * 3} [4] [5]"
+
+TEST(coxeter, simple) {
+    auto g = tc::coxeter("5 3 3");
+
+    ASSERT_EQ(g.rank(), 4);
+
+    EXPECT_EQ(g.get(0, 0), 1);
+    EXPECT_EQ(g.get(0, 1), 5);
+    EXPECT_EQ(g.get(0, 2), 2);
+    EXPECT_EQ(g.get(0, 3), 2);
+    
+    EXPECT_EQ(g.get(1, 0), 5);
+    EXPECT_EQ(g.get(1, 1), 1);
+    EXPECT_EQ(g.get(1, 2), 3);
+    EXPECT_EQ(g.get(1, 3), 2);
+    
+    EXPECT_EQ(g.get(2, 0), 2);
+    EXPECT_EQ(g.get(2, 1), 3);
+    EXPECT_EQ(g.get(2, 2), 1);
+    EXPECT_EQ(g.get(2, 3), 3);
+    
+    EXPECT_EQ(g.get(3, 0), 2);
+    EXPECT_EQ(g.get(3, 1), 2);
+    EXPECT_EQ(g.get(3, 2), 3);
+    EXPECT_EQ(g.get(3, 3), 1);
+}
+
+TEST(coxeter, looping) {
+    auto g = tc::coxeter("{5 3 4}");
+
+    ASSERT_EQ(g.rank(), 3);
+
+    EXPECT_EQ(g.get(0, 1), 5);
+    EXPECT_EQ(g.get(1, 2), 3);
+    EXPECT_EQ(g.get(2, 0), 4);
+}

tc/test/test_solve.cpp

@@ -0,0 +1,225 @@
+#include <ctime>
+#include <vector>
+
+#include <tc/groups.hpp>
+#include <tc/core.hpp>
+
+#include <gtest/gtest.h>
+
+/// helper for testing solve and testing speed
+testing::AssertionResult AssertSolveOrder(
+    const char *group_expr,
+    const char *sub_gens_expr,
+    const char *expected_order_expr,
+    const tc::Group<> &group,
+    const std::vector<size_t> &sub_gens,
+    size_t expected_order
+) {
+    auto start = std::clock();
+    auto cosets = group.solve(sub_gens);
+    auto end = std::clock();
+
+    size_t actual_order = cosets.order();
+
+    auto total_sec = (double) (end - start) / CLOCKS_PER_SEC;
+    auto cosets_per_sec = (double) actual_order / total_sec;
+
+    bool order_good = actual_order == expected_order;
+    bool speed_good = cosets_per_sec >= MIN_COS_PER_SEC || total_sec < 0.0001;
+    // extremely short times cause false negatives. ex. A2 can be solved in only 3 clocks.
+
+    if (order_good && speed_good) {
+        return testing::AssertionSuccess();
+    }
+
+    testing::AssertionResult res = testing::AssertionFailure();
+    res << group_expr << " / " << sub_gens_expr << " :";
+    if (!order_good) {
+        res << " Gave order " << actual_order << " but expected order " << expected_order << ".";
+    }
+    if (!speed_good) {
+        res << " Solution too slow (" << cosets_per_sec
+            << " cos/s < " << MIN_COS_PER_SEC << ")."
+            << " " << std::fixed << total_sec << " s.";
+    }
+
+    return res;
+}
+
+#define EXPECT_SOLVE_ORDER(group, sub_gens, expected_order) EXPECT_PRED_FORMAT3(AssertSolveOrder, group, sub_gens, expected_order);
+
+using v = std::vector<size_t>;
+
+tc::Group<> A(unsigned int n) {
+    return tc::vcoxeter("3 * {}", {n - 1});
+}
+
+tc::Group<> B(unsigned int n) {
+    return tc::vcoxeter("4 3 * {}", {n - 2});
+}
+
+tc::Group<> D(unsigned int n) {
+    return tc::vcoxeter("3 * [1 1 {}]", {n - 3});
+}
+
+tc::Group<> E(unsigned int n) {
+    return tc::vcoxeter("3 * [1 2 {}]", {n - 4});
+}
+
+tc::Group<> F4() {
+    return tc::coxeter("3 4 3");
+}
+
+tc::Group<> G2() {
+    return tc::coxeter("6");
+}
+
+tc::Group<> H(unsigned int n) {
+    return tc::vcoxeter("5 3 * {}", {n - 2});
+}
+
+tc::Group<> I2(unsigned int n) {
+    return tc::vcoxeter("{}", {n});
+}
+
+tc::Group<> T(unsigned int m, unsigned int n) {
+    return tc::vcoxeter("{} 2 {}", {m, n});
+}
+
+tc::Group<> T(unsigned int n) {
+    return T(n, n);
+}
+
+// See the group orders here https://en.wikipedia.org/wiki/Coxeter_group#Properties
+
+TEST(solve, A) {
+    EXPECT_SOLVE_ORDER(A(1), v({}), 2);
+    EXPECT_SOLVE_ORDER(A(2), v({}), 6);
+    EXPECT_SOLVE_ORDER(A(3), v({}), 24);
+    EXPECT_SOLVE_ORDER(A(3), v({0}), 12);
+    EXPECT_SOLVE_ORDER(A(3), v({0, 1}), 4);
+    EXPECT_SOLVE_ORDER(A(3), v({0, 2}), 6);
+    EXPECT_SOLVE_ORDER(A(3), v({2}), 12);
+    EXPECT_SOLVE_ORDER(A(4), v({}), 120);
+    EXPECT_SOLVE_ORDER(A(4), v({0}), 60);
+    EXPECT_SOLVE_ORDER(A(4), v({0, 1}), 20);
+    EXPECT_SOLVE_ORDER(A(4), v({2}), 60);
+    EXPECT_SOLVE_ORDER(A(4), v({0, 2}), 30);
+}
+
+TEST(solve, B) {
+    EXPECT_SOLVE_ORDER(B(2), v({}), 8);
+    EXPECT_SOLVE_ORDER(B(3), v({}), 48);
+    EXPECT_SOLVE_ORDER(B(3), v({0}), 24);
+    EXPECT_SOLVE_ORDER(B(3), v({0, 2}), 12);
+    EXPECT_SOLVE_ORDER(B(4), v({}), 384);
+    EXPECT_SOLVE_ORDER(B(4), v({0}), 192);
+    EXPECT_SOLVE_ORDER(B(4), v({0, 2}), 96);
+    EXPECT_SOLVE_ORDER(B(5), v({}), 3840);
+    EXPECT_SOLVE_ORDER(B(5), v({0}), 1920);
+    EXPECT_SOLVE_ORDER(B(5), v({0, 2}), 960);
+    EXPECT_SOLVE_ORDER(B(5), v({0, 2, 3}), 320);
+    EXPECT_SOLVE_ORDER(B(6), v({}), 46080);
+    EXPECT_SOLVE_ORDER(B(6), v({0}), 23040);
+    EXPECT_SOLVE_ORDER(B(6), v({0, 2}), 11520);
+    EXPECT_SOLVE_ORDER(B(6), v({0, 2, 3}), 3840);
+    EXPECT_SOLVE_ORDER(B(6), v({0, 2, 3, 5}), 1920);
+}
+
+TEST(solve, D) {
+    EXPECT_SOLVE_ORDER(D(3), v({}), 24);
+    EXPECT_SOLVE_ORDER(D(4), v({}), 192);
+    EXPECT_SOLVE_ORDER(D(4), v({0, 1}), 32);
+    EXPECT_SOLVE_ORDER(D(4), v({0, 1, 3}), 8);
+    EXPECT_SOLVE_ORDER(D(5), v({}), 1920);
+    EXPECT_SOLVE_ORDER(D(5), v({0, 1}), 320);
+    EXPECT_SOLVE_ORDER(D(5), v({0, 1, 3}), 80);
+    EXPECT_SOLVE_ORDER(D(5), v({0, 1, 3, 4}), 16);
+    EXPECT_SOLVE_ORDER(D(6), v({}), 23040);
+    EXPECT_SOLVE_ORDER(D(6), v({0, 1}), 3840);
+    EXPECT_SOLVE_ORDER(D(6), v({0, 1, 3}), 960);
+    EXPECT_SOLVE_ORDER(D(6), v({0, 1, 3, 5}), 480);
+}
+
+TEST(solve, E) {
+    EXPECT_SOLVE_ORDER(E(4), v({}), 120);
+    EXPECT_SOLVE_ORDER(E(4), v({2}), 60);
+    EXPECT_SOLVE_ORDER(E(4), v({2, 1}), 30);
+    EXPECT_SOLVE_ORDER(E(4), v({2, 1, 3}), 10);
+    EXPECT_SOLVE_ORDER(E(5), v({}), 1920);
+    EXPECT_SOLVE_ORDER(E(5), v({2}), 960);
+    EXPECT_SOLVE_ORDER(E(5), v({2, 1}), 480);
+    EXPECT_SOLVE_ORDER(E(5), v({2, 1, 3}), 160);
+    EXPECT_SOLVE_ORDER(E(6), v({}), 51840);
+    EXPECT_SOLVE_ORDER(E(6), v({2}), 25920);
+    EXPECT_SOLVE_ORDER(E(6), v({2, 1}), 12960);
+    EXPECT_SOLVE_ORDER(E(6), v({2, 1, 3}), 4320);
+}
+
+TEST(solve, F) {
+    EXPECT_SOLVE_ORDER(F4(), v({}), 1152);
+    EXPECT_SOLVE_ORDER(F4(), v({0}), 576);
+    EXPECT_SOLVE_ORDER(F4(), v({0, 2}), 288);
+    EXPECT_SOLVE_ORDER(F4(), v({1, 3}), 288);
+    EXPECT_SOLVE_ORDER(F4(), v({1, 2, 3}), 24);
+}
+
+TEST(solve, G) {
+    EXPECT_SOLVE_ORDER(G2(), v({}), 12);
+    EXPECT_SOLVE_ORDER(G2(), v({0}), 6);
+    EXPECT_SOLVE_ORDER(G2(), v({1}), 6);
+}
+
+TEST(solve, H) {
+    EXPECT_SOLVE_ORDER(H(2), v({}), 10);
+    EXPECT_SOLVE_ORDER(H(2), v({0}), 5);
+    EXPECT_SOLVE_ORDER(H(2), v({1}), 5);
+    EXPECT_SOLVE_ORDER(H(3), v({}), 120);
+    EXPECT_SOLVE_ORDER(H(3), v({0}), 60);
+    EXPECT_SOLVE_ORDER(H(3), v({0, 1}), 12);
+    EXPECT_SOLVE_ORDER(H(3), v({0, 2}), 30);
+    EXPECT_SOLVE_ORDER(H(3), v({1, 2}), 20);
+    EXPECT_SOLVE_ORDER(H(4), v({}), 14400);
+    EXPECT_SOLVE_ORDER(H(4), v({0}), 7200);
+    EXPECT_SOLVE_ORDER(H(4), v({1}), 7200);
+    EXPECT_SOLVE_ORDER(H(4), v({1, 3}), 3600);
+}
+
+TEST(solve, I) {
+    EXPECT_SOLVE_ORDER(I2(2), v({}), 4);
+    EXPECT_SOLVE_ORDER(I2(3), v({}), 6);
+    EXPECT_SOLVE_ORDER(I2(3), v({0}), 3);
+    EXPECT_SOLVE_ORDER(I2(3), v({1}), 3);
+    EXPECT_SOLVE_ORDER(I2(4), v({}), 8);
+    EXPECT_SOLVE_ORDER(I2(4), v({0}), 4);
+    EXPECT_SOLVE_ORDER(I2(4), v({1}), 4);
+    EXPECT_SOLVE_ORDER(I2(5), v({}), 10);
+    EXPECT_SOLVE_ORDER(I2(5), v({0}), 5);
+    EXPECT_SOLVE_ORDER(I2(5), v({1}), 5);
+}
+
+TEST(solve, T) {
+    EXPECT_SOLVE_ORDER(T(3), v({}), 36);
+    EXPECT_SOLVE_ORDER(T(4), v({}), 64);
+    EXPECT_SOLVE_ORDER(T(400), v({}), 640000);
+    EXPECT_SOLVE_ORDER(T(400), v({0}), 320000);
+    EXPECT_SOLVE_ORDER(T(400), v({0, 2}), 160000);
+    EXPECT_SOLVE_ORDER(T(400, 300), v({}), 480000);
+    EXPECT_SOLVE_ORDER(T(400, 300), v({0}), 240000);
+    EXPECT_SOLVE_ORDER(T(400, 300), v({0, 2}), 120000);
+}
+
+TEST(solve_large, B) {
+    EXPECT_SOLVE_ORDER(B(7), v({}), 645120);
+    EXPECT_SOLVE_ORDER(B(8), v({}), 10321920);
+}
+
+TEST(solve_large, E) {
+    EXPECT_SOLVE_ORDER(E(6), v({}), 51840);
+    EXPECT_SOLVE_ORDER(E(7), v({}), 2903040);
+}
+
+TEST(solve_large, T) {
+    EXPECT_SOLVE_ORDER(T(500), v({}), 1000000);
+    EXPECT_SOLVE_ORDER(T(1000), v({}), 4000000);
+}

vendor/fmt.cmake

@@ -0,0 +1,7 @@
+FetchContent_Declare(
+    fmt
+    GIT_REPOSITORY https://github.com/fmtlib/fmt.git
+    GIT_TAG 9.1.0
+    GIT_PROGRESS TRUE
+)
+FetchContent_MakeAvailable(fmt)

vendor/gtest.cmake

@@ -0,0 +1,8 @@
+FetchContent_Declare(
+    gtest
+    GIT_REPOSITORY https://github.com/google/googletest.git
+    GIT_TAG release-1.12.1
+    GIT_PROGRESS TRUE
+)
+set(gtest_force_shared_crt ON CACHE BOOL "" FORCE)
+FetchContent_MakeAvailable(gtest)

vendor/peglib.cmake

@@ -0,0 +1,16 @@
+FetchContent_Declare(
+    peglib
+    GIT_REPOSITORY https://github.com/yhirose/cpp-peglib
+    GIT_TAG v1.8.2
+    GIT_PROGRESS TRUE
+)
+set(PEGLIB_BUILD_TESTS OFF CACHE INTERNAL "")
+FetchContent_GetProperties(peglib)
+if(NOT ${peglib}_POPULATED)
+    FetchContent_Populate(peglib)
+    find_package(Threads)
+    add_library(peglib INTERFACE ${peglib_SOURCE_DIR}/peglib.h)
+    target_include_directories(peglib INTERFACE ${peglib_SOURCE_DIR})
+    target_link_libraries(peglib INTERFACE Threads::Threads)
+    add_library(peglib::peglib ALIAS peglib)
+endif()

vis/include/geometry.hpp

@@ -32,8 +32,8 @@ struct Primitive {
         if (N > 1) std::swap(inds[0], inds[1]);
     }
 
-    void apply(const tc::Cosets &table, int gen) {
+    void apply(const tc::Cosets<> &table, int gen) {
-        for (auto &ind : inds) {
+        for (auto &ind: inds) {
             ind = table.get(ind, gen);
         }
         flip();
@@ -43,8 +43,8 @@ struct Primitive {
 /**
  * Produce a list of all generators for the group context. The range [0..group.ngens).
  */
-std::vector<int> generators(const tc::Group &context) {
+std::vector<size_t> generators(const tc::Group<> &context) {
-    std::vector<int> g_gens(context.ngens);
+    std::vector<size_t> g_gens(context.rank());
     std::iota(g_gens.begin(), g_gens.end(), 0);
     return g_gens;
 }
@@ -52,19 +52,19 @@ std::vector<int> generators(const tc::Group &context) {
 /**
  * Determine which of g_gens are the correct names for sg_gens within the current context
  */
-std::vector<int> recontext_gens(
+std::vector<size_t> recontext_gens(
-    const tc::Group &context,
+    const tc::Group<> &context,
-    std::vector<int> g_gens,
+    std::vector<size_t> g_gens,
-    std::vector<int> sg_gens) {
+    std::vector<size_t> sg_gens) {
 
     std::sort(g_gens.begin(), g_gens.end());
 
-    int inv_gen_map[context.ngens];
+    int inv_gen_map[context.rank()];
     for (size_t i = 0; i < g_gens.size(); i++) {
         inv_gen_map[g_gens[i]] = i;
     }
 
-    std::vector<int> s_sg_gens;
+    std::vector<size_t> s_sg_gens;
     s_sg_gens.reserve(sg_gens.size());
     for (const auto gen : sg_gens) {
         s_sg_gens.push_back(inv_gen_map[gen]);
@@ -78,9 +78,9 @@ std::vector<int> recontext_gens(
  * Determine whether the orientation of the group sg_gens is reversed from the group g_gens within group context
  */
 int get_parity(
-    const tc::Group &context,
+    const tc::Group<> &context,
-    const std::vector<int> &g_gens,
+    const std::vector<size_t> &g_gens,
-    const std::vector<int> &sg_gens
+    const std::vector<size_t> &sg_gens
 ) {
     if (g_gens.size() != sg_gens.size() + 1) return 0;
 
@@ -99,21 +99,21 @@ int get_parity(
 /**
  * Solve the cosets generated by sg_gens within the subgroup generated by g_gens of the group context
  */
-tc::Cosets solve(
+tc::Cosets<> solve(
-    const tc::Group &context,
+    const tc::Group<> &context,
-    const std::vector<int> &g_gens,
+    const std::vector<size_t> &g_gens,
-    const std::vector<int> &sg_gens
+    const std::vector<size_t> &sg_gens
 ) {
     const auto proper_sg_gens = recontext_gens(context, g_gens, sg_gens);
-    return context.subgroup(g_gens).solve(proper_sg_gens);
+    return context.sub(g_gens).solve(proper_sg_gens);
 }
 
 /**
  * Apply some context transformation to all primitives of this mesh.
  */
 template<unsigned N>
-std::vector<Primitive<N>> apply(std::vector<Primitive<N>> prims, const tc::Cosets &table, int gen) {
+std::vector<Primitive<N>> apply(std::vector<Primitive<N>> prims, const tc::Cosets<> &table, int gen) {
-    for (auto &prim : prims) {
+    for (auto &prim: prims) {
         prim.apply(table, gen);
     }
     return prims;
@@ -136,17 +136,20 @@ template<unsigned N>
 [[nodiscard]]
 std::vector<Primitive<N>> recontext(
     std::vector<Primitive<N>> prims,
-    const tc::Group &context,
+    const tc::Group<> &context,
-    const std::vector<int> &g_gens,
+    const std::vector<size_t> &g_gens,
-    const std::vector<int> &sg_gens
+    const std::vector<size_t> &sg_gens
 ) {
     const auto proper_sg_gens = recontext_gens(context, g_gens, sg_gens);
     const auto table = solve(context, g_gens, {});
-    const auto path = solve(context, sg_gens, {}).path;
+    const auto &cosets = solve(context, sg_gens, {});
 
-    auto map = path.template walk<int, int>(0, proper_sg_gens, [table](int coset, int gen) {
+    tc::Path<size_t> path(cosets, proper_sg_gens);
+
+    std::vector<size_t> map(path.order());
+    path.walk(0, [&table](size_t coset, size_t gen){
         return table.get(coset, gen);
-    });
+    }, map.begin());
 
     std::vector<Primitive<N>> res(prims);
     for (Primitive<N> &prim : res) {
@@ -184,21 +187,22 @@ template<unsigned N>
 [[nodiscard]]
 std::vector<std::vector<Primitive<N>>> each_tile(
     std::vector<Primitive<N>> prims,
-    const tc::Group &context,
+    const tc::Group<> &context,
-    const std::vector<int> &g_gens,
+    const std::vector<size_t> &g_gens,
-    const std::vector<int> &sg_gens
+    const std::vector<size_t> &sg_gens
 ) {
     std::vector<Primitive<N>> base = recontext(prims, context, g_gens, sg_gens);
     const auto proper_sg_gens = recontext_gens(context, g_gens, sg_gens);
 
     const auto table = solve(context, g_gens, {});
-    const auto path = solve(context, g_gens, sg_gens).path;
+    tc::Path path(solve(context, g_gens, sg_gens));
 
     auto _gens = generators(context);
 
-    auto res = path.walk<std::vector<Primitive<N>>, int>(base, generators(context), [&](auto from, auto gen) {
+    std::vector<std::vector<Primitive<N>>> res(path.order());
-        return apply(from, table, gen);
+    path.walk(base, [&](auto from, auto to) {
-    });
+        return apply(from, table, to);
+    }, res.begin());
 
     return res;
 }
@@ -207,9 +211,9 @@ template<unsigned N>
 [[nodiscard]]
 std::vector<Primitive<N>> tile(
     std::vector<Primitive<N>> prims,
-    const tc::Group &context,
+    const tc::Group<> &context,
-    const std::vector<int> &g_gens,
+    const std::vector<size_t> &g_gens,
-    const std::vector<int> &sg_gens
+    const std::vector<size_t> &sg_gens
 ) {
     auto res = each_tile<N>(prims, context, g_gens, sg_gens);
 
@@ -221,12 +225,12 @@ std::vector<Primitive<N>> tile(
  */
 template<unsigned N>
 [[nodiscard]]
-std::vector<Primitive<N + 1>> fan(std::vector<Primitive<N>> prims, int root) {
+std::vector<Primitive<N + 1>> fan(std::vector<Primitive<N>> prims, size_t root) {
     std::vector<Primitive<N + 1>> res(prims.size());
     std::transform(prims.begin(), prims.end(), res.begin(),
-        [root](const Primitive<N> &prim) {
+                   [root](const Primitive<N> &prim) {
-            return Primitive<N + 1>(prim, root);
+                       return Primitive<N + 1>(prim, root);
-        }
+                   }
     );
     return res;
 }
@@ -236,8 +240,8 @@ std::vector<Primitive<N + 1>> fan(std::vector<Primitive<N>> prims, int root) {
  */
 template<unsigned N>
 std::vector<Primitive<N>> triangulate(
-    const tc::Group &context,
+    const tc::Group<> &context,
-    const std::vector<int> &g_gens
+    const std::vector<size_t> &g_gens
 ) {
     if (g_gens.size() + 1 != N) // todo make static assert
         throw std::logic_error("g_gens size must be one less than N");
@@ -261,8 +265,8 @@ std::vector<Primitive<N>> triangulate(
  */
 template<>
 std::vector<Primitive<1>> triangulate(
-    const tc::Group &context,
+    const tc::Group<> &context,
-    const std::vector<int> &g_gens
+    const std::vector<size_t> &g_gens
 ) {
     if (not g_gens.empty()) // todo make static assert
         throw std::logic_error("g_gens must be empty for a trivial Mesh");
@@ -273,12 +277,12 @@ std::vector<Primitive<1>> triangulate(
 }
 
 template<unsigned N, class T>
-auto hull(const tc::Group &group, T all_sg_gens, const std::vector<std::vector<int>> &exclude) {
+auto hull(const tc::Group<> &group, T all_sg_gens, const std::vector<std::vector<size_t>> &exclude) {
     std::vector<std::vector<Primitive<N>>> parts;
     auto g_gens = generators(group);
-    for (const std::vector<int> &sg_gens : all_sg_gens) {
+    for (std::vector<size_t> sg_gens: all_sg_gens) {
         bool excluded = false;
-        for (const auto &test : exclude) {
+        for (const auto &test: exclude) {
             if (sg_gens == test) {
                 excluded = true;
                 break;

vis/include/mirror.hpp

@@ -94,10 +94,10 @@ float dot(const V &a, const V &b) {
 }
 
 template<unsigned N>
-std::vector<vec<N>> mirror(const tc::Group &group) {
+std::vector<vec<N>> mirror(const tc::Group<> &group) {
     std::vector<std::vector<float>> mirrors;
 
-    for (int p = 0; p < group.ngens; ++p) {
+    for (int p = 0; p < group.rank(); ++p) {
         std::vector<float> vp;
         for (int m = 0; m < p; ++m) {
             auto &vq = mirrors[m];

vis/src/main.cpp

@@ -78,17 +78,22 @@ Matrices build(GLFWwindow *window, State &state) {
 }
 
 template<class C>
-std::vector<vec4> points(const tc::Group &group, const C &coords) {
+std::vector<vec4> points(const tc::Group<> &group, const C &coords) {
     auto cosets = group.solve();
     auto mirrors = mirror<5>(group);
 
+    tc::Path<vec5> path(cosets, mirrors);
+
     auto corners = plane_intersections(mirrors);
 
     auto start = barycentric(corners, coords);
 
-    const auto &higher = cosets.path.walk<vec5, vec5>(start, mirrors, reflect<vec5>);
+    std::vector<vec5> higher(path.order());
+    path.walk(start, reflect<vec5>, higher.begin());
+
     std::vector<vec4> lower(higher.size());
     std::transform(higher.begin(), higher.end(), lower.begin(), stereo<4>);
+
     return lower;
 }
 
@@ -132,11 +137,11 @@ struct SliceProp : public Prop<N> {
 
     template<class T, class C>
     static SliceProp<N> build(
-        const tc::Group &g,
+        const tc::Group<> &g,
         const C &coords,
         vec3 color,
         T all_sg_gens,
-        const std::vector<std::vector<int>> &exclude
+        const std::vector<std::vector<size_t>> &exclude
     ) {
         SliceProp<N> res(color);
 
@@ -225,13 +230,13 @@ struct WireframeProp : public Prop<2> {
     WireframeProp(WireframeProp &&) noexcept = default;
 
     template<class T, class C>
-    static WireframeProp build(const tc::Group &g,
+    static WireframeProp build(const tc::Group<> &g,
         const C &coords,
         bool curve,
         bool ortho,
         vec3 color,
         T all_sg_gens,
-        const std::vector<std::vector<int>> &exclude
+        const std::vector<std::vector<size_t>> &exclude
     ) {
         WireframeProp res(color);
 
@@ -274,10 +279,10 @@ void run(const std::string &config_file, GLFWwindow *window) {
     State state{};
     glfwSetWindowUserPointer(window, &state);
 
-    state.dimension = scene["dimension"].as<int>();
+    state.dimension = scene["dimension"].as<size_t>();
 
     for (const auto &group_info : scene["groups"]) {
-        auto symbol = group_info["symbol"].as<std::vector<int>>();
+        auto symbol = group_info["symbol"].as<std::vector<unsigned int>>();
         auto group = tc::schlafli(symbol);
         auto gens = generators(group);
 
@@ -285,19 +290,19 @@ void run(const std::string &config_file, GLFWwindow *window) {
             for (const auto &slice_info : group_info["slices"]) {
                 auto root = slice_info["root"].as<vec5>();
                 auto color = slice_info["color"].as<vec3>();
-                auto exclude = std::vector<std::vector<int>>();
+                auto exclude = std::vector<std::vector<size_t>>();
 
                 if (slice_info["exclude"].IsDefined()) {
-                    exclude = slice_info["exclude"].as<std::vector<std::vector<int>>>();
+                    exclude = slice_info["exclude"].as<std::vector<std::vector<size_t>>>();
                 }
 
                 if (slice_info["subgroups"].IsDefined()) {
-                    auto subgroups = slice_info["subgroups"].as<std::vector<std::vector<int>>>();
+                    auto subgroups = slice_info["subgroups"].as<std::vector<std::vector<size_t>>>();
                     sRen.props.push_back(SliceProp<4>::build(
                         group, root, color, subgroups, exclude
                     ));
                 } else {
-                    auto combos = Combos<int>(gens, 3);
+                    auto combos = Combos<size_t>(gens, 3);
                     sRen.props.push_back(SliceProp<4>::build(
                         group, root, color, combos, exclude
                     ));
@@ -309,16 +314,16 @@ void run(const std::string &config_file, GLFWwindow *window) {
             for (const auto &wire_info : group_info["wires"]) {
                 auto root = wire_info["root"].as<vec5>();
                 auto color = wire_info["color"].as<vec3>();
-                auto exclude = std::vector<std::vector<int>>();
+                auto exclude = std::vector<std::vector<size_t>>();
                 auto curve = wire_info["curve"].IsDefined() && wire_info["curve"].as<bool>();
                 auto ortho = wire_info["ortho"].IsDefined() && wire_info["ortho"].as<bool>();
 
                 if (wire_info["exclude"].IsDefined()) {
-                    exclude = wire_info["exclude"].as<std::vector<std::vector<int>>>();
+                    exclude = wire_info["exclude"].as<std::vector<std::vector<size_t>>>();
                 }
 
                 if (wire_info["subgroups"].IsDefined()) {
-                    auto subgroups = wire_info["subgroups"].as<std::vector<std::vector<int>>>();
+                    auto subgroups = wire_info["subgroups"].as<std::vector<std::vector<size_t>>>();
 
                     if (ortho && curve) {
                         wocRen.props.push_back(WireframeProp::build(
@@ -338,7 +343,7 @@ void run(const std::string &config_file, GLFWwindow *window) {
                         ));
                     }
                 } else {
-                    auto combos = Combos<int>(gens, 1);
+                    auto combos = Combos<size_t>(gens, 1);
 
                     if (ortho && curve) {
                         wocRen.props.push_back(WireframeProp::build(