Introduce tc::subgroup and tc::subgroups

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