Merge pull request #1 from allemangD/slices

Slices

.gitmodules

@@ -4,3 +4,6 @@
 [submodule "vendor/glm"]
 	path = vendor/glm
 	url = https://github.com/g-truc/glm.git
+[submodule "vendor/toddcox"]
+	path = vendor/toddcox
+	url = https://github.com/JCRaymond/toddcox-faster.git

.idea/misc.xml

@@ -0,0 +1,18 @@
+<?xml version="1.0" encoding="UTF-8"?>
+<project version="4">
+  <component name="CMakeWorkspace" PROJECT_DIR="$PROJECT_DIR$" />
+  <component name="CidrRootsConfiguration">
+    <sourceRoots>
+      <file path="$PROJECT_DIR$/examples" />
+      <file path="$PROJECT_DIR$/vis/include" />
+      <file path="$PROJECT_DIR$/vis/shaders" />
+      <file path="$PROJECT_DIR$/vis/src" />
+    </sourceRoots>
+    <libraryRoots>
+      <file path="$PROJECT_DIR$/vendor" />
+    </libraryRoots>
+  </component>
+  <component name="JavaScriptSettings">
+    <option name="languageLevel" value="ES6" />
+  </component>
+</project>

.idea/toddcox-visualize.iml

@@ -0,0 +1,2 @@
+<?xml version="1.0" encoding="UTF-8"?>
+<module classpath="CMake" type="CPP_MODULE" version="4" />

.idea/vcs.xml

@@ -0,0 +1,6 @@
+<?xml version="1.0" encoding="UTF-8"?>
+<project version="4">
+  <component name="VcsDirectoryMappings">
+    <mapping directory="$PROJECT_DIR$" vcs="Git" />
+  </component>
+</project>

CMakeLists.txt

@@ -3,10 +3,11 @@ project(toddcox-faster)
 
 set(CMAKE_CXX_STANDARD 17)
 
+add_subdirectory(vendor/toddcox)
 add_subdirectory(vendor/glad)
 add_subdirectory(vendor/glfw)
 add_subdirectory(vendor/glm)
 
-add_subdirectory(tc)
 add_subdirectory(vis)
-add_subdirectory(example)
+
+add_subdirectory(examples)

README.md

@@ -0,0 +1,9 @@
+# toddcox-visualize
+A new (basic) implementation of the Todd-Coxeter algorithm for Coxeter groups in C++ that beats GAP.
+
+---
+
+Example output with group F4.
+
+![Example render](screen.png)
+ 

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,31 +0,0 @@
-#include "tc/solver.h"
-#include "tc/groups.h"
-
-#include <ctime>
-#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),
-    };
-
-    for (const auto &group : groups) {
-        auto s = std::clock(); // to measure CPU time
-        auto cosets = tc::solve(group);
-        auto e = std::clock();
-
-        double diff = (double) (e - s) / CLOCKS_PER_SEC;
-        int order = cosets.len;
-
-        std::cout << group.name << "," << order << "," << diff << std::endl;
-    }
-
-    return 0;
-}

example/path.cpp

@@ -1,17 +0,0 @@
-#include "tc/solver.h"
-#include "tc/groups.h"
-
-#include <ctime>
-#include <iostream>
-
-int main() {
-    auto cube = tc::group::B(3);
-    auto vars = tc::solve(cube, {});
-
-    for (int target = 1; target < vars.len; target++) {
-        auto &action = vars.path[target];
-        std::cout << action.coset << " " << action.gen << " " << target << std::endl;
-    }
-
-    return 0;
-}

examples/CMakeLists.txt

@@ -0,0 +1,8 @@
+add_executable(memo memotest.cpp)
+target_link_libraries(memo PRIVATE tc vis-util)
+
+add_executable(geom geomtest.cpp)
+target_link_libraries(geom PRIVATE tc vis-util)
+
+add_executable(sub subtest.cpp)
+target_link_libraries(sub PRIVATE vis-util)

examples/geomtest.cpp

@@ -0,0 +1,46 @@
+//
+// Created by raymo on 1/20/2020.
+//
+
+#include <tc/groups.hpp>
+#include <geometry.hpp>
+#include <iostream>
+
+int main() {
+    auto g = tc::schlafli({3, 2});
+    GeomGen gg(g);
+
+    auto path = gg.solve().path;
+
+    //std::vector<std::string>  = {"a", "b", "c"};
+    std::string base = "";
+    auto words = path.walk<std::string, std::string>(base, {"a", "b", "c"}, [](auto s1, auto g) { return s1 + g; });
+    for (const auto word : words) {
+        std::cout << word << std::endl;
+    }
+
+    std::vector<int> gens = {0, 1, 2};
+    auto s = gg.triangulate(gens);
+    s.print();
+    return 0;
+
+    auto g_gens = gg.group_gens();
+    std::vector<int> sg_gens = {1, 2};
+    auto ns = gg.tile(g_gens, sg_gens, s);
+
+    std::cout << "Before: " << std::endl;
+    std::cout << '\t';
+    for (int val : s.vals) {
+        std::cout << val << " ";
+    }
+    std::cout << std::endl;
+
+    std::cout << " After: " << std::endl;
+    std::cout << '\t';
+    for (int val : ns.vals) {
+        std::cout << val << " ";
+    }
+    std::cout << std::endl;
+
+    return 0;
+}

examples/memotest.cpp

@@ -0,0 +1,70 @@
+#include <geometry.hpp>
+#include <tc/groups.hpp>
+#include <iostream>
+#include <chrono>
+
+int main() {
+    tc::Group g = tc::group::B(3);
+    GeomGen m(g);
+
+    m.solve({}, {});
+    m.solve({0}, {});
+    m.solve({0}, {0});
+    m.solve({1}, {});
+    m.solve({1}, {1});
+    m.solve({2}, {});
+    m.solve({2}, {2});
+    m.solve({0, 1}, {});
+    m.solve({0, 1}, {0});
+    m.solve({0, 1}, {1});
+    m.solve({0, 1}, {0, 1});
+    m.solve({1, 2}, {});
+    m.solve({1, 2}, {1});
+    m.solve({1, 2}, {2});
+    m.solve({1, 2}, {1, 2});
+    m.solve({0, 2}, {});
+    m.solve({0, 2}, {0});
+    m.solve({0, 2}, {2});
+    m.solve({0, 2}, {0, 2});
+    m.solve({0, 1, 2}, {});
+    m.solve({0, 1, 2}, {0});
+    m.solve({0, 1, 2}, {1});
+    m.solve({0, 1, 2}, {2});
+    m.solve({0, 1, 2}, {0, 1});
+    m.solve({0, 1, 2}, {1, 2});
+    m.solve({0, 1, 2}, {0, 2});
+    m.solve({0, 1, 2}, {0, 1, 2});
+
+    tc::Group big = tc::group::B(8);
+    GeomGen mbig(big);
+
+    auto s1 = std::chrono::system_clock::now();
+    auto res1 = mbig.solve({0, 1, 2, 3, 4, 7}, {2, 4, 7});
+    auto e1 = std::chrono::system_clock::now();
+
+    std::chrono::duration<double> t1 = e1 - s1;
+    std::cout << t1.count() << ": " << res1.size() << std::endl;
+
+    auto s2 = std::chrono::system_clock::now();
+    auto res2 = mbig.solve({0, 2, 4, 7, 1, 3}, {4, 7, 2});
+    auto e2 = std::chrono::system_clock::now();
+
+    std::chrono::duration<double> t2 = e2 - s2;
+    std::cout << t2.count() << ": " << res2.size() << std::endl;
+
+    std::vector<int> gens = {0, 1, 2, 3, 4, 5};
+    auto s3 = std::chrono::system_clock::now();
+    auto res3 = mbig.triangulate(gens);
+    auto e3 = std::chrono::system_clock::now();
+
+    std::chrono::duration<double> t3 = e3 - s3;
+    std::cout << t3.count() << ": " << res3.size() << std::endl;
+
+    auto s4 = std::chrono::system_clock::now();
+    auto res4 = mbig.triangulate(gens);
+    auto e4 = std::chrono::system_clock::now();
+
+    std::chrono::duration<double> t4 = e4 - s4;
+    std::cout << t4.count() << ": " << res4.size() << std::endl;
+
+}

