diff --git a/example/bench.cpp b/example/bench.cpp
index 6744ee0..5f993ef 100644
--- a/example/bench.cpp
+++ b/example/bench.cpp
@@ -15,24 +15,24 @@ void test(const G &group) {
     int order = cosets.size();
 
     std::cout
-        << std::setw(7) << group.name << ", "
+        << std::setw(7) << group.name() << ", "
         << std::setw(7) << order << ", "
         << std::fixed << std::setprecision(6) << diff << "s"
         << std::endl;
 }
 
 int main() {
-    test(tc::group::H(2));
-    test(tc::group::H(3));
-    test(tc::group::H(4));
+    test(tc::group::H<2>());
+    test(tc::group::H<3>());
+    test(tc::group::H<4>());
     test(tc::group::T(100));
     test(tc::group::T(500));
     test(tc::group::T(1000));
-    test(tc::group::E(6));
-    test(tc::group::E(7));
-    test(tc::group::B(6));
-    test(tc::group::B(7));
-    test(tc::group::B(8));
+    test(tc::group::E<6>());
+    test(tc::group::E<7>());
+    test(tc::group::B<6>());
+    test(tc::group::B<7>());
+    test(tc::group::B<8>());
 
     return 0;
 }
diff --git a/include/tc/core.hpp b/include/tc/core.hpp
index 73b3a7a..25f2fb9 100644
--- a/include/tc/core.hpp
+++ b/include/tc/core.hpp
@@ -2,7 +2,6 @@
 
 #include "tc/core/action.hpp"
 #include "tc/core/path.hpp"
-#include "tc/core/cosets.hpp"
 #include "tc/core/rel.hpp"
 #include "tc/core/group.hpp"
 #include "tc/core/solver.hpp"
diff --git a/include/tc/core/cosets.hpp b/include/tc/core/cosets.hpp
deleted file mode 100644
index d028019..0000000
--- a/include/tc/core/cosets.hpp
+++ /dev/null
@@ -1,39 +0,0 @@
-#pragma once
-
-#include <vector>
-
-#include "path.hpp"
-
-struct Cosets {
-    int ngens;
-    std::vector<int> data;
-    Path path;
-
-    Cosets(const Cosets &) = default;
-
-    explicit Cosets(int ngens)
-        : ngens(ngens) {
-    }
-
-    void add_row() {
-        data.resize(data.size() + ngens, -1);
-        path.add_row();
-    }
-
-    void put(int coset, int gen, int target) {
-        data[coset * ngens + gen] = target;
-        data[target * ngens + gen] = coset;
-
-        if (path.get(target).from_idx == -1) {
-            path.put(coset, gen, target);
-        }
-    }
-
-    [[nodiscard]] int get(int coset, int gen) const {
-        return data[coset * ngens + gen];
-    }
-
-    [[nodiscard]] size_t size() const {
-        return path.size();
-    }
-};
diff --git a/include/tc/core/group.hpp b/include/tc/core/group.hpp
index e89bd4a..cc8629c 100644
--- a/include/tc/core/group.hpp
+++ b/include/tc/core/group.hpp
@@ -8,126 +8,112 @@
 #include <queue>
 
 #include "rel.hpp"
