point generation working; still needs broken into functions though.

tc/include/tc/cosets.h

@@ -26,5 +26,7 @@ namespace tc {
         [[nodiscard]] int get(int coset, int gen) const;
 
         [[nodiscard]] int get(int idx) const;
+
+        [[nodiscard]] int size() const;
     };
 }

tc/include/tc/groups.h

@@ -21,6 +21,8 @@ namespace tc {
 
         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;

tc/src/cosets.cpp

@@ -38,4 +38,7 @@ namespace tc {
         return data[idx];
     }
 
+    int Cosets::size() const {
+        return len;
+    }
 }

tc/src/groups.cpp

@@ -11,20 +11,19 @@ namespace tc {
         }
 
         for (Rel rel : rels) {
-            if (rel.gens[0] < rel.gens[1])
-                _mults[rel.gens[0]][rel.gens[1]] = rel.mult;
-            else
-                _mults[rel.gens[1]][rel.gens[0]] = rel.mult;
+            setmult(rel);
         }
     }
 
     void Group::setmult(Rel rel) {
-        if (rel.gens[0] < rel.gens[1])
         _mults[rel.gens[0]][rel.gens[1]] = rel.mult;
-        else
         _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++) {
@@ -75,7 +74,7 @@ namespace tc {
         return g.power(p);
     }
 
-    Group schlafli(const std::vector<int> &mults, const std::string& name) {
+    Group schlafli(const std::vector<int> &mults, const std::string &name) {
         int ngens = (int) mults.size() + 1;
 
         Group g(ngens, {}, name);

vis/include/geom.h

@@ -0,0 +1,88 @@
+#pragma once
+
+#include <tc/cosets.h>
+#include <glm/glm.hpp>
+#include <cmath>
+#include <vector>
+#include <algorithm>
+
+glm::vec4 round(const glm::vec4 &f, int prec) {
+    auto dec = (float) pow(10, prec);
+    auto res = glm::trunc(f * dec + 0.5f) / dec;
+    return res;
+}
+
+float dot(int n, const glm::vec4 &a, const glm::vec4 &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;
+
+    for (int p = 0; p < group.ngens; ++p) {
+        glm::vec4 vp{};
+        for (int m = 0; m < p; ++m) {
+            glm::vec4 vq = mirrors[m];
+            std::cout << "p " << p << " m " << m << " - " << group.rel(p, m).mult << "; ";
+            vp[m] = (cos(M_PI / group.rel(p, m).mult) - dot(m, vp, vq)) / vq[m];
+        }
+        std::cout << std::endl;
+        vp[p] = std::sqrt(1 - glm::dot(vp, vp));
+
+        for (const auto &v : mirrors) {
+            if (glm::dot(vp, v) > 0) {
+                vp *= -1;
+                break;
+            }
+        }
+
+        mirrors.push_back(round(vp, 15));
+    }
+
+    return mirrors;
+}
+
+glm::vec4 project(const glm::vec4 &vec, const glm::vec4 &target) {
+    return glm::dot(vec, target) / glm::dot(target, target) * target;
+}
+
+glm::vec4 reflect(const glm::vec4 &vec, const glm::vec4 axis) {
+    return vec - 2.f * project(vec, axis);
+}
+
+glm::vec4 gram_schmidt_last(std::vector<glm::vec4> vecs) {
+    int N = vecs.size();
+    for (int i = 0; i < N; ++i) {
+        for (int j = 0; j < i; ++j) {
+            vecs[i] -= project(vecs[i], vecs[j]);
+        }
+    }
+
+    return glm::normalize(vecs[N - 1]);
+}
+
+glm::vec4 barycentric(std::vector<glm::vec4> basis, std::vector<float> coords) {
+    glm::vec4 res{};
+
+    int N = std::min(basis.size(), coords.size());
+    for (int i = 0; i < N; ++i) {
+        res += basis[i] * coords[i];
+    }
+    return glm::normalize(res);
+}
+
+std::vector<glm::vec4> plane_intersections(std::vector<glm::vec4> normals) {
+    int N = normals.size();
+    std::vector<glm::vec4> results(N);
+
+    for (int i = 0; i < N; ++i) {
+        std::rotate(normals.begin(), normals.begin() + 1, normals.end());
+        results[i] = gram_schmidt_last(normals);
+    }
+
+    return results;
+}

vis/src/main.cpp

@@ -5,6 +5,10 @@
 #include <tc/groups.h>
 #include <tc/solver.h>
 
+#include "geom.h"
+
+#include <glm/gtx/string_cast.hpp>
+
 #ifdef _WIN32
 extern "C" {
 __attribute__((unused)) __declspec(dllexport) int NvOptimusEnablement = 0x00000001;
@@ -39,13 +43,34 @@ int main(int argc, char *argv[]) {
         << "    OpenGL version:  " << glGetString(GL_VERSION) << std::endl
         << "    Shading version: " << glGetString(GL_SHADING_LANGUAGE_VERSION) << std::endl;
 
-    auto group = tc::group::A(5);
+    auto group = tc::group::B(4);
     auto res = tc::solve(group);
     std::cout
         << "Coset Solution Test:" << std::endl
         << "    Group: " << group.name << std::endl
         << "    Order: " << res.len << std::endl;
 
+    auto mirrors = mirror(group);
+    std::cout << "Mirrors:" << std::endl;
+    for (const auto &mirror : mirrors) {
+        std::cout << "    " << glm::to_string(mirror) << std::endl;
+    }
+
+    auto corners = plane_intersections(mirrors);
+
+    auto points = std::vector<glm::vec4>(res.size());
+    points[0] = barycentric(corners, {1, 0, 0, 0});
+
+    for (int i = 1; i < res.size(); i++) {
+        auto action = res.path[i];
+        points[i] = reflect(points[action.coset], mirrors[action.gen]);
+    }
+
+    std::cout << "Points:" << std::endl;
+    for (const auto &point : points) {
+        std::cout << "    " << to_string(point) << std::endl;
+    }
+
     while (!glfwWindowShouldClose(window)) {
         glClear(GL_COLOR_BUFFER_BIT);
         glfwSwapBuffers(window);