WIP: Copy tc logic from toddcox-faster

This commit is contained in:
David Allemang
2022-02-23 18:36:51 -05:00
parent 22976a9778
commit 63eb9e47b4
8 changed files with 730 additions and 1 deletions

161
tc/src/core.cpp Normal file
View File

@@ -0,0 +1,161 @@
#include "tc/core.hpp"
#include <utility>
#include <sstream>
#include <algorithm>
namespace tc {
Action::Action(int from_idx, int gen)
: from_idx(from_idx), gen(gen) {
}
void Path::add_row() {
path.resize(path.size() + 1);
}
Action Path::get(int to_idx) const {
return path[to_idx];
}
void Path::put(int from_idx, int gen, int to_idx) {
path[to_idx] = Action(from_idx, gen);
}
size_t Path::size() const {
return path.size();
}
Cosets::Cosets(int ngens)
: ngens(ngens) {
}
void Cosets::add_row() {
data.resize(data.size() + ngens, -1);
path.add_row();
}
void Cosets::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);
}
}
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.get(target).from_idx == -1) {
path.put(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;
}
SubGroup Group::subgroup(const std::vector<int> &gens) const {
return SubGroup(*this, gens);
}
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;
}
SubGroup::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.get(gen_map[i], gen_map[j]);
set(Rel(i, j, mult));
}
}
}
Group operator*(const Group &g, const Group &h) {
return g.product(h);
}
Group operator^(const Group &g, int p) {
return g.power(p);
}
}

120
tc/src/groups.cpp Normal file
View File

@@ -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());
}
}
}

189
tc/src/solve.cpp Normal file
View File

@@ -0,0 +1,189 @@
#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();
if (ngens == 0) {
return cosets;
}
for (int g : sub_gens) {
if (g < ngens)
cosets.put(0, 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;
if (target == coset)
for (int table_idx : gen_map[gen])
if (target_row.lst_ptrs[table_idx] == nullptr)
target_row.gnrs[table_idx] = -1;
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;
}
}