-#include "cosets.hpp"
 
 #include <Eigen/Eigen>
 #include <iostream>
 
 namespace tc {
+    template<unsigned int Rank>
     struct Group;
+
+    template<unsigned int Rank, unsigned int PRank>
     struct SubGroup;
 
-    struct Group {
-        using Matrix = Eigen::MatrixXi;
+    template<unsigned int Rank>
+    class Group {
+    public:
+        using Matrix = Eigen::Matrix<unsigned int, Rank, Rank>;
 
-        int ngens;
-        std::string name;
+    private:
+    public:
+        std::string _name;
         Matrix _data;
         Eigen::SelfAdjointView<Matrix, Eigen::Upper> _mults;
 
-        Group(const Group &g)
-            : ngens(g.ngens),
-            name(g.name),
+    public:
+        Group(const Group<Rank> &g) :
+            _name(g._name),
             _data(g._data),
             _mults(_data) {
         }
 
-        Group(Group &&g) noexcept
-            : ngens(g.ngens),
-            name(std::move(g.name)),
+        Group(Group &&g) noexcept:
+            _name(std::move(g._name)),
             _data(std::move(g._data)),
             _mults(_data) {
         }
 
-        explicit Group(
-            int ngens,
-            std::string name = "G"
-        ) : ngens(ngens),
-            name(std::move(name)),
-            _data(ngens, ngens),
+        explicit Group(std::string name = "G") :
+            _name(std::move(name)),
+            _data(),
             _mults(_data) {
             _data.fill(2);
         }
 
-        Matrix::Scalar &operator()(int a, int b) {
+        unsigned int rank() const {
+            return _data.rows();
+        }
+
+        std::string name() const {
+            return _name;
+        }
+
+        typename Matrix::Scalar &operator()(int a, int b) {
             return _mults(a, b);
         }
 
-        Matrix::Scalar operator()(int a, int b) const {
+        typename Matrix::Scalar operator()(int a, int b) const {
             return _mults(a, b);
         }
 
         [[nodiscard]] std::vector<Rel> get_rels() const {
             std::vector<Rel> res;
-            for (int i = 0; i < ngens - 1; ++i) {
-                for (int j = i + 1; j < ngens; ++j) {
-                    res.emplace_back(i, j, operator()(i, j));
+            for (int i = 0; i < Rank - 1; ++i) {
+                for (int j = i + 1; j < Rank; ++j) {
+                    res.emplace_back(i, j, _mults(i, j));
                 }
             }
             return res;
         }
-
-        [[nodiscard]] SubGroup subgroup(
-            const std::vector<int> &gens
-        ) const;
     };
 
-    struct SubGroup : public Group {
-        std::vector<int> gen_map;
-        const Group &parent;
-
-        SubGroup(const Group &parent, std::vector<int> gen_map) : Group(gen_map.size()), parent(parent) {
-
-            std::sort(gen_map.begin(), gen_map.end());
-            this->gen_map = gen_map;
-
-            for (size_t i = 0; i < gen_map.size(); ++i) {
-                for (size_t j = 0; j < gen_map.size(); ++j) {
-                    int mult = parent(gen_map[i], gen_map[j]);
-                    operator()(i, j) = mult;
-                }
-            }
-        }
-    };
-
-    SubGroup Group::subgroup(const std::vector<int> &gens) const {
-        return SubGroup(*this, gens);
-    }
-
-    Group product(const Group &g, const Group &h) {
+    template<unsigned int GR, unsigned int HR>
+    Group<GR + HR> product(const Group<GR> &g, const Group<HR> &h) {
         std::stringstream ss;
         ss << g.name << "*" << h.name;
 
-        Group res(g.ngens + h.ngens, ss.str());
+        Group<GR + HR> res(ss.str());
 
         int off = 0;
-
-        for (int i = 0; i < g.ngens; ++i) {
-            for (int j = i; j < g.ngens; ++j) {
+        for (int i = 0; i < GR; ++i) {
+            for (int j = i; j < GR; ++j) {
                 res(i + off, j + off) = g(i, j);
             }
         }
-        off += g.ngens;
+        off += GR;
 
-        for (int i = 0; i < h.ngens; ++i) {
-            for (int j = i; j < h.ngens; ++j) {
+        for (int i = 0; i < HR; ++i) {
+            for (int j = i; j < HR; ++j) {
                 res(i + off, j + off) = h(i, j);
             }
         }
+        off += HR;
 
         return res;
     }
 
-    Group power(const Group &g, int p) {
+    template<unsigned int GR, unsigned int P>
+    Group<GR * P> power(const Group<GR> &g) {
         std::stringstream ss;
-        ss << g.name << "^" << p;
+        ss << g.name << "^" << P;
 
-        Group res(g.ngens * p, ss.str());
+        Group<GR * P> res(ss.str());
 
-        for (int i = 0; i < g.ngens; ++i) {
-            for (int j = i; j < g.ngens; ++j) {
-                for (int k = 0; k < p; ++k) {
-                    int off = k * g.ngens;
+        for (int k = 0; k < P; ++k) {
+            int off = k * GR;
+
+            for (int i = 0; i < GR; ++i) {
+                for (int j = i; j < GR; ++j) {
                     res(i + off, j + off) = g(i, j);
                 }
             }
@@ -136,11 +122,11 @@ namespace tc {
         return res;
     }
 
-    Group operator*(const Group &g, const Group &h) {
-        return product(g, h);
-    }
-
-    Group operator^(const Group &g, int p) {
-        return power(g, p);
-    }
+//    Group operator*(const Group &g, const Group &h) {
+//        return product(g, h);
+//    }
+//
+//    Group operator^(const Group &g, int p) {
+//        return power(g, p);
+//    }
 }
diff --git a/include/tc/core/solver.hpp b/include/tc/core/solver.hpp
index 5470ae1..736792e 100644
--- a/include/tc/core/solver.hpp
+++ b/include/tc/core/solver.hpp
@@ -3,6 +3,42 @@
 #include <memory>
 #include <vector>
 
+namespace tc {
+    struct Cosets {
+        int ngens;
+        std::vector<int> data;
+        Path path;
+
+        Cosets(const Cosets &) = default;
+
+        explicit Cosets(int ngens)
+            : ngens(ngens) {
+        }
+
+        void add_row() {
+            data.resize(data.size() + ngens, -1);
+            path.add_row();
+        }
+
+        void put(int coset, int gen, int target) {
+            data[coset * ngens + gen] = target;
+            data[target * ngens + gen] = coset;
+
+            if (path.get(target).from_idx == -1) {
+                path.put(coset, gen, target);
+            }
+        }
+
+        [[nodiscard]] int get(int coset, int gen) const {
+            return data[coset * ngens + gen];
+        }
+
+        [[nodiscard]] size_t size() const {
+            return path.size();
+        }
+    };
+}
+
 namespace {
     struct Row {
         int gnr;
@@ -20,20 +56,24 @@ namespace {
         }
     };
 
-    struct Tables {
+    template<unsigned int Rank>
+    class Tables {
+    public:
+        static constexpr unsigned int Rels = Rank * (Rank + 1) / 2 - Rank;
+
     private:
         std::shared_ptr<int> null_lst_ptr = std::make_shared<int>();
-        int ngens;
 
-        std::vector<std::shared_ptr<Table>> tables;
-        std::vector<std::vector<std::shared_ptr<Table>>> deps;
+        std::array<std::shared_ptr<Table>, Rels> tables;
+        std::array<std::vector<std::shared_ptr<Table>>, Rank> deps;
 
     public:
-        explicit Tables(const tc::Group &group)
-            : ngens(group.ngens), deps(ngens) {
-            for (const auto &rel: group.get_rels()) {
+        explicit Tables(const tc::Group<Rank> &group) {
+            const auto &rels = group.get_rels();
+            for (int i = 0; i < Rels; ++i) {
+                const auto &rel = rels[i];
                 auto table = std::make_shared<Table>(rel);
-                tables.push_back(table);
+                tables[i] = table;
                 deps[rel.gens[0]].push_back(table);
                 deps[rel.gens[1]].push_back(table);
             }
@@ -46,7 +86,7 @@ namespace {
             }
         }
 
-        void initialize(int target, const Cosets &cosets) {
+        void initialize(int target, const tc::Cosets &cosets) {
             for (auto &table: tables) {
                 const Rel &rel = table->rel;
                 Row &row = table->rows[target];
@@ -64,7 +104,7 @@ namespace {
             }
         }
 
-        void learn(int coset, int gen, int target, const Cosets &cosets, std::priority_queue<int> &facts) {
+        void learn(int coset, int gen, int target, const tc::Cosets &cosets, std::priority_queue<int> &facts) {
             if (target == coset) {
                 for (auto &table: deps[gen]) {
                     Row &target_row = table->rows[target];
@@ -92,11 +132,11 @@ namespace {
                         // forward learn
                         int lst = *target_row.lst;
                         int gen_ = rel.gens[rel.gens[0] == gen];
-                        facts.push(lst * ngens + gen_);
+                        facts.push(lst * Rank + gen_);
                     } else if (target_row.gnr == -rel.mult) {
                         // stationary learn
                         int gen_ = rel.gens[rel.gens[0] == gen];
-                        facts.push(target * ngens + gen_);
+                        facts.push(target * Rank + gen_);
                     } else if (target_row.gnr == rel.mult - 1) {
                         // determined family
                         *target_row.lst = target;
@@ -108,28 +148,29 @@ namespace {
 }
 
 namespace tc {
-    Cosets solve(const Group &g, const std::vector<int> &sub_gens = {}) {
-        Cosets cosets(g.ngens);
+    template<unsigned int Rank>
+    tc::Cosets solve(const Group <Rank> &g, const std::vector<int> &sub_gens = {}) {
+        tc::Cosets cosets(Rank);
         cosets.add_row();
 
-        if (g.ngens == 0) {
+        if (Rank == 0) {
             return cosets;
         }
 
         for (int gen: sub_gens) {
-            if (gen < g.ngens)
+            if (gen < Rank)
                 cosets.put(0, gen, 0);
         }
 
-        Tables tables(g);
+        Tables<Rank> tables(g);
         tables.add_row();
         tables.initialize(0, cosets);
 
         std::priority_queue<int> facts;
 
         for (int idx = 0; idx < cosets.data.size(); idx++) {
-            int coset = idx / g.ngens;
-            int gen = idx % g.ngens;
+            int coset = idx / Rank;
+            int gen = idx % Rank;
             if (cosets.get(coset, gen) >= 0) continue;
 
             int target = cosets.size();
@@ -147,8 +188,8 @@ namespace tc {
                 int fact_idx = facts.top();
                 facts.pop();
 
-                int coset_ = fact_idx / g.ngens;
-                int gen_ = fact_idx % g.ngens;
+                int coset_ = fact_idx / Rank;
+                int gen_ = fact_idx % Rank;
 
                 if (cosets.get(coset_, gen_) != -1)
                     continue;
diff --git a/include/tc/groups.hpp b/include/tc/groups.hpp
index 8f60941..2ab3625 100644
--- a/include/tc/groups.hpp
+++ b/include/tc/groups.hpp
@@ -7,12 +7,11 @@ 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<int> &mults, const std::string &name) {
-        int ngens = (int) mults.size() + 1;
+    template<unsigned int Rank>
+    Group <Rank> schlafli(const std::array<unsigned int, Rank - 1> &mults, const std::string &name) {
+        Group<Rank> g(name);
 
-        Group g(ngens, name);
-
-        for (int i = 0; i < (int) mults.size(); i++) {
+        for (int i = 0; i < Rank - 1; i++) {
             g(i, i + 1) = mults[i];
         }
 
@@ -23,61 +22,68 @@ 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<int> &mults) {
+    template<unsigned int Rank>
+    Group <Rank> schlafli(const std::array<unsigned int, Rank - 1> &mults) {
         std::stringstream ss;
         ss << "[";
-        if (!mults.empty()) {
-            for (size_t i = 0; i < mults.size() - 1; ++i) {
+        if (Rank) {
+            for (size_t i = 0; i < Rank - 2; ++i) {
                 ss << mults[i] << ",";
             }
-            ss << mults.back();
+            ss << mults[Rank - 1];
         }
         ss << "]";
 
-        return schlafli(mults, ss.str());
+        return schlafli<Rank>(mults, ss.str());
     }
 
     namespace group {
         /**
          * Simplex
          */
-        Group A(int dim) {
+        template<unsigned int Rank>
+        Group <Rank> A() {
             std::stringstream ss;
-            ss << "A(" << dim << ")";
+            ss << "A(" << Rank << ")";
 
-            if (dim == 0)
-                return Group(0, ss.str());
+            if (Rank == 0)
+                return Group<Rank>(ss.str());
 
-            const std::vector<int> &mults = std::vector<int>(dim - 1, 3);
+            std::array<unsigned int, Rank - 1> mults;
+            mults.fill(3);
 
-            return schlafli(mults, ss.str());
+            return schlafli<Rank>(mults, ss.str());
         }
 
         /**
          * Cube, Orthoplex
          */
-        Group B(int dim) {
+        template<unsigned int Rank>
+        Group <Rank> B() {
             std::stringstream ss;
-            ss << "B(" << dim << ")";
+            ss << "B(" << Rank << ")";
 
-            std::vector<int> mults(dim - 1, 3);
+            std::array<unsigned int, Rank - 1> mults;
+            mults.fill(3);
             mults[0] = 4;
 
-            return schlafli(mults, ss.str());
+            return schlafli<Rank>(mults, ss.str());
         }
 
         /**
          * Demicube, Orthoplex
          */
-        Group D(int dim) {
+        template<unsigned int Rank>
+        Group <Rank> D() {
             std::stringstream ss;
-            ss << "D(" << dim << ")";
+            ss << "D(" << Rank << ")";
 
-            std::vector<int> mults(dim - 1, 3);
-            mults[dim - 2] = 2;
+            std::array<unsigned int, Rank - 1> mults;
+            mults.fill(3);
+            mults[Rank - 2] = 2;
 
-            Group g = schlafli(mults, ss.str());
-            g(1, dim - 1) = 3;
+            Group<Rank> g = schlafli<Rank>(mults, ss.str());
+            g(1, Rank - 1) = 3;
 
             return g;
         }
@@ -85,15 +91,17 @@ namespace tc {
         /**
          * E groups
          */
-        Group E(int dim) {
+        template<unsigned int Rank>
+        Group <Rank> E() {
             std::stringstream ss;
-            ss << "E(" << dim << ")";
+            ss << "E(" << Rank << ")";
 
-            std::vector<int> mults(dim - 1, 3);
-            mults[dim - 2] = 2;
+            std::array<unsigned int, Rank - 1> mults;
+            mults.fill(3);
+            mults[Rank - 2] = 2;
 
-            Group g = schlafli(mults, ss.str());
-            g(2, dim - 1) = 3;
+            Group<Rank> g = schlafli<Rank>(mults, ss.str());
+            g(2, Rank - 1) = 3;
 
             return g;
         }
@@ -101,58 +109,60 @@ namespace tc {
         /**
          * 24 Cell
          */
-        Group F4() {
-            return schlafli({3, 4, 3}, "F4");
+        Group<4> F4() {
+            return schlafli<4>({3, 4, 3}, "F4");
         }
 
         /**
          * Hexagon
          */
-        Group G2() {
-            return schlafli({6}, "G2");
+        Group<2> G2() {
+            return schlafli<2>({6}, "G2");
         }
 
         /**
          * Icosahedron
          */
-        Group H(int dim) {
+        template<unsigned int Rank>
+        Group <Rank> H() {
             std::stringstream ss;
-            ss << "H(" << dim << ")";
+            ss << "H(" << Rank << ")";
 
-            std::vector<int> mults(dim - 1, 3);
+            std::array<unsigned int, Rank - 1> mults;
+            mults.fill(3);
             mults[0] = 5;
 
-            return schlafli(mults, ss.str());
+            return schlafli<Rank>(mults, ss.str());
         }
 
         /**
          * Polygonal
          */
-        Group I2(int n) {
+        Group<2> I2(unsigned int n) {
             std::stringstream ss;
             ss << "I2(" << n << ")";
 
-            return schlafli({n}, ss.str());
+            return schlafli<2>({n}, ss.str());
         }
 
         /**
          * Toroidal. I2(n) * I2(m)
          */
-        Group T(int n, int m) {
+        Group<4> T(unsigned int n, unsigned int m) {
             std::stringstream ss;
             ss << "T(" << n << "," << m << ")";
 
-            return schlafli({n, 2, m}, ss.str());
+            return schlafli<4>({n, 2, m}, ss.str());
         }
 
         /**
          * Toroidal. T(n, n)
          */
-        Group T(int n) {
+        Group<4> T(unsigned int n) {
             std::stringstream ss;
             ss << "T(" << n << ")";
 
-            return schlafli({n, 2, n}, ss.str());
+            return schlafli<4>({n, 2, n}, ss.str());
         }
     }
 }