From 2e4e84ab80ec86a56325df69b34130cd0ffed402 Mon Sep 17 00:00:00 2001
From: allem <allemang.d@gmail.com>
Date: Sat, 14 Sep 2019 19:39:35 -0400
Subject: [PATCH 1/2] correctly generate vertices in N dimensions

---
 cosets/src/main.cpp         | 13 +++----
 cosets/src/mesh.cpp         | 20 ++++++-----
 cosets/src/mirrors.cpp      |  2 +-
 cosets/src/util/coxeter.hpp | 70 +++++++++++++++++++++++++------------
 cosets/src/util/mesh.hpp    | 68 ++++++++++++++---------------------
 cosets/src/util/mirrors.hpp |  2 +-
 6 files changed, 92 insertions(+), 83 deletions(-)

diff --git a/cosets/src/main.cpp b/cosets/src/main.cpp
index 0ac8492..5a46978 100644
--- a/cosets/src/main.cpp
+++ b/cosets/src/main.cpp
@@ -32,12 +32,9 @@ public:
 
       program = build_program("main", vs, fs);
 
-      verts_data = vertices<3>({}, Multiplicites<3>({
-         {0, 1, 5},
-         {1, 2, 3}
-      }), {
-         10, 1, 1
-      });
+      verts_data = vertices<3>(
+         schlafli<3>({5, 3}),
+         {10, 1, 1});
 
       glGenBuffers(1, &verts);
       glBindBuffer(GL_ARRAY_BUFFER, verts);
@@ -81,9 +78,7 @@ public:
 
       glBindVertexArray(tri);
       glUniform4f(0, 1, 1, 1, 1);
-      glDrawArrays(GL_POINTS, 0, verts_data.size() - 3);
-      glUniform4f(0, 1, 1, 1, .1);
-      glDrawArrays(GL_TRIANGLES, verts_data.size() - 3, 3);
+      glDrawArrays(GL_POINTS, 0, verts_data.size());
 
       swapbuffers();
    }