examples/subtest.cpp

@@ -0,0 +1,24 @@
+#include <iostream>
+#include <vector>
+
+#include <combo_iterator.hpp>
+#include <numeric>
+
+template<class T>
+std::ostream &operator<<(std::ostream &o, const std::vector<T> &v) {
+    for (const auto &e : v) o << e << " ";
+    return o;
+}
+
+int main() {
+    std::vector<int> gens(5);
+    std::iota(gens.begin(), gens.end(), 0);
+
+    const Combos<int> &combos = Combos(gens, 2);
+
+    for (const auto &e : combos) {
+        std::cout << e << std::endl;
+    }
+
+    return EXIT_SUCCESS;
+}

tc/CMakeLists.txt

@@ -1,6 +0,0 @@
-add_library(tc STATIC
-        src/groups.cpp
-        src/solver.cpp
-        src/cosets.cpp)
-
-target_include_directories(tc PUBLIC include)

tc/include/tc/cosets.h

@@ -1,33 +0,0 @@
-#pragma once
-
-#include "groups.h"
-#include <vector>
-
-namespace tc {
-    struct Action {
-        int coset = -1;
-        int gen = -1;
-        int target = -1;
-    };
-
-    struct Cosets {
-        int ngens;
-        std::vector<int> data;
-        std::vector<Action> path;
-        int len;
-
-        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]] int size() const;
-    };
-}

tc/include/tc/groups.h

@@ -1,96 +0,0 @@
-#pragma once
-
-#include <array>
-#include <vector>
-
-namespace tc {
-    struct Rel {
-        std::array<int, 2> gens;
-        int mult;
-    };
-
-    /**
-     * A presentation of a coxeter group. Contains a number of generators and some relations of the form (ab)^n = e
-     */
-    struct Group {
-        const int ngens;
-        std::vector<std::vector<int>> _mults;  // lookup table for multiplicities
-        std::string name;
-
-        explicit Group(int ngens, const std::vector<Rel> &rels = {}, std::string name = "G");
-
-        void setmult(Rel rel);
-
-        [[nodiscard]] Rel rel(int a, int b) const;
-
-        [[nodiscard]] std::vector<Rel> get_rels() const;
-
-        [[nodiscard]] Group product(const Group &other) const;
-
-        [[nodiscard]] Group power(int p) const;
-
-        [[nodiscard]] bool trivial() const;
-    };
-
-    Group operator*(const Group &g, const Group &h);
-
-    Group operator^(const Group &g, const int &p);
-
-    /**
-     * 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<int> &mults);
-
-    namespace group {
-        /**
-         * Simplex
-         */
-        Group A(int dim);
-
-        /**
-         * Cube, Orthoplex
-         */
-        Group B(int dim);
-
-        /**
-         * Demicube, Orthoplex
-         */
-        Group D(int dim);
-
-        /**
-         * E groups
-         */
-        Group E(int dim);
-
-        /**
-         * 24 Cell
-         */
-        Group F4();
-
-        /**
-         * Hexagon
-         */
-        Group G2();
-
-        /**
-         * Icosahedron
-         */
-        Group H(int dim);
-
-        /**
-         * Polygonal
-         */
-        Group I2(int n);
-
-        /**
-         * Toroidal. I2(n) * I2(m)
-         */
-        Group T(int n, int m);
-
-        /**
-         * Toroidal. T(n, n)
-         */
-        Group T(int n);
-    }
-}

tc/include/tc/solver.h

@@ -1,7 +0,0 @@
-#pragma once
-
-#include "cosets.h"
-
-namespace tc {
-    Cosets solve(const Group &g, const std::vector<int> &sub_gens = {});
-}

tc/src/cosets.cpp

@@ -1,44 +0,0 @@
-#include "tc/cosets.h"
-
-namespace tc {
-    Cosets::Cosets(int ngens) : ngens(ngens), len(0) {
-    }
-
-    void Cosets::add_row() {
-        len++;
-        data.resize(data.size() + ngens, -1);
-        path.resize(path.size() + 1);
-    }
-
-    void Cosets::put(int coset, int gen, int target) {
-        data[coset * ngens + gen] = target;
-        data[target * ngens + gen] = coset;
-
-        if (path[target].coset == -1) {
-            path[target] = {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[target].coset == -1) {
-            path[target] = {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];
-    }
-
-    int Cosets::size() const {
-        return len;
-    }
-}

tc/src/groups.cpp

@@ -1,207 +0,0 @@
-#include "tc/groups.h"
-
-#include <iterator>
-#include <sstream>
-
-namespace tc {
-    Group::Group(int ngens, const std::vector<Rel> &rels, std::string name) : ngens(ngens), name(std::move(name)) {
-        _mults.resize(ngens);
-        for (int i = 0; i < ngens; i++) {
-            _mults[i].resize(ngens, 2);
-        }
-
-        for (Rel rel : rels) {
-            setmult(rel);
-        }
-    }
-
-    void Group::setmult(Rel rel) {
-        _mults[rel.gens[0]][rel.gens[1]] = rel.mult;
-        _mults[rel.gens[1]][rel.gens[0]] = rel.mult;
-    }
-
-    Rel Group::rel(int a, int b) const {
-        return {a, b, _mults[a][b]};
-    }
-
-    std::vector<Rel> Group::get_rels() const {
-        std::vector<Rel> rels;
-        for (int i = 0; i < ngens - 1; i++) {
-            for (int j = i + 1; j < ngens; j++) {
-                rels.push_back({i, j, _mults[i][j]});
-            }
-        }
-        return rels;
-    }
-
-    Group Group::product(const Group &other) const {
-        int off = ngens;
-
-        Group g(ngens + other.ngens, get_rels());
-
-        for (Rel rel : other.get_rels()) {
-            g.setmult({off + rel.gens[0], off + rel.gens[1], rel.mult});
-        }
-
-        std::stringstream ss;
-        ss << name << "*" << other.name;
-        g.name = ss.str();
-
-        return g;
-    }
-
-    Group Group::power(int p) const {
-        Group g(ngens * p);
-
-        for (Rel rel : get_rels()) {
-            for (int off = 0; off < g.ngens; off += ngens) {
-                g.setmult({off + rel.gens[0], off + rel.gens[1], rel.mult});
-            }
-        }
-
-        std::stringstream ss;
-        ss << name << "^" << p;
-        g.name = ss.str();
-
-        return g;
-    }
-
-    bool Group::trivial() const {
-        for (int i = 0; i < ngens; ++i) {
-            for (int j = 0; j < ngens; ++j) {
-                if (_mults[i][j] != 2) {
-                    return false;
-                }
-            }
-        }
-
-        return true;
-    }
-
-    Group operator*(const Group &g, const Group &h) {
-        return g.product(h);
-    }
-
-    Group operator^(const Group &g, const int &p) {
-        return g.power(p);
-    }
-
-    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.setmult({i, i + 1, mults[i]});
-        }
-
-        return g;
-    }
-
-    Group schlafli(const std::vector<int> &mults) {
-        std::stringstream ss;
-        ss << "[";
-        if (!mults.empty()) {
-            copy(mults.begin(), mults.end() - 1, std::ostream_iterator<int>(ss, ","));
-            ss << mults.back();
-        }
-        ss << "]";
-
-        return schlafli(mults, ss.str());
-    }
-
-    namespace group {
-        Group A(const int dim) {
-            if (dim == 0)
-                return Group(0, {}, "A(0)");
-
-            const std::vector<int> &mults = std::vector<int>(dim - 1, 3);
-
-            std::stringstream ss;
-            ss << "A(" << dim << ")";
-
-            return schlafli(mults, ss.str());
-        }
-
-        Group B(const int dim) {
-            std::vector<int> mults(dim - 1, 3);
-            mults[0] = 4;
-
-            std::stringstream ss;
-            ss << "B(" << dim << ")";
-
-            return schlafli(mults, ss.str());
-        }
-
-        Group D(const int dim) {
-            std::vector<int> mults(dim - 1, 3);
-            mults[dim - 2] = 2;
-            Group g = schlafli(mults);
-            g.setmult({1, dim - 1, 3});
-
-            std::stringstream ss;
-            ss << "D(" << dim << ")";
-            g.name = ss.str();
-
-            return g;
-        }
-
-        Group E(const int dim) {
-            std::vector<int> mults(dim - 1, 3);
-            mults[dim - 2] = 2;
-            Group g = schlafli(mults);
-            g.setmult({2, dim - 1, 3});
-
-            std::stringstream ss;
-            ss << "E(" << dim << ")";
-            g.name = ss.str();
-
-            return g;
-        }
-
-        Group F4() {
-            return schlafli({3, 4, 3}, "F4");
-        }
-
-        Group G2() {
-            return schlafli({6}, "G2");
-        }
-
-        Group H(const int dim) {
-            std::vector<int> mults(dim - 1, 3);
-            mults[0] = 5;
-
-            std::stringstream ss;
-            ss << "H(" << dim << ")";
-
-            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) {
-            Group g = I2(n) * I2(m);
-
-            std::stringstream ss;
-            ss << "T(" << n << "," << m << ")";
-            g.name = ss.str();
-
-            return g;
-        }
-
-        Group T(const int n) {
-            Group g = I2(n) ^2;
-
-            std::stringstream ss;
-            ss << "T(" << n << ")";
-            g.name = ss.str();
-
-            return g;
-        }
-    }
-}

