diff --git a/.gitignore b/.gitignore
new file mode 100644
index 0000000..6447641
--- /dev/null
+++ b/.gitignore
@@ -0,0 +1,4 @@
+.idea
+
+*.[oa]
+cmake-build*
diff --git a/CMakeLists.txt b/CMakeLists.txt
new file mode 100644
index 0000000..376e5e5
--- /dev/null
+++ b/CMakeLists.txt
@@ -0,0 +1,11 @@
+cmake_minimum_required(VERSION 3.10)
+project(toddcox-faster)
+
+add_library(tc STATIC
+        src/groups.cpp
+        src/solve.cpp
+        src/core.cpp)
+
+target_include_directories(tc PUBLIC include)
+
+add_subdirectory(example)
diff --git a/example/CMakeLists.txt b/example/CMakeLists.txt
new file mode 100644
index 0000000..8d87503
--- /dev/null
+++ b/example/CMakeLists.txt
@@ -0,0 +1,5 @@
+add_executable(bench bench.cpp)
+target_link_libraries(bench PRIVATE tc)
+
+add_executable(path path.cpp)
+target_link_libraries(path PRIVATE tc)
diff --git a/example/bench.cpp b/example/bench.cpp
new file mode 100644
index 0000000..4be5049
--- /dev/null
+++ b/example/bench.cpp
@@ -0,0 +1,35 @@
+#include "tc/core.hpp"
+#include "tc/groups.hpp"
+
+#include <ctime>
+#include <vector>
+#include <iostream>
+
+int main() {
+    std::vector<tc::Group> groups = {
+        tc::group::H(2),
+        tc::group::H(3),
+        tc::group::H(4),
+        tc::group::T(100),
+        tc::group::T(500),
+        tc::group::T(1000),
+        tc::group::E(6),
+        tc::group::E(7),
+        tc::group::B(6),
+        tc::group::B(7),
+        tc::group::B(8),
+    };
+
+    for (const auto &group : groups) {
+        auto s = std::clock(); // to measure CPU time
+        auto cosets = group.solve();
+        auto e = std::clock();
+
+        double diff = (double) (e - s) / CLOCKS_PER_SEC;
+        int order = cosets.size();
+
+        std::cout << group.name << "," << order << "," << diff << std::endl;
+    }
+
+    return 0;
+}
diff --git a/example/path.cpp b/example/path.cpp
new file mode 100644
index 0000000..53029c0
--- /dev/null
+++ b/example/path.cpp
@@ -0,0 +1,17 @@
+#include "tc/solver.h"
+#include "tc/groups.h"
+
+#include <ctime>
+#include <iostream>
+
+int main() {
+    auto cube = tc::group::B(3);
+    auto vars = tc::solve(cube, {});
+
+    for (int target = 1; target < vars.len; target++) {
+        auto &action = vars.path[target];
+        std::cout << action.coset << " " << action.gen << " " << target << std::endl;
+    }
+
+    return 0;
+}
diff --git a/groups.cpp b/groups.cpp
deleted file mode 100644
index db140e5..0000000
--- a/groups.cpp
+++ /dev/null
@@ -1,52 +0,0 @@
-#include "solver.cpp"
-
-Group A(const int n) {
-    if (n == 0)
-        return Group(0);
-
-    return Group::schlafli(std::vector<int>(n-1,3));
-}
-
-Group B(const int n) {
-    std::vector<int> mults(n-1,3);
-    mults[0] = 4;
-    return Group::schlafli(mults);
-}
-
-Group D(const int n) {
-    std::vector<int> mults(n-1,3);
-    mults[n-2] = 2;
-    Group g = Group::schlafli(mults);
-    g.setmult({1,n-1,3});
-    return g;
-}
-
-Group E(const int n) {
-    std::vector<int> mults(n-1,3);
-    mults[n-2] = 2;
-    Group g = Group::schlafli(mults);
-    g.setmult({2,n-1,3});
-    return g;
-}
-
-Group F4() {
-    return Group::schlafli({3,4,3});
-}
-
-Group G2() {
-    return Group::schlafli({6});
-}
-
-Group H(const int n) {
-    std::vector<int> mults(n-1,3);
-    mults[0] = 5;
-    return Group::schlafli(mults);
-}
-
-Group I2(const int n) {
-    return Group::schlafli({n});
-}
-
-Group T(const int n) {
-    return I2(n)^2;
-}
diff --git a/include/tc/core.hpp b/include/tc/core.hpp
new file mode 100644
index 0000000..ef196c5
--- /dev/null
+++ b/include/tc/core.hpp
@@ -0,0 +1,80 @@
+#pragma once
+
+#include<array>
+#include<vector>
+#include<string>
+
+namespace tc {
+    struct Action {
+        int coset = -1;
+        int gen = -1;
+        int target = -1;
+
+        Action() = default;
+
+        Action(const Action &) = default;
+
+        Action(int coset, int gen, int target);
+    };
+
+    struct Cosets {
+        int ngens;
+        std::vector<int> data;
+        std::vector<Action> path;
+
+        Cosets(const Cosets &) = default;
+
+        explicit Cosets(int ngens);
+
+        void add_row();
+
+        void put(int coset, int gen, int target);
+
+        void put(int idx, int target);
+
+        [[nodiscard]] int get(int coset, int gen) const;
+
+        [[nodiscard]] int get(int idx) const;
+
+        [[nodiscard]] size_t size() const;
+    };
+
+    struct Rel {
+        std::array<int, 2> gens;
+        int mult;
+
+        Rel() = default;
+
+        Rel(const Rel &) = default;
+
+        Rel(int a, int b, int m);
+
+        [[nodiscard]] Rel shift(int off) const;
+    };
+
+    struct Group {
+        const int ngens;
+        std::vector<std::vector<int>> _mults;
+        std::string name;
+
+        Group(const Group &) = default;
+
+        explicit Group(int ngens, const std::vector<Rel> &rels = {}, std::string name = "G");
+
+        void set(const Rel &r);
+
+        [[nodiscard]] int get(int a, int b) const;
+
+        [[nodiscard]] std::vector<Rel> rels() const;
+
+        [[nodiscard]] Group product(const Group &other) const;
+
+        [[nodiscard]] Group power(int p) const;
+
+        [[nodiscard]] Cosets solve(const std::vector<int>& sub_gens = {}) const;
+    };
+
+    Group operator*(const Group &g, const Group &h);
+
+    Group operator^(const Group &g, int p);
+}
diff --git a/include/tc/groups.hpp b/include/tc/groups.hpp
new file mode 100644
index 0000000..7a3e421
--- /dev/null
+++ b/include/tc/groups.hpp
@@ -0,0 +1,70 @@
+#pragma once
+
+#include "core.hpp"
+
+
+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);
+
+    /**
+     * 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);
+
+    namespace group {
+        /**
+         * Simplex
+         */
+        Group A(int dim);
+
+        /**
+         * Cube, Orthoplex
+         */
+        Group B(int dim);
+
+        /**
+         * Demicube, Orthoplex
+         */
+        Group D(int dim);
+
+        /**
+         * E groups
+         */
+        Group E(int dim);
+
+        /**
+         * 24 Cell
+         */
+        Group F4();
+
+        /**
+         * Hexagon
+         */
+        Group G2();
+
+        /**
+         * Icosahedron
+         */
+        Group H(int dim);
+
+        /**
+         * Polygonal
+         */
+        Group I2(int n);
+
+        /**
+         * Toroidal. I2(n) * I2(m)
+         */
+        Group T(int n, int m);
+
+        /**
+         * Toroidal. T(n, n)
+         */
+        Group T(int n);
+    }
+}
diff --git a/solver.cpp b/solver.cpp
deleted file mode 100644
index 4a69d16..0000000
--- a/solver.cpp
+++ /dev/null
@@ -1,290 +0,0 @@
-#include <vector>
-#include <algorithm>
-#include <queue>
-#include <utility>
-#include <functional>
-#include <iostream>
-
-struct Cosets {
-    int ngens;
-    std::vector<int> data;
-    int len;
-
-    Cosets(int ngens, std::vector<int> data): ngens(ngens), data(data) {
-        len = data.size() / ngens;
-    }
-
-    void add_row() {
-        len++;
-        data.resize(data.size()+ngens, -1);
-    }
-
-    void put(int coset, int gen, int target) {
-        data[coset * ngens + gen] = target;
-        data[target * ngens + gen] = coset;
-    }
-
-    void put(int idx, int target) {
-        int coset = idx / ngens;
-        int gen = idx % ngens;
-        data[idx] = target;
-        data[target * ngens + gen] = coset;
-    }
-
-    int get(int coset, int gen) {
-        return data[coset * ngens + gen];
-    }
-
-    int get(int idx) {
-        return data[idx];
-    }
-};
-
-struct Mult {
-    int gen0, gen1, mult;
-};
-
-struct RelTablesRow {
-    int* gnrs;
-    int** lst_ptrs;
-    RelTablesRow (int N, int* gnrs, int** lst_ptrs): gnrs(gnrs), lst_ptrs(lst_ptrs) {
-        for (int i = 0; i < N; i++) {
-            lst_ptrs[i] = nullptr;
-        }
-    }
-};
-
-struct TableInfo {
-    int gens[2];
-    int mult;
-    TableInfo(Mult m) {
-        gens[0] = m.gen0;
-        gens[1] = m.gen1;
-        mult = m.mult;
-    }
-};
-
-struct RelTables {
-    static const int ROW_BLOCK_SIZE = 128;
-    std::vector<TableInfo> table_info;
-    std::vector<RelTablesRow*> rows;
-    int start = 0;
-    int num_tables;
-    int buffer_rows = 0;
-
-    RelTables (std::vector<Mult> mults): num_tables(mults.size()) {
-        for (Mult m : mults) {
-            table_info.emplace_back(m);
-        }
-    }
-
-    void add_row() {
-        if (buffer_rows == 0) {
-            int* gnrs_alloc = new int[num_tables*RelTables::ROW_BLOCK_SIZE];
-            int** lst_ptrs_alloc = new int*[num_tables*RelTables::ROW_BLOCK_SIZE];
-            for (int i = 0; i < RelTables::ROW_BLOCK_SIZE; i++) {
-                rows.push_back(new RelTablesRow(num_tables, &gnrs_alloc[i*num_tables], &lst_ptrs_alloc[i*num_tables]));
-            }
-            buffer_rows = RelTables::ROW_BLOCK_SIZE;
-        }
-
-        buffer_rows--;
-    }
-
-    void del_rows_to(int idx) {
-        const int del_to = (idx/RelTables::ROW_BLOCK_SIZE)*RelTables::ROW_BLOCK_SIZE;
-        for (int i = start; i < del_to; i += RelTables::ROW_BLOCK_SIZE) {
-            delete[] rows[i]->gnrs;
-            delete[] rows[i]->lst_ptrs;
-            for (int j = 0; j < RelTables::ROW_BLOCK_SIZE; j++) {
-                delete rows[i+j];
-            }
-            start += RelTables::ROW_BLOCK_SIZE;
-        }
-    }
-};
-
-struct Group {
-    int ngens;
-    std::vector<std::vector<int>> _mults;
-
-    Group(int ngens, std::vector<Mult> rels = {}): ngens(ngens) {
-        _mults.resize(ngens);
-        for (int i = 0; i < ngens; i++) {
-            _mults[i].resize(ngens, 2);
-        }
-
-        for (Mult m : rels) {
-            if (m.gen0 < m.gen1)
-                _mults[m.gen0][m.gen1] = m.mult;
-            else
-                _mults[m.gen1][m.gen0] = m.mult;
-        }
-    }
-
-    void setmult(Mult m) {
-        if (m.gen0 < m.gen1)
-            _mults[m.gen0][m.gen1] = m.mult;
-        else
-            _mults[m.gen1][m.gen0] = m.mult;
-    }
-
-    std::vector<Mult> get_mults() const {
-        std::vector<Mult> mults; 
-        for (int i = 0; i < ngens - 1; i++) {
-            for (int j = i+1; j < ngens; j++) {
-                mults.push_back({i,j,_mults[i][j]});
-            }
-        }
-        return mults;
-    }
-
-    Group operator*(const Group &other) {
-        int off = ngens;
-
-        Group g(ngens + other.ngens, get_mults());
-        
-        for (Mult m : other.get_mults()) {
-            g.setmult({off + m.gen0, off + m.gen1, m.mult});
-        }
-
-        return g;
-    }
-
-    Group operator^(int p) {
-        Group g(ngens * p);
-
-        for (Mult m : get_mults()) {
-            for (int off = 0; off < g.ngens; off += ngens) {
-                g.setmult({off + m.gen0, off + m.gen1, m.mult});
-            }
-        }
-
-        return g;
-    }
-
-    static Group schlafli(std::vector<int> mults) {
-        int ngens = mults.size() + 1;
-        Group g(ngens);
-        for (int i = 0; i < mults.size(); i++) {
-            g.setmult({i, i+1, mults[i]});
-        }
-        return g;
-    }
-
-    Cosets solve(std::vector<int> sub_gens = {}) {
-        std::vector<int> init_row(ngens, -1);
-        for (int i : sub_gens) {
-            init_row[i] = 0;
-        }
-
-        Cosets cosets(ngens, init_row);
-        RelTables rel_tables(get_mults());
-        std::vector<std::vector<int>> gen_map(ngens);
-        int rel_idx = 0;
-        for (Mult m : get_mults()) {
-            gen_map[m.gen0].push_back(rel_idx);
-            gen_map[m.gen1].push_back(rel_idx);
-            rel_idx++;
-        }
-
-        int null_lst_ptr;
-        rel_tables.add_row();
-        RelTablesRow &row = *(rel_tables.rows[0]);
-        for (int table_idx = 0; table_idx < rel_tables.num_tables; table_idx++) {
-            TableInfo &ti = rel_tables.table_info[table_idx];
-            
-            if (cosets.get(ti.gens[0]) + cosets.get(ti.gens[1]) == -2) {
-                row.lst_ptrs[table_idx] = new int;
-                row.gnrs[table_idx] = 0;
-            }
-            else {
-                row.lst_ptrs[table_idx] = &null_lst_ptr;
-                row.gnrs[table_idx] = -1;
-            }
-        }
-
-        int idx = 0;
-        int coset, gen, target, fact_idx, lst, gen_;
-        while (true) {
-            while (idx < cosets.data.size() and cosets.get(idx) >= 0)
-                idx++;
-
-            if (idx == cosets.data.size()) {
-                rel_tables.del_rows_to(idx / ngens);
-                break;
-            }
-
-            target = cosets.len;
-            cosets.add_row();
-            rel_tables.add_row();
-            
-            std::vector<int> facts;
-            facts.push_back(idx);
-
-            coset = idx / ngens;
-            gen = idx % ngens;
-
-            rel_tables.del_rows_to(coset);
-
-            RelTablesRow &target_row = *(rel_tables.rows[target]);
-            while (!facts.empty()) {
-                fact_idx = facts.back();
-                facts.pop_back();
-
-                if (cosets.get(fact_idx) != -1)
-                    continue;
-
-                cosets.put(fact_idx, target);
-
-                coset = fact_idx / ngens;
-                gen = fact_idx % ngens;
-
-                RelTablesRow &coset_row = *(rel_tables.rows[coset]);
-                for (int table_idx : gen_map[gen]) {
-                    if ( target_row.lst_ptrs[table_idx] == nullptr ) {
-                        TableInfo &ti = rel_tables.table_info[table_idx];
-                        target_row.lst_ptrs[table_idx] = coset_row.lst_ptrs[table_idx];
-                        target_row.gnrs[table_idx] = coset_row.gnrs[table_idx] + 1;
-                        
-                        if (coset_row.gnrs[table_idx] < 0)
-                            target_row.gnrs[table_idx] -= 2;
-
-                        if (target_row.gnrs[table_idx] == ti.mult) {
-                            lst = *(target_row.lst_ptrs[table_idx]);
-                            delete target_row.lst_ptrs[table_idx];
-                            gen_ = ti.gens[(int)(ti.gens[0] == gen)];
-                            facts.push_back(lst*ngens + gen_);
-                        }
-                        else if (target_row.gnrs[table_idx] == -ti.mult) {
-                            gen_ = ti.gens[ti.gens[0] == gen];
-                            facts.push_back(target*ngens + gen_);
-                        }
-                        else if (target_row.gnrs[table_idx] == ti.mult - 1) {
-                            *(target_row.lst_ptrs[table_idx]) = target;
-                        }
-                    }
-                }
-
-                std::sort(facts.begin(), facts.end(), std::greater<int>());
-            }
-
-            for (int table_idx = 0; table_idx < rel_tables.num_tables; table_idx++) {
-                TableInfo &ti = rel_tables.table_info[table_idx];
-                if (target_row.lst_ptrs[table_idx] == nullptr) {
-                    if ( (cosets.get(target, ti.gens[0]) != target) and 
-                         (cosets.get(target, ti.gens[1]) != target) ) {
-                        target_row.lst_ptrs[table_idx] = new int;
-                        target_row.gnrs[table_idx] = 0;
-                    }
-                    else {
-                        target_row.lst_ptrs[table_idx] = &null_lst_ptr;
-                        target_row.gnrs[table_idx] = -1;
-                    }
-                }
-            }
-        }
-        
-        return cosets;
-    }
-};
diff --git a/src/core.cpp b/src/core.cpp
new file mode 100644
index 0000000..dcffb3b
--- /dev/null
+++ b/src/core.cpp
@@ -0,0 +1,118 @@
+#include "tc/core.hpp"
+
+#include <utility>
+#include <sstream>
+
+namespace tc {
+    Action::Action(int coset, int gen, int target)
+        : coset(coset), gen(gen), target(target) {
+    }
+
+    Cosets::Cosets(int ngens)
+        : ngens(ngens) {
+    }
+
+    void Cosets::add_row() {
+        data.resize(data.size() + ngens, -1);
+        path.resize(path.size() + 1);
+    }
+
+    void Cosets::put(int coset, int gen, int target) {
+        data[coset * ngens + gen] = target;
+        data[target * ngens + gen] = coset;
+
+        if (path[target].coset == -1) {
+            path[target] = Action(coset, gen, target);
+        }
+    }
+
+    void Cosets::put(int idx, int target) {
+        int coset = idx / ngens;
+        int gen = idx % ngens;
+        data[idx] = target;
+        data[target * ngens + gen] = coset;
+
+        if (path[target].coset == -1) {
+            path[target] = Action(coset, gen, target);
+        }
+    }
+
+    int Cosets::get(int coset, int gen) const {
+        return data[coset * ngens + gen];
+    }
+
+    int Cosets::get(int idx) const {
+        return data[idx];
+    }
+
+    size_t Cosets::size() const {
+        return path.size();
+    }
+
+    Rel::Rel(int a, int b, int m)
+        : gens({a, b}), mult(m) {
+    }
+
+    Rel Rel::shift(int off) const {
+        return Rel(gens[0] + off, gens[1] + off, mult);
+    }
+
+    Group::Group(int ngens, const std::vector<Rel> &rels, std::string name)
+        : ngens(ngens), name(std::move(name)) {
+        _mults.resize(ngens);
+
+        for (auto &mult : _mults) {
+            mult.resize(ngens, 2);
+        }
+
+        for (const auto &rel : rels) {
+            set(rel);
+        }
+    }
+
+    void Group::set(const Rel &r) {
+        _mults[r.gens[0]][r.gens[1]] = r.mult;
+        _mults[r.gens[1]][r.gens[0]] = r.mult;
+    }
+
+    int Group::get(int a, int b) const {
+        return _mults[a][b];
+    }
+
+    std::vector<Rel> Group::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, get(i, j));
+            }
+        }
+        return res;
+    }
+
+    Group Group::product(const Group &other) const {
+        std::stringstream ss;
+        ss << name << "*" << other.name;
+
+        Group g(ngens + other.ngens, rels(), ss.str());
+
+        for (const auto &rel : other.rels()) {
+            g.set(rel.shift(ngens));
+        }
+
+        return g;
+    }
+
+    Group Group::power(int p) const {
+        std::stringstream ss;
+        ss << name << "^" << p;
+
+        Group g(ngens * p, {}, ss.str());
+        for (const auto &rel : rels()) {
+            for (int off = 0; off < g.ngens; off += ngens) {
+                g.set(rel.shift(off));
+            }
+        }
+
+        return g;
+    }
+}
\ No newline at end of file
diff --git a/src/groups.cpp b/src/groups.cpp
new file mode 100644
index 0000000..1cebc12
--- /dev/null
+++ b/src/groups.cpp
@@ -0,0 +1,120 @@
+#include "tc/groups.hpp"
+
+#include <sstream>
+
+namespace tc {
+    Group schlafli(const std::vector<int> &mults, const std::string &name) {
+        int ngens = (int) mults.size() + 1;
+
+        Group g(ngens, {}, name);
+
+        for (int i = 0; i < (int) mults.size(); i++) {
+            g.set(Rel(i, i + 1, mults[i]));
+        }
+
+        return g;
+    }
+
+    Group schlafli(const std::vector<int> &mults) {
+        std::stringstream ss;
+        ss << "[";
+        if (!mults.empty()) {
+            for (size_t i = 0; i < mults.size() - 1; ++i) {
+                ss << mults[i] << ",";
+            }
+            ss << mults.back();
+        }
+        ss << "]";
+
+        return schlafli(mults, ss.str());
+    }
+
+    namespace group {
+        Group A(const int dim) {
+            std::stringstream ss;
+            ss << "A(" << dim << ")";
+
+            if (dim == 0)
+                return Group(0, {}, ss.str());
+
+            const std::vector<int> &mults = std::vector<int>(dim - 1, 3);
+
+            return schlafli(mults, ss.str());
+        }
+
+        Group B(const int dim) {
+            std::stringstream ss;
+            ss << "B(" << dim << ")";
+
+            std::vector<int> mults(dim - 1, 3);
+            mults[0] = 4;
+
+            return schlafli(mults, ss.str());
+        }
+
+        Group D(const int dim) {
+            std::stringstream ss;
+            ss << "D(" << dim << ")";
+
+            std::vector<int> mults(dim - 1, 3);
+            mults[dim - 2] = 2;
+
+            Group g = schlafli(mults, ss.str());
+            g.set(Rel(1, dim - 1, 3));
+
+            return g;
+        }
+
+        Group E(const int dim) {
+            std::stringstream ss;
+            ss << "E(" << dim << ")";
+
+            std::vector<int> mults(dim - 1, 3);
+            mults[dim - 2] = 2;
+
+            Group g = schlafli(mults, ss.str());
+            g.set(Rel(2, dim - 1, 3));
+
+            return g;
+        }
+
+        Group F4() {
+            return schlafli({3, 4, 3}, "F4");
+        }
+
+        Group G2() {
+            return schlafli({6}, "G2");
+        }
+
+        Group H(const int dim) {
+            std::stringstream ss;
+            ss << "H(" << dim << ")";
+
+            std::vector<int> mults(dim - 1, 3);
+            mults[0] = 5;
+
+            return schlafli(mults, ss.str());
+        }
+
+        Group I2(const int n) {
+            std::stringstream ss;
+            ss << "I2(" << n << ")";
+
+            return schlafli({n}, ss.str());
+        }
+
+        Group T(const int n, const int m) {
+            std::stringstream ss;
+            ss << "T(" << n << "," << m << ")";
+
+            return schlafli({n, 2, m}, ss.str());
+        }
+
+        Group T(const int n) {
+            std::stringstream ss;
+            ss << "T(" << n << ")";
+
+            return schlafli({n, 2, n}, ss.str());
+        }
+    }
+}
diff --git a/src/solve.cpp b/src/solve.cpp
new file mode 100644
index 0000000..f56f51a
--- /dev/null
+++ b/src/solve.cpp
@@ -0,0 +1,178 @@
+#include "tc/core.hpp"
+
+#include <algorithm>
+
+namespace tc {
+    struct RelTablesRow {
+        int *gnrs;
+        int **lst_ptrs;
+
+        RelTablesRow(int N, int *gnrs, int **lst_ptrs) : gnrs(gnrs), lst_ptrs(lst_ptrs) {
+            for (int i = 0; i < N; i++) {
+                lst_ptrs[i] = nullptr;
+            }
+        }
+    };
+
+    struct RelTables {
+        static const int ROW_BLOCK_SIZE = 64;
+        std::vector<Rel> rels;
+        std::vector<RelTablesRow *> rows;
+        int start = 0;
+        int num_tables;
+        int buffer_rows = 0;
+
+        explicit RelTables(const std::vector<Rel> &rels)
+            : num_tables(rels.size()), rels(rels) {
+        }
+
+        void add_row() {
+            if (buffer_rows == 0) {
+                int *gnrs_alloc = new int[num_tables * RelTables::ROW_BLOCK_SIZE];
+                int **lst_ptrs_alloc = new int *[num_tables * RelTables::ROW_BLOCK_SIZE];
+                for (int i = 0; i < RelTables::ROW_BLOCK_SIZE; i++) {
+                    rows.push_back(
+                        new RelTablesRow(num_tables, &gnrs_alloc[i * num_tables], &lst_ptrs_alloc[i * num_tables]));
+                }
+                buffer_rows = RelTables::ROW_BLOCK_SIZE;
+            }
+
+            buffer_rows--;
+        }
+
+        void del_rows_to(int idx) {
+            const int del_to = (idx / RelTables::ROW_BLOCK_SIZE) * RelTables::ROW_BLOCK_SIZE;
+            for (int i = start; i < del_to; i += RelTables::ROW_BLOCK_SIZE) {
+                delete[] rows[i]->gnrs;
+                delete[] rows[i]->lst_ptrs;
+                for (int j = 0; j < RelTables::ROW_BLOCK_SIZE; j++) {
+                    delete rows[i + j];
+                }
+                start += RelTables::ROW_BLOCK_SIZE;
+            }
+        }
+
+        ~RelTables() {
+            while (start < rows.size()) {
+                delete[] rows[start]->gnrs;
+                delete[] rows[start]->lst_ptrs;
+                for (int j = 0; j < RelTables::ROW_BLOCK_SIZE; j++) {
+                    delete rows[start + j];
+                }
+                start += RelTables::ROW_BLOCK_SIZE;
+            }
+        }
+    };
+
+    Cosets Group::solve(const std::vector<int> &sub_gens) const {
+        Cosets cosets(ngens);
+        cosets.add_row();
+        for (int g : sub_gens) {
+            cosets.put(g, 0);
+        }
+
+        RelTables rel_tables(rels());
+        std::vector<std::vector<int>> gen_map(ngens);
+        int rel_idx = 0;
+        for (Rel m : rels()) {
+            gen_map[m.gens[0]].push_back(rel_idx);
+            gen_map[m.gens[1]].push_back(rel_idx);
+            rel_idx++;
+        }
+
+        int null_lst_ptr;
+        rel_tables.add_row();
+        RelTablesRow &row = *(rel_tables.rows[0]);
+        for (int table_idx = 0; table_idx < rel_tables.num_tables; table_idx++) {
+            Rel &ti = rel_tables.rels[table_idx];
+
+            if (cosets.get(ti.gens[0]) + cosets.get(ti.gens[1]) == -2) {
+                row.lst_ptrs[table_idx] = new int;
+                row.gnrs[table_idx] = 0;
+            } else {
+                row.lst_ptrs[table_idx] = &null_lst_ptr;
+                row.gnrs[table_idx] = -1;
+            }
+        }
+
+        int idx = 0;
+        int coset, gen, target, fact_idx, lst, gen_;
+        while (true) {
+            while (idx < cosets.data.size() and cosets.get(idx) >= 0)
+                idx++;
+
+            if (idx == cosets.data.size()) {
+                rel_tables.del_rows_to(idx / ngens);
+                break;
+            }
+
+            target = cosets.size();
+            cosets.add_row();
+            rel_tables.add_row();
+
+            std::vector<int> facts;
+            facts.push_back(idx);
+
+            coset = idx / ngens;
+            gen = idx % ngens;
+
+            rel_tables.del_rows_to(coset);
+
+            RelTablesRow &target_row = *(rel_tables.rows[target]);
+            while (!facts.empty()) {
+                fact_idx = facts.back();
+                facts.pop_back();
+
+                if (cosets.get(fact_idx) != -1)
+                    continue;
+
+                cosets.put(fact_idx, target);
+
+                coset = fact_idx / ngens;
+                gen = fact_idx % ngens;
+
+                RelTablesRow &coset_row = *(rel_tables.rows[coset]);
+                for (int table_idx : gen_map[gen]) {
+                    if (target_row.lst_ptrs[table_idx] == nullptr) {
+                        Rel &ti = rel_tables.rels[table_idx];
+                        target_row.lst_ptrs[table_idx] = coset_row.lst_ptrs[table_idx];
+                        target_row.gnrs[table_idx] = coset_row.gnrs[table_idx] + 1;
+
+                        if (coset_row.gnrs[table_idx] < 0)
+                            target_row.gnrs[table_idx] -= 2;
+
+                        if (target_row.gnrs[table_idx] == ti.mult) {
+                            lst = *(target_row.lst_ptrs[table_idx]);
+                            delete target_row.lst_ptrs[table_idx];
+                            gen_ = ti.gens[(int) (ti.gens[0] == gen)];
+                            facts.push_back(lst * ngens + gen_);
+                        } else if (target_row.gnrs[table_idx] == -ti.mult) {
+                            gen_ = ti.gens[ti.gens[0] == gen];
+                            facts.push_back(target * ngens + gen_);
+                        } else if (target_row.gnrs[table_idx] == ti.mult - 1) {
+                            *(target_row.lst_ptrs[table_idx]) = target;
+                        }
+                    }
+                }
+
+                std::sort(facts.begin(), facts.end(), std::greater<>());
+            }
+
+            for (int table_idx = 0; table_idx < rel_tables.num_tables; table_idx++) {
+                Rel &ti = rel_tables.rels[table_idx];
+                if (target_row.lst_ptrs[table_idx] == nullptr) {
+                    if ((cosets.get(target, ti.gens[0]) != target) and
+                        (cosets.get(target, ti.gens[1]) != target)) {
+                        target_row.lst_ptrs[table_idx] = new int;
+                        target_row.gnrs[table_idx] = 0;
+                    } else {
+                        target_row.lst_ptrs[table_idx] = &null_lst_ptr;
+                        target_row.gnrs[table_idx] = -1;
+                    }
+                }
+            }
+        }
+
+        return cosets;
+    }
+}
\ No newline at end of file
diff --git a/test.cpp b/test.cpp
deleted file mode 100644
index 114f7c5..0000000
--- a/test.cpp
+++ /dev/null
@@ -1,19 +0,0 @@
-#include "groups.cpp"
-#include <chrono>
-#include <iostream>
-
-int main() {
-    Group g = B(8);
-
-    auto s = std::chrono::system_clock::now();
-    auto cosets = g.solve();
-    auto e = std::chrono::system_clock::now();
-
-    std::chrono::duration<double> diff = e - s;
-    int order = cosets.len;
-
-    std::cout << order << std::endl;
-    std::cout << diff.count() << std::endl;
-
-    return 0;
-}