diff --git a/cosets/src/mesh.cpp b/cosets/src/mesh.cpp
index eb26070..4eb26ab 100644
--- a/cosets/src/mesh.cpp
+++ b/cosets/src/mesh.cpp
@@ -1,16 +1,20 @@
 #include "util/mesh.hpp"
 
 int main(int argc, char *argv[]) {
-   auto vs = vertices<3>({}, Multiplicites<3>({
-      {0, 1, 4},
-      {1, 2, 3}
-   }), {
-      1.0f, 1.0f, 1.0f
-   });
+   const auto &mults = schlafli<4>({4, 3, 3});
+   const std::vector<glm::vec4> &vs = vertices<4>(mults, {10, 1, 1, 1});
 
-   for (const auto &v : vs) {
+   for (const auto &v:vs) {
       std::cout << glm::to_string(v) << std::endl;
    }
 
-   return 0;
+//   Table *table = solve({}, mults);
+//   std::cout << table->size() << std::endl;
+//
+//   for (const auto &rel : mults.relations()) {
+//      for (const auto &e : rel) {
+//         std::cout << e << " ";
+//      }
+//      std::cout << std::endl;
+//   }
 }
\ No newline at end of file
diff --git a/cosets/src/mirrors.cpp b/cosets/src/mirrors.cpp
index 51ba1d6..d6281bf 100644
--- a/cosets/src/mirrors.cpp
+++ b/cosets/src/mirrors.cpp
@@ -1,7 +1,7 @@
 #include "util/mirrors.hpp"
 
 int main(int argc, char *argv[]) {
-   auto normals = mirror<3>(Multiplicites<3>({
+   auto normals = mirror<3>(Multiplicities<3>({
       {0, 1, 4},
       {1, 2, 3}
    }));
diff --git a/cosets/src/util/coxeter.hpp b/cosets/src/util/coxeter.hpp
index 054a49f..54d39d6 100644
--- a/cosets/src/util/coxeter.hpp
+++ b/cosets/src/util/coxeter.hpp
@@ -5,39 +5,60 @@
 #include <iomanip>
 #include <algorithm>
 #include <cmath>
-
-struct Mult {
-   int a, b, mult;
-};
+#include <map>
+#include <tuple>
 
 template<int N>
-struct Multiplicites {
-   std::vector<int> mults;
-
-   explicit Multiplicites(const std::vector<Mult> &vals) {
-      mults = std::vector<int>(N * (N - 1) / 2, 2);
-      for (const auto &mult : vals) {
-         set(mult.a, mult.b, mult.mult);
-      }
-   }
+struct Mults {
+   std::map<std::tuple<int, int>, int> mults{};
 
    void set(int a, int b, int mult) {
-      if (a > N or b > N) throw std::logic_error("mirror does not exist");
-      if (a == b) throw std::logic_error("cannot compare mirror to itself");
-      if (a < b) std::swap(a, b); // a is bigger
-
-      mults[a + b] = mult;
+      if (a > b) std::swap(a, b);
+      mults[std::make_tuple(a, b)] = mult;
    }
 
    [[nodiscard]] int get(int a, int b) const {
-      if (a > N or b > N) throw std::logic_error("mirror does not exist");
-      if (a == b) throw std::logic_error("cannot compare mirror to itself");
-      if (a < b) std::swap(a, b); // a is bigger
+      if (a == b) return 1;
+      if (a > b) std::swap(a, b);
 
-      return mults[a + b];
+      const auto &tup = std::make_tuple(a, b);
+      const auto &res = mults.find(tup);
+
+      if (res == mults.end()) return 2;
+      return res->second;
+   }
+
+   [[nodiscard]] std::vector<int> relation(int a, int b) const {
+      std::vector<int> res{};
+      int mult = get(a, b);
+
+      for (int i = 0; i < mult; ++i) {
+         res.push_back(a);
+         res.push_back(b);
+      }
+
+      return res;
+   }
+
+   [[nodiscard]] std::vector<std::vector<int>> relations() const {
+      std::vector<std::vector<int>> res{};
+      for (int a = 0; a < N; ++a) {
+         for (int b = a; b < N; ++b) {
+            res.push_back(relation(a, b));
+         }
+      }
+      return res;
    }
 };
 
+template<int N>
+Mults<N> schlafli(const int (&symbol)[N - 1]) {
+   Mults<N> mults{};
+   for (int i = 0; i < N; ++i) {
+      mults.set(i, i + 1, symbol[i]);
+   }
+   return mults;
+}
 
 struct Table {
    int N;
@@ -171,7 +192,10 @@ std::ostream &operator<<(std::ostream &out, const Row &row) {
    return out;
 }
 
-Table *solve(int gens, const std::vector<int> &subgens, const std::vector<std::vector<int>> &rels) {
+template<int N>
+Table *solve(const std::vector<int> &subgens, const Mults<N> &mults) {
+   const auto rels = mults.relations();
+   int gens = N;
    auto *table = new Table(gens);
 
    for (int gen : subgens)
diff --git a/cosets/src/util/mesh.hpp b/cosets/src/util/mesh.hpp
index 9380bb6..23197d4 100644
--- a/cosets/src/util/mesh.hpp
+++ b/cosets/src/util/mesh.hpp
@@ -7,58 +7,44 @@
 #include <vector>
 
 template<int N>
-std::vector<glm::vec4>
-vertices(const std::vector<int> &subgens, const Multiplicites<N> &mults, const float (&coords)[N]) {
+std::vector<std::vector<int>> coxeter_rels() {
    std::vector<std::vector<int>> rels{};
-   rels.reserve(N + N * (N - 1) / 2);
 
-   for (int i = 0; i < N; ++i) {
-      rels.push_back({i, i});
-   }
-   for (int a = 0; a < N; ++a) {
-      for (int b = 0; b < a; ++b) {
-         int mult = mults.get(a, b);
+   return rels;
+}
 
-         std::vector<int> rel{};
-         for (int i = 0; i < mult; ++i) {
-            rel.push_back(a);
-            rel.push_back(b);
-         }
-         rels.push_back(rel);
-      }
-   }
-
-   Table *table = solve(N, subgens, rels);
-
-//   std::cout << table->size() << std::endl;
-//   std::cout << *table << std::endl;
-
-//   for (const auto &v : table->words()) {
-//      std::cout << "[ ";
-//      for (auto e : v) {
-//         std::cout << e << " ";
-//      }
-//      std::cout << "]\n";
-//   }
-
-   const std::vector<glm::vec4> &normals = mirror<3>(mults);
-   const std::vector<glm::vec4> &corners = plane_intersections(normals);
+template<int N>
+glm::vec4 identity(const std::vector<glm::vec4> &normals, const float(&coords)[N]) {
+   const std::vector<glm::vec4> corners = plane_intersections(normals);
    const std::vector<float> coords_vec(coords, coords + N);
-   const glm::vec4 &identity = glm::normalize(barycentric(corners, coords_vec));
+   const glm::vec4 identity = barycentric(corners, coords_vec);
+   return glm::normalize(identity);
+}
+
+template<int N>
+std::vector<glm::vec4>
+vertices(const Mults<N> &mults, const float (&coords)[N]) {
+   Table *table = solve({}, mults);
+
+   const std::vector<glm::vec4> normals = mirror<N>(mults);
+   glm::vec4 ident = identity(normals, coords);
 
    std::vector<glm::vec4> verts{};
-
    for (const auto &word : table->words()) {
-      glm::vec4 vert = identity;
+      glm::vec4 vert = ident;
       for (const auto &gen : word) {
          vert = reflect(vert, normals[gen]);
       }
       verts.push_back(vert);
    }
 
-   for (const auto &e : corners) {
-      verts.push_back(e * 1.1f);
-   }
-
    return verts;
-}
\ No newline at end of file
+}
+
+template<int N>
+std::vector<unsigned> lines(const Mults<N> &mults) {
+//   Table *table = solve(N, {}, coxeter_rels(mults));
+
+   return {};
+}
+
diff --git a/cosets/src/util/mirrors.hpp b/cosets/src/util/mirrors.hpp
index 6a629bd..ce107bd 100644
--- a/cosets/src/util/mirrors.hpp
+++ b/cosets/src/util/mirrors.hpp
@@ -11,7 +11,7 @@
 #include "coxeter.hpp"
 
 template<int N>
-std::vector<glm::vec4> mirror(const Multiplicites<N> &mults) {
+std::vector<glm::vec4> mirror(const Mults<N> &mults) {
    static_assert(1 <= N and N <= 4, "Vector size is unsupported");
 
    std::vector<glm::vec4> mirrors{};

From a3573acca4da426863db01bde7d53e57ae6bcbca Mon Sep 17 00:00:00 2001
From: allem <allemang.d@gmail.com>
Date: Sat, 14 Sep 2019 21:42:12 -0400
Subject: [PATCH 2/2] correctly generate and render edges

---
 cosets/src/main.cpp             | 48 ++++++++++++------
 cosets/src/mesh.cpp             | 24 ++++-----
 cosets/src/shaders/main.vs.glsl |  6 ++-
 cosets/src/util/coxeter.hpp     | 87 ++++++++++++++++++++++++---------
 cosets/src/util/mesh.hpp        | 28 ++++++++---
 cosets/src/util/mirrors.hpp     |  2 +-
 cosets/src/util/numeric.hpp     | 20 ++++++++
 7 files changed, 156 insertions(+), 59 deletions(-)

diff --git a/cosets/src/main.cpp b/cosets/src/main.cpp
index 5a46978..15b8336 100644
--- a/cosets/src/main.cpp
+++ b/cosets/src/main.cpp
@@ -13,11 +13,14 @@ using namespace std;
 
 class CosetsWindow : public Window {
    GLint program;
-   GLuint tri;
+   GLuint points_vao, edges_vao;
 
-   GLuint verts;
+   GLuint verts_buf;
    std::vector<glm::vec4> verts_data;
 
+   GLuint edges_buf;
+   std::vector<int> edge_data;
+
 public:
    void init() override {
       auto vs = build_shader_file(
@@ -32,21 +35,34 @@ public:
 
       program = build_program("main", vs, fs);
 
-      verts_data = vertices<3>(
-         schlafli<3>({5, 3}),
-         {10, 1, 1});
+      const Mults &mults = schlafli<4>({3, 4, 3});
+      verts_data = vertices<4>(mults, {10, 1, 1, 1});
+      edge_data = edges<4>(mults);
 
-      glGenBuffers(1, &verts);
-      glBindBuffer(GL_ARRAY_BUFFER, verts);
+      glGenBuffers(1, &verts_buf);
+      glBindBuffer(GL_ARRAY_BUFFER, verts_buf);
       glBufferData(GL_ARRAY_BUFFER, sizeof(glm::vec4) * verts_data.size(), &verts_data[0], GL_STATIC_DRAW);
 
+      glGenBuffers(1, &edges_buf);
+      glBindBuffer(GL_ARRAY_BUFFER, edges_buf);
+      glBufferData(GL_ARRAY_BUFFER, sizeof(int) * edge_data.size(), &edge_data[0], GL_STATIC_DRAW);
 
-      glGenVertexArrays(1, &tri);
-      glBindVertexArray(tri);
-      glBindBuffer(GL_ARRAY_BUFFER, verts);
+      glGenVertexArrays(1, &points_vao);
+      glBindVertexArray(points_vao);
+      glBindBuffer(GL_ARRAY_BUFFER, verts_buf);
       glEnableVertexAttribArray(0);
       glVertexAttribPointer(0, 4, GL_FLOAT, false, 0, nullptr);
       glBindBuffer(GL_ARRAY_BUFFER, 0);
+
+      glGenVertexArrays(1, &edges_vao);
+      glBindVertexArray(edges_vao);
+      glBindBuffer(GL_ARRAY_BUFFER, verts_buf);
+      glEnableVertexAttribArray(0);
+      glVertexAttribPointer(0, 4, GL_FLOAT, false, 0, nullptr);
+      glBindBuffer(GL_ARRAY_BUFFER, 0);
+      glBindBuffer(GL_ELEMENT_ARRAY_BUFFER, edges_buf);
+
+      glBindVertexArray(0);
    }
 
    void render() override {
@@ -66,26 +82,30 @@ public:
 
       const glm::mat4 proj = glm::ortho(-ar * sc, ar * sc, -sc, sc, -10.f, 10.f);
       const glm::mat4 view = glm::rotate(id, angle, ax_1) * glm::rotate(id, angle, ax_2);
-      const glm::mat4 mat = proj * view;
 
       glUseProgram(program);
-      glUniformMatrix4fv(1, 1, false, glm::value_ptr(mat));
+      glUniformMatrix4fv(1, 1, false, glm::value_ptr(proj));
+      glUniformMatrix4fv(2, 1, false, glm::value_ptr(view));
       glEnable(GL_DEPTH_TEST);
       glEnable(GL_POINT_SMOOTH);
       glEnable(GL_BLEND);
       glBlendFunc(GL_SRC_ALPHA, GL_ONE_MINUS_SRC_ALPHA);
       glPointSize(10.f);
 
-      glBindVertexArray(tri);
+      glBindVertexArray(points_vao);
       glUniform4f(0, 1, 1, 1, 1);
       glDrawArrays(GL_POINTS, 0, verts_data.size());
 
+      glBindVertexArray(edges_vao);
+      glUniform4f(0, 1, 1, 0, 1);
+      glDrawElements(GL_LINES, edge_data.size(), GL_UNSIGNED_INT, 0);
+
       swapbuffers();
    }
 
    void deinit() override {
       glDeleteProgram(program);
-      glDeleteVertexArrays(1, &tri);
+      glDeleteVertexArrays(1, &points_vao);
    }
 };
 
diff --git a/cosets/src/mesh.cpp b/cosets/src/mesh.cpp
index 4eb26ab..2ef34ed 100644
--- a/cosets/src/mesh.cpp
+++ b/cosets/src/mesh.cpp
@@ -1,20 +1,14 @@
 #include "util/mesh.hpp"
 
 int main(int argc, char *argv[]) {
-   const auto &mults = schlafli<4>({4, 3, 3});
-   const std::vector<glm::vec4> &vs = vertices<4>(mults, {10, 1, 1, 1});
+   const int N = 3;
+   const Mults &mults = schlafli<N>({4, 3});
+   const std::vector<glm::vec4> &vs = vertices<N>(mults, {10, 1, 1});
 
-   for (const auto &v:vs) {
-      std::cout << glm::to_string(v) << std::endl;
-   }
+   std::cout << "# verts: " << vs.size() << std::endl;
 
-//   Table *table = solve({}, mults);
-//   std::cout << table->size() << std::endl;
-//
-//   for (const auto &rel : mults.relations()) {
-//      for (const auto &e : rel) {
-//         std::cout << e << " ";
-//      }
-//      std::cout << std::endl;
-//   }
-}
\ No newline at end of file
+   const auto es = edges<N>(mults);
+
+
+   return 0;
+}
diff --git a/cosets/src/shaders/main.vs.glsl b/cosets/src/shaders/main.vs.glsl
index a9f2ad3..a9269cc 100644
--- a/cosets/src/shaders/main.vs.glsl
+++ b/cosets/src/shaders/main.vs.glsl
@@ -1,6 +1,7 @@
 #version 440 core
 
 uniform mat4 proj;
+uniform mat4 view;
 
 in vec4 pos;
 out vec4 v;
@@ -9,7 +10,10 @@ out vec4 screen;
 void main(){
   v = pos;
 
+  vec4 vert = view * pos;
+
   /* stereographic projection */
-  vec4 vert = vec4(pos.xyz / (1 - pos.w), 1);
+  vert = vec4(vert.xyw / (1 - vert.z), 1);
+
   screen = gl_Position = proj * vert;
 }
diff --git a/cosets/src/util/coxeter.hpp b/cosets/src/util/coxeter.hpp
index 54d39d6..f385992 100644
--- a/cosets/src/util/coxeter.hpp
+++ b/cosets/src/util/coxeter.hpp
@@ -1,6 +1,7 @@
 #pragma once
 
 #include <iostream>
+#include <utility>
 #include <vector>
 #include <iomanip>
 #include <algorithm>
@@ -8,7 +9,6 @@
 #include <map>
 #include <tuple>
 
-template<int N>
 struct Mults {
    std::map<std::tuple<int, int>, int> mults{};
 
@@ -40,11 +40,11 @@ struct Mults {
       return res;
    }
 
-   [[nodiscard]] std::vector<std::vector<int>> relations() const {
+   [[nodiscard]] std::vector<std::vector<int>> relations(std::vector<int> gens) const {
       std::vector<std::vector<int>> res{};
-      for (int a = 0; a < N; ++a) {
-         for (int b = a; b < N; ++b) {
-            res.push_back(relation(a, b));
+      for (int a = 0; a < gens.size(); ++a) {
+         for (int b = a; b < gens.size(); ++b) {
+            res.push_back(relation(gens[a], gens[b]));
          }
       }
       return res;
@@ -52,8 +52,8 @@ struct Mults {
 };
 
 template<int N>
-Mults<N> schlafli(const int (&symbol)[N - 1]) {
-   Mults<N> mults{};
+Mults schlafli(const int (&symbol)[N - 1]) {
+   Mults mults{};
    for (int i = 0; i < N; ++i) {
       mults.set(i, i + 1, symbol[i]);
    }
@@ -61,23 +61,31 @@ Mults<N> schlafli(const int (&symbol)[N - 1]) {
 }
 
 struct Table {
-   int N;
+   const std::vector<int> gens;
    std::vector<std::vector<int>> fwd{};
    std::vector<std::vector<int>> rev{};
 
-   explicit Table(int N) {
-      this->N = N;
+   explicit Table(const std::vector<int> gens) : gens(gens) {
       add_row();
    }
 
+   [[nodiscard]] int gen_index(int gen) const {
+      for (unsigned i = 0; i < gens.size(); ++i)
+         if (gens[i] == gen)
+            return i;
+
+      throw std::logic_error("Referencing nonexistant generator");
+   }
+
    void add_row() {
+      int N = gens.size();
       fwd.emplace_back(N, -1);
       rev.emplace_back(N, -1);
    }
 
    int add_coset() {
       for (int from = 0; from < (int) size(); ++from) {
-         for (int gen = 0; gen < N; ++gen) {
+         for (int gen : gens) {
             if (get(from, gen) < 0) {
                int to = (int) size();
                add_row();
@@ -94,16 +102,19 @@ struct Table {
    }
 
    void set(int from, int gen, int to) {
-      fwd[from][gen] = to;
-      rev[to][gen] = from;
+      int i = gen_index(gen);
+      fwd[from][i] = to;
+      rev[to][i] = from;
    }
 
    [[nodiscard]] int get(int from, int gen) const {
-      return fwd[from][gen];
+      int i = gen_index(gen);
+      return fwd[from][i];
    }
 
    [[nodiscard]] int rget(int gen, int to) const {
-      return rev[to][gen];
+      int i = gen_index(gen);
+      return rev[to][i];
    }
 
    std::vector<std::vector<int>> words() {
@@ -117,7 +128,7 @@ struct Table {
                continue;
             }
 
-            for (int gen = 0; gen < (int) N; ++gen) {
+            for (int gen : gens) {
                int to = get(from, gen);
                if (vecs[to] != nullptr) {
                   continue;
@@ -137,6 +148,31 @@ struct Table {
 
       return res;
    }
+
+   [[nodiscard]] int apply(int e, const std::vector<int> &word) const {
+      for (const auto &gen : word) {
+         e = get(e, gen);
+      }
+      return e;
+   }
+
+   [[nodiscard]] int apply(const std::vector<int> &word) const {
+      return apply(0, word);
+   }
+
+   [[nodiscard]] std::vector<int> apply_each(int e, const std::vector<std::vector<int>> &words) const {
+      std::vector<int> res{};
+
+      for (const auto &word:words) {
+         res.push_back(apply(e, word));
+      }
+
+      return res;
+   }
+
+   [[nodiscard]] std::vector<int> apply_each(const std::vector<std::vector<int>> &words) const {
+      return apply_each(0, words);
+   }
 };
 
 struct Row {
@@ -192,10 +228,8 @@ std::ostream &operator<<(std::ostream &out, const Row &row) {
    return out;
 }
 
-template<int N>
-Table *solve(const std::vector<int> &subgens, const Mults<N> &mults) {
-   const auto rels = mults.relations();
-   int gens = N;
+Table *solve(const std::vector<int> &gens, const std::vector<int> &subgens, const Mults &mults) {
+   const auto rels = mults.relations(gens);
    auto *table = new Table(gens);
 
    for (int gen : subgens)
@@ -231,6 +265,15 @@ Table *solve(const std::vector<int> &subgens, const Mults<N> &mults) {
    return table;
 }
 
+template<int N>
+Table *solve(const std::vector<int> &subgens, const Mults &mults) {
+   std::vector<int> gens{};
+   for (int i = 0; i < N; ++i) {
+      gens.push_back(i);
+   }
+   return solve(gens, subgens, mults);
+}
+
 std::ostream &operator<<(std::ostream &out, const Table &table) {
    int k = ceil(log10(table.size()));
 
@@ -238,10 +281,10 @@ std::ostream &operator<<(std::ostream &out, const Table &table) {
    for (unsigned j = 0; j < table.size(); ++j) {
       auto arr = table.fwd[j];
       out << " " << std::setw(k) << j << " [";
-      for (int i = 0; i < table.N; ++i) {
+      for (int i = 0; i < (int) table.gens.size(); ++i) {
          out << arr[i];
 
-         if (i < table.N - 1)
+         if (i < table.gens.size() - 1)
             out << " ";
       }
 
diff --git a/cosets/src/util/mesh.hpp b/cosets/src/util/mesh.hpp
index 23197d4..01a666b 100644
--- a/cosets/src/util/mesh.hpp
+++ b/cosets/src/util/mesh.hpp
@@ -23,8 +23,8 @@ glm::vec4 identity(const std::vector<glm::vec4> &normals, const float(&coords)[N
 
 template<int N>
 std::vector<glm::vec4>
-vertices(const Mults<N> &mults, const float (&coords)[N]) {
-   Table *table = solve({}, mults);
+vertices(const Mults &mults, const float (&coords)[N]) {
+   Table *table = solve<N>({}, mults);
 
    const std::vector<glm::vec4> normals = mirror<N>(mults);
    glm::vec4 ident = identity(normals, coords);
@@ -42,9 +42,25 @@ vertices(const Mults<N> &mults, const float (&coords)[N]) {
 }
 
 template<int N>
-std::vector<unsigned> lines(const Mults<N> &mults) {
-//   Table *table = solve(N, {}, coxeter_rels(mults));
+std::vector<int> edges(const Mults &mults) {
+   std::vector<int> res{};
 
-   return {};
+   Table *verts = solve<N>({}, mults);
+
+   int K = 1;
+   for (const auto &subgens : combinations(N, K)) {
+      Table *edge = solve(subgens, {}, mults);
+
+      std::vector<int> primitive = verts->apply_each(edge->words());
+
+      Table *cosets = solve<N>(subgens, mults);
+
+      for (const auto &coset : cosets->words()) {
+         for (const auto &e : primitive) {
+            res.push_back(verts->apply(e, coset));
+         }
+      }
+   }
+
+   return res;
 }
-
diff --git a/cosets/src/util/mirrors.hpp b/cosets/src/util/mirrors.hpp
index ce107bd..9551a26 100644
--- a/cosets/src/util/mirrors.hpp
+++ b/cosets/src/util/mirrors.hpp
@@ -11,7 +11,7 @@
 #include "coxeter.hpp"
 
 template<int N>
-std::vector<glm::vec4> mirror(const Mults<N> &mults) {
+std::vector<glm::vec4> mirror(const Mults &mults) {
    static_assert(1 <= N and N <= 4, "Vector size is unsupported");
 
    std::vector<glm::vec4> mirrors{};
diff --git a/cosets/src/util/numeric.hpp b/cosets/src/util/numeric.hpp
index 663be4a..064892c 100644
--- a/cosets/src/util/numeric.hpp
+++ b/cosets/src/util/numeric.hpp
@@ -57,3 +57,23 @@ std::vector<glm::vec4> plane_intersections(std::vector<glm::vec4> normals) {
 
    return results;
 }
+
+
+std::vector<std::vector<int>> combinations(int N, int K) {
+   std::string bitmask(K, 1); // K leading 1's
+   bitmask.resize(N, 0); // N-K trailing 0's
+   std::vector<std::vector<int>> combos{};
+
+   do {
+      std::vector<int> combo{};
+      for (int i = 0; i < N; ++i) // [0..N-1] integers
+      {
+         if (bitmask[i]) {
+            combo.push_back(i);
+         }
+      }
+      combos.push_back(combo);
+   } while (std::prev_permutation(bitmask.begin(), bitmask.end()));
+
+   return combos;
+}