tc/src/solver.cpp

@@ -1,135 +0,0 @@
-#include "tc/solver.h"
-
-#include <algorithm>
-
-namespace tc {
-    struct RelTable {
-        Rel rel;
-
-        std::vector<int *> lst_ptr;
-        std::vector<int> gen;
-
-        int &mult = rel.mult;
-        std::array<int, 2> &gens = rel.gens;
-
-        explicit RelTable(Rel rel) : rel(rel) {
-        }
-
-        int add_row() {
-            int idx = lst_ptr.size();
-            lst_ptr.push_back(nullptr);
-            gen.push_back(-1);
-            return idx;
-        }
-    };
-
-    struct RelationSet {
-        const Cosets &cosets;
-
-        std::vector<RelTable> tables;
-        std::vector<std::vector<RelTable *>> gen_map; // which relations involve which generators
-
-        explicit RelationSet(const Group &g, const Cosets &cosets) : gen_map(g.ngens), cosets(cosets) {
-            const std::vector<Rel> &rels = g.get_rels();
-            tables.reserve(rels.size());
-            for (const auto &rel : rels) {
-                RelTable &table = tables.emplace_back(rel);
-                gen_map[rel.gens[0]].push_back(&table);
-                gen_map[rel.gens[1]].push_back(&table);
-            }
-        }
-
-        void add_row() {
-            for (auto &table : tables) {
-                table.add_row();
-            }
-        }
-
-        void fill_row(int idx) {
-            for (auto &table : tables) {
-                if (table.lst_ptr[idx] != nullptr) continue;
-
-                table.lst_ptr[idx] = new int;
-
-                if ((cosets.get(idx, table.gens[0]) != idx) and
-                    (cosets.get(idx, table.gens[1]) != idx)) {
-                    table.gen[idx] = 0;
-                } else {
-                    table.gen[idx] = -1;
-                }
-            }
-        }
-    };
-
-    Cosets solve(const Group &group, const std::vector<int> &sub_gens) {
-        Cosets cosets(group.ngens);
-        cosets.add_row();
-        for (const auto &i : sub_gens) {
-            cosets.put(0, i, 0);
-        }
-
-        RelationSet rels(group, cosets);
-        rels.add_row();
-        rels.fill_row(0);
-
-        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())
-                break;
-
-            target = cosets.len;
-
-            cosets.add_row();
-            rels.add_row();
-
-            std::vector<int> facts = {idx};
-
-            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 / group.ngens;
-                gen = fact_idx % group.ngens;
-
-                for (RelTable *pTable : rels.gen_map[gen]) {
-                    RelTable &table = *pTable;
-
-                    if (table.lst_ptr[target] == nullptr) {
-                        table.lst_ptr[target] = table.lst_ptr[coset];
-                        table.gen[target] = table.gen[coset] + 1;
-
-                        if (table.gen[coset] < 0)
-                            table.gen[target] -= 2;
-
-                        if (table.gen[target] == table.rel.mult) {
-                            lst = *(table.lst_ptr[target]);
-                            delete table.lst_ptr[target];
-                            gen_ = table.gens[table.gens[0] == gen];
-                            facts.push_back(lst * group.ngens + gen_);
-                        } else if (table.gen[target] == -table.mult) {
-                            gen_ = table.gens[table.gens[0] == gen];
-                            facts.push_back(target * group.ngens + gen_);
-                        } else if (table.gen[target] == table.mult - 1) {
-                            *(table.lst_ptr[target]) = target;
-                        }
-                    }
-                }
-
-                std::sort(facts.begin(), facts.end(), std::greater<>());
-            }
-
-            rels.fill_row(target);
-        }
-
-        return cosets;
-    }
-}

vis/CMakeLists.txt

