diff --git a/example/CMakeLists.txt b/example/CMakeLists.txt
index 8d87503..63a2aa8 100644
--- a/example/CMakeLists.txt
+++ b/example/CMakeLists.txt
@@ -3,3 +3,6 @@ target_link_libraries(bench PRIVATE tc)
 
 add_executable(path path.cpp)
 target_link_libraries(path PRIVATE tc)
+
+add_executable(group group.cpp)
+target_link_libraries(group PRIVATE tc)
diff --git a/example/group.cpp b/example/group.cpp
new file mode 100644
index 0000000..8adb148
--- /dev/null
+++ b/example/group.cpp
@@ -0,0 +1,34 @@
+#include <iostream>
+
+#include <tc/solver.hpp>
+#include <tc/groups.hpp>
+
+int main() {
+    constexpr unsigned int Rank = 5;
+    constexpr unsigned int SRank = 3;
+
+    tc::Group<Rank> group = tc::schlafli<Rank>({5, 3, 2, 3});
+//    tc::Group<Rank> group = tc::group::E<Rank>();
+
+    tc::SubGroups<Rank, SRank> subs = tc::subgroups<Rank, SRank>(group);
+
+    std::cout << "Group " << group.name << " (" << subs.size() << " subgroups)" << std::endl;
+    std::cout << group << std::endl;
+
+    for (const auto &sub: subs) {
+        for (int i = 0; i < SRank; ++i) {
+            for (int j = 0; j < SRank; ++j) {
+                auto sub_mult = sub(i, j);
+                auto src_mult = group(sub.gens(i), sub.gens(j));
+
+                if (sub_mult != src_mult) {
+                    std::cout << "Incorrect subgroup " << sub.name << std::endl;
+                    std::cout << sub << std::endl;
+                    return 1;
+                }
+            }
+        }
+    }
+
+    return 0;
+}
diff --git a/include/tc/group.hpp b/include/tc/group.hpp
index 092858f..5a6b17b 100644
--- a/include/tc/group.hpp
+++ b/include/tc/group.hpp
@@ -2,6 +2,7 @@
 
 #include <string>
 #include <sstream>
+#include <numeric>
 
 #include <Eigen/Eigen>
 
@@ -9,25 +10,95 @@ namespace {
     template<class T>
     std::string stringify(const T &vec) {
         std::stringstream ss;
-        ss << "[" << vec << "]";
+        ss << "[" << vec.transpose() << "]";
         return ss.str();
     }
 }
 
 namespace tc {
+    template<unsigned int Rank>
+    using Symbol = Eigen::Vector<unsigned int, Rank>;
+
+    template<unsigned int Rank>
+    Symbol<Rank> iota() {
+        Symbol<Rank> res;
+        for (int i = 0; i < Rank; ++i) {
+            res(i) = i;
+        }
+        return res;
+    }
+
     /// A Coxeter Matrix
     template<unsigned int Rank>
     class Group : public Eigen::Matrix<unsigned int, Rank, Rank> {
     public:
         using Base = Eigen::Matrix<unsigned int, Rank, Rank>;
 
-        std::string name;
+        std::string name = "G";
+        Symbol<Rank> gens = iota<Rank>();
 
         using Base::Base;
     };
 
-    template<unsigned int Rank>
-    using Symbol = Eigen::Vector<unsigned int, Rank>;
+    template<unsigned int Rank, unsigned int SRank>
+    Group<SRank> subgroup(const Group<Rank> &group, Symbol<SRank> gens) {
+        Group<SRank> res;
+        res.name = group.name + ":" + stringify(gens);
+        res.gens = gens;
+
+        for (int i = 0; i < SRank; ++i) {
+            for (int j = 0; j < SRank; ++j) {
+                res(i, j) = group(gens[i], gens[j]);
+            }
+        }
+
+        return res;
+    }
+
+    template<unsigned int Rank, unsigned int SRank>
+    Symbol<Rank> inverse(Symbol<SRank> gens) {
+        Symbol<Rank> res;
+        res.fill(0);
+        for (int i = 0; i < SRank; ++i) {
+            res(gens(i)) = i;
+        }
+        return res;
+    }
+
+    constexpr unsigned int Factorial(unsigned int n) {
+        unsigned int res = 1;
+        for (int i = 1; i <= n; ++i) {
+            res *= i;
+        }
+        return res;
+    }
+
+    constexpr unsigned int Choose(unsigned int n, unsigned int k) {
+        return Factorial(n) / Factorial(k) / Factorial(n - k);
+    }
+
+    template<unsigned int Rank, unsigned int SRank>
+    using SubGroups = std::array<Group<SRank>, Choose(Rank, SRank)>;
+
+    template<unsigned int Rank, unsigned int SRank>
+    SubGroups<Rank, SRank> subgroups(const Group<Rank> &group) {
+        std::vector<bool> mask(Rank, false);
+        std::fill_n(mask.begin(), SRank, true);
+
+        SubGroups<Rank, SRank> res;
+
+        size_t i = 0;
+        Symbol<SRank> row;
+        do {
+            for (int j = 0, k = 0; j < Rank; ++j) {
+                if (mask[j])
+                    row(k++) = j;
+            }
+            res[i++] = subgroup<Rank, SRank>(group, row);
+        } while (std::prev_permutation(mask.begin(), mask.end()));
+
+        return res;
+    }
 
     /**
      * Create a named coxeter matrix from a simplified schlafli symbol