@@ -1,7 +1,3 @@
-add_executable(vis src/main.cpp)
-target_include_directories(vis PRIVATE include)
-target_link_libraries(vis PRIVATE tc glad glm glfw)
-
 add_custom_target(shaders ALL DEPENDS shader_output)
 add_custom_command(
         OUTPUT shader_output
@@ -9,4 +5,10 @@ add_custom_command(
         COMMENT "copied shaders"
 )
 
+add_library(vis-util INTERFACE)
+target_include_directories(vis-util INTERFACE include)
+
+add_executable(vis src/main.cpp)
+target_include_directories(vis PRIVATE include)
+target_link_libraries(vis PRIVATE tc glad glm glfw)
 add_dependencies(vis shaders)

vis/include/combo_iterator.hpp

@@ -0,0 +1,94 @@
+#pragma once
+
+#include <vector>
+#include <algorithm>
+#include <stdexcept>
+
+template<class T>
+struct ComboIterator {
+    typedef ComboIterator<T> self_type;
+    typedef const std::vector<T> value_type;;
+    typedef const std::vector<T> *pointer;
+    typedef const std::vector<T> &reference;
+    typedef size_t difference_type;
+    typedef std::forward_iterator_tag iterator_category;
+
+    const std::vector<T> &vals;
+    const size_t k;
+    size_t n;
+    std::vector<bool> bits;
+    std::vector<T> curr;
+
+    ComboIterator(ComboIterator &) = default;
+
+    ComboIterator(const std::vector<T> &vals, const size_t k, const size_t n)
+        : vals(vals), k(k), n(n), curr(k), bits(vals.size()) {
+        for (size_t i = 0; i < vals.size(); ++i) {
+            bits[i] = i < k;
+        }
+        std::reverse(bits.begin(), bits.end());
+        set_curr();
+    }
+
+    ~ComboIterator() = default;
+
+    void set_curr() {
+        for (size_t i = 0, j = 0; i < vals.size(); ++i) {
+            if (bits[i]) curr[j++] = vals[i];
+        }
+    }
+
+    [[nodiscard]] bool operator==(const ComboIterator<T> &o) const {
+        return n == o.n;
+    }
+
+    [[nodiscard]] bool operator!=(const ComboIterator<T> &o) const {
+        return n != o.n;
+    }
+
+    reference operator*() const {
+        return curr;
+    }
+
+    pointer operator->() const {
+        return &curr;
+    }
+
+    self_type operator++(int) {
+        std::next_permutation(bits.begin(), bits.end());
+        set_curr();
+        ++n;
+        return *this;
+    }
+
+    self_type operator++() {
+        self_type r = *this;
+        std::next_permutation(bits.begin(), bits.end());
+        set_curr();
+        ++n;
+        return r;
+    }
+};
+
+size_t choose(size_t n, size_t k) {
+    if (k == 0) return 1;
+    return n * choose(n - 1, k - 1) / k;
+}
+
+template<class T>
+struct Combos {
+    const std::vector<T> &vals;
+    size_t k;
+
+    Combos(const std::vector<T> &vals, size_t k) : vals(vals), k(k) {
+    }
+
+    [[nodiscard]] ComboIterator<T> begin() const {
+        return ComboIterator(vals, k, 0);
+    }
+
+    [[nodiscard]] ComboIterator<T> end() const {
+        int j = choose(vals.size(), k);
+        return ComboIterator(vals, k, j);
+    }
+};

vis/include/geometry.hpp

@@ -0,0 +1,284 @@
+#pragma once
+
+#include <tc/core.hpp>
+#include <cmath>
+#include <optional>
+#include <numeric>
+#include <iostream>
+#include "combo_iterator.hpp"
+
+size_t get_key_from_gens(std::vector<int> &gens) {
+    size_t key = 0;
+    for (const auto gen : gens) {
+        key += (1u << (unsigned) gen);
+    }
+    return key;
+}
+
+size_t num_gens_from_key(size_t key) {
+    size_t mask = 1;
+    size_t count = 0;
+    while (mask <= key) {
+        if (key & mask)
+            count++;
+        mask <<= 1u;
+    }
+    return count;
+}
+
+template<unsigned N>
+struct Primitive {
+    std::array<unsigned, N> inds;
+
+    Primitive() = default;
+
+    Primitive(const Primitive<N> &) = default;
+
+    Primitive(const Primitive<N - 1> &sub, unsigned root) {
+        std::copy(sub.inds.begin(), sub.inds.end(), inds.begin());
+        inds[N - 1] = root;
+    }
+
+    ~Primitive() = default;
+
+    inline void flip() {
+        if (N > 1) std::swap(inds[0], inds[1]);
+    }
+
+    void apply(const tc::Cosets &table, int gen) {
+        for (auto &ind : inds) {
+            ind = table.get(ind, gen);
+        }
+        flip();
+    }
+};
+
+std::vector<int> gens(const tc::Group &context) {
+    std::vector<int> g_gens(context.ngens);
+    std::iota(g_gens.begin(), g_gens.end(), 0);
+    return g_gens;
+}
+
+std::vector<int> recontext_gens(
+    const tc::Group &context,
+    std::vector<int> g_gens,
+    std::vector<int> sg_gens) {
+
+    std::sort(g_gens.begin(), g_gens.end());
+
+    int inv_gen_map[context.ngens];
+    for (size_t i = 0; i < g_gens.size(); i++) {
+        inv_gen_map[g_gens[i]] = i;
+    }
+
+    std::vector<int> 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]);
+    }
+    std::sort(s_sg_gens.begin(), s_sg_gens.end());
+
+    return s_sg_gens;
+
+//    std::sort(g_gens.begin(), g_gens.end());
+//
+//    std::vector<int> inv_g_map(g_gens.size());
+//    for (int i = 0; i < g_gens.size(); ++i) {
+//        inv_g_map[g_gens[i]] = i;
+//    }
+//
+//    std::transform(sg_gens.begin(), sg_gens.end(), sg_gens.begin(),
+//        [inv_g_map](const auto &gen) {
+//            return inv_g_map[gen];
+//        }
+//    );
+//
+//    std::sort(sg_gens.begin(), sg_gens.end());
+//    return sg_gens;
+}
+
+int get_parity(
+    const tc::Group &context,
+    const std::vector<int> &g_gens,
+    const std::vector<int> &sg_gens
+) {
+    if (g_gens.size() != sg_gens.size() + 1) return 0;
+
+    const auto proper_sg_gens = recontext_gens(context, g_gens, sg_gens);
+
+    int i = 0;
+    for (; i < sg_gens.size(); ++i) {
+        if (proper_sg_gens[i] != i) {
+            break;
+        }
+    }
+
+    return i & 1;
+}
+
+tc::Cosets solve(
+    const tc::Group &context,
+    const std::vector<int> &g_gens,
+    const std::vector<int> &sg_gens
+) {
+    const auto proper_sg_gens = recontext_gens(context, g_gens, sg_gens);
+    return context.subgroup(g_gens).solve(proper_sg_gens);
+}
+
+tc::Cosets solve_sg(
+    const tc::Group &context,
+    const std::vector<int> &sg_gens
+) {
+    return solve(context, gens(context), sg_gens);
+}
+
+tc::Cosets solve_g(
+    const tc::Group &context,
+    const std::vector<int> &g_gens
+) {
+    std::vector<int> sg_gens;
+    return solve(context, g_gens, sg_gens);
+}
+
+tc::Cosets solve(
+    const tc::Group &context
+) {
+    std::vector<int> sg_gens;
+    return solve_sg(context, sg_gens);
+}
+
+template<unsigned N>
+struct Mesh {
+    std::vector<Primitive<N>> prims;
+
+    Mesh() : prims() {}
+
+    Mesh(const Mesh<N> &) = default;
+
+    explicit Mesh(std::vector<Primitive<N>> &prims) : prims(prims) {}
+
+    [[nodiscard]] size_t size() const {
+        return prims.size();
+    }
+
+    void apply(const tc::Cosets &table, int gen) {
+        for (auto &prim : prims) {
+            prim.apply(table, gen);
+        }
+    }
+
+    void flip() {
+        for (auto &prim : prims) {
+            prim.flip();
+        }
+    }
+
+    [[nodiscard]]
+    Mesh<N> recontext(
+        const tc::Group &context,
+        const std::vector<int> &g_gens,
+        const std::vector<int> &sg_gens
+    ) const {
+        const auto proper_sg_gens = recontext_gens(context, g_gens, sg_gens);
+        const auto table = solve_g(context, g_gens);
+        const auto path = solve_g(context, sg_gens).path;
+
+        auto map = path.template walk<int, int>(0, proper_sg_gens, [table](int coset, int gen) {
+            return table.get(coset, gen);
+        });
+
+        Mesh<N> res = *this;
+        for (Primitive<N> &prim : res.prims) {
+            for (auto &ind : prim.inds) {
+                ind = map[ind];
+            }
+        }
+
+        if (get_parity(context, g_gens, sg_gens) == 1)
+            res.flip();
+
+        return res;
+    }
+
+    [[nodiscard]]
+    Mesh<N> tile(
+        const tc::Group &context,
+        const std::vector<int> &g_gens,
+        const std::vector<int> &sg_gens
+    ) const {
+        Mesh<N> base = recontext(context, g_gens, sg_gens);
+        const auto proper_sg_gens = recontext_gens(context, g_gens, sg_gens);
+
+        const auto table = solve_g(context, g_gens);
+        const auto path = solve(context, g_gens, sg_gens).path;
+
+        const auto all = path.template walk<Mesh<N>, int>(base, gens(context), [table](Mesh<N> from, int gen) {
+            from.apply(table, gen);
+            return from;
+        });
+
+        return merge(all);
+    }
+
+    [[nodiscard]]
+    Mesh<N + 1> fan(int root) const {
+        std::vector<Primitive<N + 1>> res(prims.size());
+        std::transform(prims.begin(), prims.end(), res.begin(),
+            [root](const Primitive<N> &prim) {
+                return Primitive<N + 1>(prim, root);
+            }
+        );
+        return Mesh<N + 1>(res);
+    }
+};
+
+template<unsigned N>
+Mesh<N> merge(const std::vector<Mesh<N>> &meshes) {
+    size_t size = 0;
+    for (const auto &mesh : meshes) {
+        size += mesh.size();
+    }
+
+    std::vector<Primitive<N>> prims;
+    prims.reserve(size);
+    for (const auto &mesh : meshes) {
+        prims.insert(prims.end(), mesh.prims.begin(), mesh.prims.end());
+    }
+
+    return Mesh(prims);
+}
+
+template<unsigned N>
+Mesh<N> triangulate(
+    const tc::Group &context,
+    const std::vector<int> &g_gens
+) {
+    if (g_gens.size() + 1 != N)
+        throw std::logic_error("g_gens size must be one less than N");
+
+    const auto &combos = Combos(g_gens, g_gens.size() - 1);
+
+    std::vector<Mesh<N>> meshes;
+    for (const auto &sg_gens : combos) {
+        Mesh<N - 1> base = triangulate<N - 1>(context, sg_gens);
+        Mesh<N - 1> raised = base.tile(context, g_gens, sg_gens);
+        raised.prims.erase(raised.prims.begin(), raised.prims.begin() + base.size());
+        Mesh<N> fan = raised.fan(0);
+        meshes.push_back(fan);
+    }
+
+    return merge(meshes);
+}
+
+template<>
+Mesh<1> triangulate<1>(
+    const tc::Group &context,
+    const std::vector<int> &g_gens
+) {
+    if (not g_gens.empty())
+        throw std::logic_error("g_gens must be empty for a trivial Mesh");
+
+    Mesh<1> res;
+    res.prims.emplace_back();
+    return res;
+}

vis/include/mirror.hpp

@@ -1,6 +1,6 @@
 #pragma once
 
-#include <tc/cosets.h>
+#include <tc/core.hpp>
 #include <glm/glm.hpp>
 #include <cmath>
 #include <vector>
@@ -20,28 +20,49 @@ float dot(int n, const glm::vec4 &a, const glm::vec4 &b) {
     return sum;
 }
 
+float dot(int n, const std::vector<float> &a, const std::vector<float> &b) {
+    float sum = 0;
+    for (int i = 0; i < n; ++i) {
+        sum += a[i] * b[i];
+    }
+    return sum;
+}
+
 std::vector<glm::vec4> mirror(const tc::Group &group) {
-    std::vector<glm::vec4> mirrors;
+    std::vector<std::vector<float>> mirrors;
 
     for (int p = 0; p < group.ngens; ++p) {
-        glm::vec4 vp{};
+        std::vector<float> vp;
         for (int m = 0; m < p; ++m) {
-            glm::vec4 vq = mirrors[m];
+            auto &vq = mirrors[m];
-            vp[m] = (cos(M_PI / group.rel(p, m).mult) - dot(m, vp, vq)) / vq[m];
+            vp.push_back((cos(M_PI / group.get(p, m)) - dot(m, vp, vq)) / vq[m]);
         }
-        vp[p] = std::sqrt(1 - glm::dot(vp, vp));
+        vp.push_back(std::sqrt(1 - dot(p, vp, vp)));
 
         for (const auto &v : mirrors) {
-            if (glm::dot(vp, v) > 0) {
+            if (dot(p, vp, vp) > 0) {
-                vp *= -1;
+                for (auto &e : vp) {
+                    e *= -1;
+                }
                 break;
             }
         }
 
-        mirrors.push_back(round(vp, 15));
+        mirrors.push_back(vp);
     }
 
-    return mirrors;
+    std::vector<glm::vec4> res;
+    for (const auto &vec : mirrors) {
+        glm::vec4 rvec{};
+
+        // ortho proj
+        for (int i = 0; i < std::min(vec.size(), (size_t) 4); ++i) {
+            rvec[i] = vec[i];
+        }
+
+        res.push_back(rvec);
+    }
+    return res;
 }
 
 glm::vec4 project(const glm::vec4 &vec, const glm::vec4 &target) {
@@ -52,6 +73,10 @@ glm::vec4 reflect(const glm::vec4 &vec, const glm::vec4 axis) {
     return vec - 2.f * project(vec, axis);
 }
 
+glm::vec4 reflect_scaled(const  glm::vec4 &vec, const glm::vec4 axis) {
+    return vec - 2.f * glm::length(axis) * project(vec, axis);
+}
+
 glm::vec4 gram_schmidt_last(std::vector<glm::vec4> vecs) {
     int N = vecs.size();
     for (int i = 0; i < N; ++i) {
@@ -84,3 +109,12 @@ std::vector<glm::vec4> plane_intersections(std::vector<glm::vec4> normals) {
 
     return results;
 }
+
+glm::mat4 utilRotate(const int u, const int v, const float theta) {
+    auto res = glm::identity<glm::mat4>();
+    res[u][u] = std::cos(theta);
+    res[u][v] = std::sin(theta);
+    res[v][u] = -std::sin(theta);
+    res[v][v] = std::cos(theta);
+    return res;
+}

vis/include/util.hpp

@@ -0,0 +1,161 @@
+#pragma once
+
+#include <cerrno>
+#include <fstream>
+#include <sstream>
+#include <string>
+#include <vector>
+
+#include <glad/glad.h>
+
+class gl_error : public std::domain_error {
+public:
+    explicit gl_error(const std::string &arg) : domain_error(arg) {}
+
+    explicit gl_error(const char *string) : domain_error(string) {}
+};
+
+class shader_error : public gl_error {
+public:
+    explicit shader_error(const std::string &arg) : gl_error(arg) {}
+
+    explicit shader_error(const char *string) : gl_error(string) {}
+};
+
+class program_error : public gl_error {
+public:
+    explicit program_error(const std::string &arg) : gl_error(arg) {}
+
+    explicit program_error(const char *string) : gl_error(string) {}
+};
+
+template<class T, GLenum prop>
+T utilGetShader(GLuint shader) {
+    GLint res;
+    glGetShaderiv(shader, prop, &res);
+    return static_cast<T>(res);
+}
+
+#define getShaderInfoLogLength utilGetShader<size_t, GL_INFO_LOG_LENGTH>
+#define getShaderCompileStatus utilGetShader<size_t, GL_COMPILE_STATUS>
+
+template<class T, GLenum prop>
+T utilGetProgram(GLuint program) {
+    GLint res;
+    glGetProgramiv(program, prop, &res);
+    return static_cast<T>(res);
+}
+
+#define getProgramInfoLogLength utilGetProgram<size_t, GL_INFO_LOG_LENGTH>
+#define getProgramLinkStatus utilGetProgram<size_t, GL_LINK_STATUS>
+
+void utilShaderSource(GLuint shader, const std::vector<std::string> &sources) {
+    char const *ptrs[sources.size()];
+    for (size_t i = 0; i < sources.size(); ++i) {
+        ptrs[i] = sources[i].c_str();
+    }
+    glShaderSource(shader, sources.size(), ptrs, nullptr);
+}
+
+std::string getShaderInfoLog(GLuint shader) {
+    int len = getShaderInfoLogLength(shader);
+    char buffer[len];
+    glGetShaderInfoLog(shader, len, nullptr, buffer);
+    return std::string(buffer);
+}
+
+std::string getProgramInfoLog(GLuint program) {
+    int len = getProgramInfoLogLength(program);
+    char buffer[len];
+    glGetProgramInfoLog(program, len, nullptr, buffer);
+    return std::string(buffer);
+}
+
+std::string utilInfo() {
+    std::stringstream ss;
+    ss
+        << "Graphics Information:" << std::endl
+        << "    Vendor:          " << glGetString(GL_VENDOR) << std::endl
+        << "    Renderer:        " << glGetString(GL_RENDERER) << std::endl
+        << "    OpenGL version:  " << glGetString(GL_VERSION) << std::endl
+        << "    Shading version: " << glGetString(GL_SHADING_LANGUAGE_VERSION) << std::endl;
+    return ss.str();
+}
+
+std::string utilReadFile(const std::string &filename) {
+    std::ifstream in(filename, std::ios::in | std::ios::binary);
+    if (in) {
+        std::ostringstream contents;
+        contents << in.rdbuf();
+        in.close();
+        return (contents.str());
+    }
+    throw std::system_error(errno, std::generic_category());
+}
+
+GLuint utilCompileFiles(const GLenum type, const std::vector<std::string> &files) {
+    std::vector<std::string> sources;
+    sources.reserve(files.size());
+    for (const auto &file : files) {
+        sources.push_back(utilReadFile(file));
+    }
+
+    GLuint shader = glCreateShader(type);
+    utilShaderSource(shader, sources);
+    glCompileShader(shader);
+
+    if (getShaderCompileStatus(shader)) return shader;
+
+    throw shader_error(getShaderInfoLog(shader));
+}
+
+GLuint utilLinkProgram(const std::vector<GLuint> &shaders) {
+    GLuint program = glCreateProgram();
+    for (const auto &shader : shaders) {
+        glAttachShader(program, shader);
+    }
+    glLinkProgram(program);
+
+    if (getProgramLinkStatus(program)) return program;
+
+    throw program_error(getProgramInfoLog(program));
+}
+
+GLuint utilCreateShaderProgram(GLenum type, const std::vector<std::string> &src) {
+    std::vector<const char *> c_str(src.size());
+    std::transform(src.begin(), src.end(), c_str.begin(), [](auto &str) {
+        return str.c_str();
+    });
+
+    GLuint program = glCreateShaderProgramv(type, src.size(), &c_str[0]);
+
+    if (getProgramLinkStatus(program)) return program;
+
+    throw program_error(getProgramInfoLog(program));
+}
+
+GLuint utilCreateShaderProgramFile(GLenum type, const std::vector<std::string> &files) {
+    std::vector<std::string> sources(files.size());
+    std::transform(files.begin(), files.end(), sources.begin(), utilReadFile);
+    return utilCreateShaderProgram(type, sources);
+}
+
+std::vector<GLuint> utilCreateVertexArrays(int n) {
+    std::vector<GLuint> res(n);
+    glCreateVertexArrays(n, &res[0]);
+    return res;
+}
+
+GLuint utilCreateVertexArray() {
+    return utilCreateVertexArrays(1)[0];
+}
+
+std::vector<GLuint> utilCreateBuffers(int n) {
+    std::vector<GLuint> res(n);
+    glCreateBuffers(n, &res[0]);
+    return res;
+}
+
+GLuint utilCreateBuffer() {
+    return utilCreateBuffers(1)[0];
+}

vis/shaders/4d/4d.gm.glsl

@@ -0,0 +1,65 @@
+#version 440 core
+
+layout(points) in;
+layout(triangle_strip, max_vertices=4) out;
+
+layout(std430, binding=1) buffer Positions {
+    vec4 verts[];
+};
+
+layout(std140, binding=1) uniform Matrices {
+    mat4 proj;
+    mat4 view;
+};
+
+in ivec4 vInds[];
+
+out vec4 pos;
+
+out gl_PerVertex {
+    vec4 gl_Position;
+};
+
+float unmix(float u, float v) {
+    return (u) / (u - v);
+}
+
+void emit(vec4 v) {
+    pos = v;
+    gl_Position = proj * vec4(v.xyz, 1);
+    EmitVertex();
+}
+
+void main() {
+    vec4 pos4[4];
+    for (int i = 0; i < 4; ++i) pos4[i] = view * verts[vInds[0][i]];
+
+    int lo[4], L = 0;
+    int hi[4], H = 0;
+
+    float x = 0.7;
+
+    for (int i = 0; i < 4; ++i) {
+        if (pos4[i].w < 0) {
+            lo[L++] = i;
+        } else {
+            hi[H++] = i;
+        }
+    }
+
+    vec4 sect[4]; int S = 0;
+    for (int l = 0; l < L; ++l)
+    for (int h = H-1; h >=0; --h) {
+        vec4 a = pos4[lo[l]];
+        vec4 b = pos4[hi[h]];
+
+        float t = unmix(a.w, b.w);
+        sect[S++] = mix(a, b, t);
+    }
+
+    for (int s = 0; s < S; ++s) {
+        emit(sect[s]);
+    }
+
+    EndPrimitive();
+}

vis/shaders/4d/4d.vs.glsl

@@ -0,0 +1,27 @@
+#version 440 core
+
+layout(std430, binding=1) buffer Positions {
+    vec4 verts[];
+};
+
+layout(std140, binding=1) uniform Matrices {
+    mat4 proj;
+    mat4 view;
+};
+
+layout(location=0) in ivec4 inds;
+
+out ivec4 vInds;
+
+out gl_PerVertex {
+    vec4 gl_Position;
+    float gl_PointSize;
+};
+
+void main() {
+    vInds = inds;
+
+    vec4 pos = view * verts[vInds.x];
+    gl_Position = proj * vec4(pos.xyz, 1);
+    gl_PointSize = 5;
+}

vis/shaders/one-color.fs.glsl

@@ -0,0 +1,12 @@
+#version 430
+
+layout(location=2) uniform vec3 c;
+
+in vec4 pos;
+
+out vec4 color;
+
+void main() {
+    float d = smoothstep(-2, 2, pos.z);
+    color = vec4(c * d, 1);
+}

vis/shaders/ortho.vs.glsl

@@ -0,0 +1,17 @@
+#version 430
+
+out gl_PerVertex {
+    vec4 gl_Position;
+};
+
+layout(location=0) uniform mat4 proj;
+layout(location=1) uniform mat4 view;
+
+layout(location=0) in vec4 pos;
+
+out vec4 vpos;
+
+void main() {
+    vpos = view * pos;
+    gl_Position = proj * vec4(vpos.xyz / (1), 1);
+}

vis/shaders/stereo-proper.gm.glsl

@@ -0,0 +1,25 @@
+#version 430
+
+#define SUBS 20
+
+layout(lines) in;
+layout(line_strip, max_vertices = SUBS) out;
+
+layout(location=0) uniform mat4 proj;
+
+in vec4 gpos[];
+out vec4 vpos;
+
+vec4 stereo(vec4 v) {
+    return vec4(v.xyz / (1 - v.w), 1);
+}
+
+void main() {
+    for (int i = 0; i < SUBS; i++) {
+        vpos = mix(gpos[0], gpos[1], i * 1.0f / (SUBS - 1));
+        vpos = normalize(vpos);
+        gl_Position = proj * stereo(vpos);
+        EmitVertex();
+    }
+    EndPrimitive();
+}

vis/shaders/stereo-proper.vs.glsl

@@ -0,0 +1,14 @@
+#version 430
+
+layout(location=0) uniform mat4 proj;
+layout(location=1) uniform mat4 view;
+
+layout(location=0) in vec4 pos;
+
+out vec4 vpos;
+
+void main() {
+    vpos = view * pos;
+    gl_Position = proj * vec4(vpos.xyz / (1), 1);
+    gl_PointSize = 5;
+}

vis/shaders/stereo.vs.glsl

@@ -0,0 +1,15 @@
+#version 430
+
+layout(location=0) uniform mat4 proj;
+layout(location=1) uniform mat4 view;
+
+layout(location=0) in vec4 pos;
+
+out vec4 vpos;
+
+void main() {
+    int i = gl_VertexID;
+    vpos = view * pos;
+     gl_Position = proj * vec4(vpos.xyz / (1 - vpos.w), 1);
+    gl_PointSize = 5 * smoothstep(-2, 2, gl_Position.z);
+}

vis/shaders/w-axis-hue.fs.glsl

@@ -0,0 +1,11 @@
+#version 430
+
+in vec4 vpos;
+
+out vec4 color;
+
+void main() {
+    float d = smoothstep(-2, 2, vpos.z);
+    vec3 off = 1.04 * vec3(0, 2, 4) + 2 * vec3(vpos.w);
+    color = vec4(d * cos(off), 1);
+}

vis/src/main.cpp

@@ -1,14 +1,15 @@
 #include <glad/glad.h>
 #include <GLFW/glfw3.h>
+#include <cmath>
 #include <iostream>
 
-#include <tc/groups.h>
-#include <tc/solver.h>
-
-#include "geom.h"
-
 #include <glm/gtc/type_ptr.hpp>
-#include <glm/gtx/string_cast.hpp>
+
+#include <tc/groups.hpp>
+
+#include "util.hpp"
+#include "mirror.hpp"
+#include "geometry.hpp"
 
 #ifdef _WIN32
 extern "C" {
@@ -16,92 +17,33 @@ __attribute__((unused)) __declspec(dllexport) int NvOptimusEnablement = 0x000000
 }
 #endif
 
-void utilShaderSource(GLuint shader, const std::vector<std::string> &sources) {
+struct Matrices {
-    char const *ptrs[sources.size()];
+    glm::mat4 proj;
-    for (size_t i = 0; i < sources.size(); ++i) {
+    glm::mat4 view;
-        ptrs[i] = sources[i].c_str();
+};
-    }
-    glShaderSource(shader, sources.size(), ptrs, nullptr);
-}
 
-std::string utilShaderInfoLog(GLuint shader) {
+Matrices build(GLFWwindow *window, float st) {
-    int len;
+    int width, height;
-    glGetShaderiv(shader, GL_INFO_LOG_LENGTH, &len);
+    glfwGetFramebufferSize(window, &width, &height);
-    char buffer[len];
-    glGetShaderInfoLog(shader, len, nullptr, buffer);
-    return std::string(buffer);
-}
 
-std::string utilProgramInfoLog(GLuint program) {
+    auto aspect = (float) width / (float) height;
-    int len;
+    auto pheight = 1.4f;
-    glGetProgramiv(program, GL_INFO_LOG_LENGTH, &len);
+    auto pwidth = aspect * pheight;
-    char buffer[len];
+    glm::mat4 proj = glm::ortho(-pwidth, pwidth, -pheight, pheight, -10.0f, 10.0f);
-    glGetProgramInfoLog(program, len, nullptr, buffer);
-    return std::string(buffer);
-}
 
+    auto view = glm::identity<glm::mat4>();
+    view *= utilRotate(0, 1, st * .40f);
+    view *= utilRotate(0, 2, st * .20f);
+    view *= utilRotate(0, 3, st * 1.30f);
+    view *= utilRotate(1, 2, st * .50f);
+    view *= utilRotate(1, 3, st * .25f);
+    view *= utilRotate(2, 3, st * 1.42f);
 
-template<int D>
+    return {proj, view};
-tc::Group shrink(const tc::Group &g, const std::array<int, D> &sub) {
-    tc::Group h(D);
-    for (int i = 0; i < D; ++i) {
-        for (int j = 0; j < D; ++j) {
-            h.setmult({i, j, g.rel(sub[i], sub[j]).mult});
 }
-    }
-    return h;
-}
-
-template<int D>
-std::vector<int> raise(
-    const tc::Cosets &g_res,
-    const tc::Cosets &h_res,
-    const std::array<int, D> &gen_map,
-    const std::vector<int> &path
-) {
-    std::vector<int> res(path.size(), 0);
-    for (size_t i = 1; i < path.size(); ++i) {
-        auto action = h_res.path[path[i]];
-
-        res[i] = g_res.get(res[action.coset], gen_map[action.gen]);
-    }
-    return res;
-}
-
-std::vector<int> targets(const std::vector<tc::Action> &path) {
-    std::vector<int> res(path.size(), 0);
-    for (size_t i = 0; i < path.size(); ++i) {
-        res[i] = path[i].target;
-    }
-    return res;
-}
-
-std::vector<int> tile(const tc::Cosets &map, std::vector<int> geom) {
-    int K = geom.size();
-    geom.resize(K * map.size());
-    for (int i = 1; i < map.size(); ++i) { // tile
-        auto gaction = map.path[i];
-        for (int j = 0; j < K; ++j) {
-            geom[i * K + j] = map.get(geom[gaction.coset * K + j], gaction.gen);
-        }
-    }
-    return geom;
-}
-
-template<int D>
-std::vector<int> build(const tc::Group &g, const std::array<int, D> &sub) {
-    tc::Group h = shrink<D>(g, sub);
-    auto hres = tc::solve(h); // recursion would happen here
-
-    auto gres = tc::solve(g);
-    std::vector<int> geom = targets(hres.path);
-    geom = raise<D>(gres, hres, sub, geom);
-
-    return tile(gres, geom);
-}
-
 
 int main(int argc, char *argv[]) {
+    //region init window
     if (!glfwInit()) {
         std::cerr << "Failed to initialize GLFW" << std::endl;
         return EXIT_FAILURE;
@@ -121,175 +63,126 @@ int main(int argc, char *argv[]) {
     glfwMakeContextCurrent(window);
     gladLoadGLLoader((GLADloadproc) glfwGetProcAddress);
     glfwSwapInterval(0);
+    glClear(GL_COLOR_BUFFER_BIT);
+    glfwSwapBuffers(window);
+    //endregion
 
-    std::cout
+    std::cout << utilInfo();
-        << "Graphics Information:" << std::endl
-        << "    Vendor:          " << glGetString(GL_VENDOR) << std::endl
-        << "    Renderer:        " << glGetString(GL_RENDERER) << std::endl
-        << "    OpenGL version:  " << glGetString(GL_VERSION) << std::endl
-        << "    Shading version: " << glGetString(GL_SHADING_LANGUAGE_VERSION) << std::endl;
 
-    auto group = tc::group::B(2);
+    glEnable(GL_PROGRAM_POINT_SIZE);
-    auto res = tc::solve(group);
+    glEnable(GL_POINT_SMOOTH);
-    auto mirrors = mirror(group);
+    glEnable(GL_DEPTH_TEST);
-    auto corners = plane_intersections(mirrors);
+//    glEnable(GL_CULL_FACE);
+    glCullFace(GL_BACK);
 
-    auto points = std::vector<glm::vec4>(res.size());
+    //region shaders
-    points[0] = barycentric(corners, {1.00f, 0.50f, 0.50f, 0.50f});
+    GLuint pipe;
-    for (int i = 1; i < res.size(); ++i) {
+    glCreateProgramPipelines(1, &pipe);
-        auto action = res.path[i];
-        points[i] = reflect(points[action.coset], mirrors[action.gen]);
-    }
 
-    auto h = shrink<2>(group, {0, 1});
+    GLuint vs, gm, fs;
-    auto i = shrink<1>(h, {0});
 
-    auto g_map = res;
+    try {
-    auto h_map = tc::solve(h);
+        vs = utilCreateShaderProgramFile(GL_VERTEX_SHADER, {"shaders/4d/4d.vs.glsl"});
-    auto i_map = tc::solve(i);
+        gm = utilCreateShaderProgramFile(GL_GEOMETRY_SHADER, {"shaders/4d/4d.gm.glsl"});
+        fs = utilCreateShaderProgramFile(GL_FRAGMENT_SHADER, {"shaders/one-color.fs.glsl"});
 
-    auto g0 = targets(i_map.path);
+        glUseProgramStages(pipe, GL_VERTEX_SHADER_BIT, vs);
-    g0 = raise<1>(h_map, i_map, {0}, g0);
+        glUseProgramStages(pipe, GL_GEOMETRY_SHADER_BIT, gm);
-    g0 = tile(h_map, g0);
+        glUseProgramStages(pipe, GL_FRAGMENT_SHADER_BIT, fs);
-    g0 = raise<2>(g_map, h_map, {0, 1}, g0);
+    } catch (const gl_error &e) {
-
+        std::cerr << e.what() << std::endl;
-//    auto g0 = build<2>(group, {0});
-//    auto g1 = build<2>(group, {1});
-//    auto g2 = build<2>(group, {2});
-//    auto g3 = build<2>(group, {3});
-
-    GLuint vs = glCreateShader(GL_VERTEX_SHADER);
-    utilShaderSource(vs, {
-        "#version 430\n",
-
-        "layout(location=0) uniform mat4 proj;"
-        "layout(location=1) uniform mat4 view;"
-        ""
-        "layout(location=0) in vec4 pos;"
-        ""
-        "out vec4 vpos;"
-        ""
-        "void main() {"
-        "   int i = gl_VertexID;"
-        "   vpos = view * pos;"
-        //        "   gl_Position = proj * vec4(vpos.xyz / (1), 1);"
-        "   gl_Position = proj * vec4(vpos.xyz / (1 - vpos.w), 1);"
-        "   gl_PointSize = 5;"
-        "}"
-    });
-    glCompileShader(vs);
-
-    GLuint fs = glCreateShader(GL_FRAGMENT_SHADER);
-    utilShaderSource(fs, {
-        "#version 430\n",
-
-        "layout(location=2) uniform vec3 c;"
-        ""
-        "in vec4 vpos;"
-        ""
-        "out vec4 color;"
-        ""
-        "void main() {"
-        "   float d = smoothstep(-2, 2, vpos.z);"
-        "   color = vec4(c * d, 1);"
-        "}"
-    });
-
-    GLuint pgm = glCreateProgram();
-    glAttachShader(pgm, vs);
-    glAttachShader(pgm, fs);
-    glLinkProgram(pgm);
-
-    GLint status;
-
-    glGetShaderiv(vs, GL_COMPILE_STATUS, &status);
-    if (!status) {
-        std::cerr << utilShaderInfoLog(vs) << "\n=========\n" << std::endl;
-    }
-
-    glGetShaderiv(fs, GL_COMPILE_STATUS, &status);
-    if (!status) {
-        std::cerr << utilShaderInfoLog(fs) << "\n=========\n" << std::endl;
-    }
-
-    glGetProgramiv(pgm, GL_LINK_STATUS, &status);
-    if (!status) {
-        std::cerr << utilProgramInfoLog(pgm) << "\n=========\n" << std::endl;
         glfwTerminate();
-        return EXIT_FAILURE;
+        exit(EXIT_FAILURE);
     }
+    //endregion
+
+    //region points
+    auto group = tc::group::H(4);
+    auto res = group.solve();
+    auto mirrors = mirror(group);
+
+    auto corners = plane_intersections(mirrors);
+//    auto start = barycentric(corners, {1.0f, 1.0f, 1.0f, 1.0f});
+    auto start = barycentric(corners, {1.00f, 0.2f, 0.1f, 0.05f});
+//    auto start = barycentric(corners, {0.05f, 0.1f, 0.2f, 1.00f});
+    auto points = res.path.walk<glm::vec4, glm::vec4>(start, mirrors, reflect);
+
+    auto g_gens = gens(group);
+
+    std::vector<GLuint> vaos;
+    std::vector<GLuint> ibos;
+    std::vector<unsigned> counts;
+
+    auto combos = Combos(g_gens, 3);
+//    std::vector<std::vector<int>> chosen = {
+//        {1, 2, 3},
+//        {0, 2, 3},
+//    };
+    auto chosen = combos;
+
+    for (const auto &sg_gens : chosen) {
+        const auto s = triangulate<4>(group, sg_gens)
+            .tile(group, g_gens, sg_gens);
+
+        GLuint vao = utilCreateVertexArray();
+        GLuint ibo = utilCreateBuffer();
+        unsigned count = s.size();
+
+        glBindVertexArray(vao);
+        glBindBuffer(GL_ARRAY_BUFFER, ibo);
+        glBufferData(GL_ARRAY_BUFFER, sizeof(Primitive<4>) * s.size(), &s.prims[0], GL_STATIC_DRAW);
+        glEnableVertexAttribArray(0);
+        glVertexAttribIPointer(0, 4, GL_INT, 0, nullptr);
+        glBindBuffer(GL_ARRAY_BUFFER, 0);
+
+        vaos.push_back(vao);
+        ibos.push_back(ibo);
+        counts.push_back(count);
+    }
+    //endregion
 
     GLuint vbo;
     glGenBuffers(1, &vbo);
     glBindBuffer(GL_ARRAY_BUFFER, vbo);
     glBufferData(GL_ARRAY_BUFFER, sizeof(glm::vec4) * points.size(), &points[0], GL_STATIC_DRAW);
+    glBindBuffer(GL_ARRAY_BUFFER, 0);
 
-    glEnableVertexAttribArray(0);
+    GLuint ubo;
-    glVertexAttribPointer(0, 4, GL_FLOAT, false, 0, nullptr);
+    glGenBuffers(1, &ubo);
 
-    GLuint ibo0;
+    glBindBufferBase(GL_SHADER_STORAGE_BUFFER, 1, vbo);
-    glGenBuffers(1, &ibo0);
+    glBindBufferBase(GL_UNIFORM_BUFFER, 1, ubo);
-    glBindBuffer(GL_ARRAY_BUFFER, ibo0);
-    glBufferData(GL_ARRAY_BUFFER, sizeof(int) * g0.size(), &g0[0], GL_STATIC_DRAW);
 
-//    GLuint ibo1;
+    glBindVertexArray(0);
-//    glGenBuffers(1, &ibo1);
-//    glBindBuffer(GL_ARRAY_BUFFER, ibo1);
-//    glBufferData(GL_ARRAY_BUFFER, sizeof(int) * g1.size(), &g1[0], GL_STATIC_DRAW);
-//
-//    GLuint ibo2;
-//    glGenBuffers(1, &ibo2);
-//    glBindBuffer(GL_ARRAY_BUFFER, ibo2);
-//    glBufferData(GL_ARRAY_BUFFER, sizeof(int) * g2.size(), &g2[0], GL_STATIC_DRAW);
-//
-//    GLuint ibo3;
-//    glGenBuffers(1, &ibo3);
-//    glBindBuffer(GL_ARRAY_BUFFER, ibo3);
-//    glBufferData(GL_ARRAY_BUFFER, sizeof(int) * g3.size(), &g3[0], GL_STATIC_DRAW);
 
     while (!glfwWindowShouldClose(window)) {
-        glClear(GL_COLOR_BUFFER_BIT | GL_DEPTH_BUFFER_BIT);
-
-        glUseProgram(pgm);
-        glEnable(GL_PROGRAM_POINT_SIZE);
-        glEnable(GL_POINT_SMOOTH);
-        glEnable(GL_DEPTH_TEST);
-
         int width, height;
         glfwGetFramebufferSize(window, &width, &height);
         glViewport(0, 0, width, height);
-        auto aspect = (float) width / (float) height;
-        auto pheight = 1.4f;
-        auto pwidth = aspect * pheight;
-        glm::mat4 proj = glm::ortho(-pwidth, pwidth, -pheight, pheight);
-        glUniformMatrix4fv(0, 1, false, glm::value_ptr(proj));
 
-        auto t = (float) glfwGetTime() / 3;
+        glClear(GL_COLOR_BUFFER_BIT | GL_DEPTH_BUFFER_BIT);
-        auto view = glm::identity<glm::mat4>();
-        view = glm::rotate(view, t / 1, glm::vec3(0, 1, 0));
-        view = glm::rotate(view, t / 3, glm::vec3(0, 0, 1));
-        view = glm::rotate(view, t / 4, glm::vec3(1, 0, 0));
-        glUniformMatrix4fv(1, 1, false, glm::value_ptr(view));
 
-        glUniform3f(2, 1.0f, 1.0f, 1.0f);
+        auto st = (float) glfwGetTime() / 8;
-        glDrawArrays(GL_POINTS, 0, points.size());
+        Matrices mats = build(window, st);
+        glBindBuffer(GL_UNIFORM_BUFFER, ubo);
+        glBufferData(GL_UNIFORM_BUFFER, sizeof(mats), &mats, GL_STATIC_DRAW);
+        glBindBuffer(GL_UNIFORM_BUFFER, 0);
 
-        glLineWidth(2.0f);
+        for (int i = 0; i < vaos.size(); ++i) {
-        glUniform3f(2, 1.0f, 0.0f, 0.0f);
+            auto c = glm::mix(
-        glBindBuffer(GL_ELEMENT_ARRAY_BUFFER, ibo0);
+                glm::vec3(.3, .2, .5),
-        glDrawElements(GL_LINES, g0.size(), GL_UNSIGNED_INT, nullptr);
+                glm::vec3(.9, .9, .95),
+                (float) (i) / (vaos.size() - 1.f)
+            );
 
-//        glUniform3f(2, 0.0f, 1.0f, 0.0f);
+            glBindProgramPipeline(pipe);
-//        glBindBuffer(GL_ELEMENT_ARRAY_BUFFER, ibo1);
+            glBindVertexArray(vaos[i]);
-//        glDrawElements(GL_LINES, g1.size(), GL_UNSIGNED_INT, nullptr);
+            glProgramUniform3f(fs, 2, c.r, c.g, c.b);
-//
+            glDrawArrays(GL_POINTS, 0, counts[i]);
-//        glUniform3f(2, 0.0f, 0.0f, 1.0f);
+        }
-//        glBindBuffer(GL_ELEMENT_ARRAY_BUFFER, ibo2);
+
-//        glDrawElements(GL_LINES, g2.size(), GL_UNSIGNED_INT, nullptr);
+        glBindProgramPipeline(0);
-//
+        glBindVertexArray(0);
-//        glUniform3f(2, 0.7f, 0.7f, 0.0f);
-//        glBindBuffer(GL_ELEMENT_ARRAY_BUFFER, ibo3);
-//        glDrawElements(GL_LINES, g3.size(), GL_UNSIGNED_INT, nullptr);
 
         glfwSwapBuffers(window);