me/toddcox-faster
1
0
Fork 0
forked from mirror/toddcox-faster
Code Pull Requests Activity

Un-template groups.

Browse Source 
Templating the groups doesn't improve performance that much, but it severely hinders usability. Switch things to be dynamically sized so that different ranked groups are the same type.

Could also now separate headers and implementation... but will probably leave it for now.
This commit is contained in:
David Allemang
2021-11-07 20:40:34 -05:00
parent 24e7ab47d1
commit fa13493569
10 changed files with 222 additions and 237 deletions
Download Patch File Download Diff File Expand all files Collapse all files

7
Testing/common.hpp
View File

@@ -21,11 +21,10 @@ std::vector<unsigned int> parse_vec(const std::string &part) {
return res;
}
template<unsigned int N>
int compute(
const tc::Group<N> &group,
size_t compute(
const tc::Group &group,
const std::vector<unsigned int> &vgens
) {
auto table = tc::solve<N>(group, vgens);
auto table = tc::solve(group, vgens);
return table.order();
}

32
Testing/schlafli.cpp
View File

@@ -7,34 +7,10 @@ int main(int argc, char **argv) {
auto vgens = parse_vec(argv[2]);
auto target = std::stoul(argv[3]);
int order;
switch (vsymbol.size()) {
case 0: {
Eigen::Vector<unsigned int, 0> symbol(vsymbol.data());
order = compute<1>(tc::schlafli<1>(symbol), vgens);
break;
}
case 1: {
Eigen::Vector<unsigned int, 1> symbol(vsymbol.data());
order = compute<2>(tc::schlafli<2>(symbol), vgens);
break;
}
case 2: {
Eigen::Vector<unsigned int, 2> symbol(vsymbol.data());
order = compute<3>(tc::schlafli<3>(symbol), vgens);
break;
}
case 3: {
Eigen::Vector<unsigned int, 3> symbol(vsymbol.data());
order = compute<4>(tc::schlafli<4>(symbol), vgens);
break;
}
case 4: {
Eigen::Vector<unsigned int, 4> symbol(vsymbol.data());
order = compute<5>(tc::schlafli<5>(symbol), vgens);
break;
}
}
tc::Symbol symbol(vsymbol.size());
symbol << Eigen::Map<tc::Symbol>(vsymbol.data(), vsymbol.size());
tc::Group group = tc::schlafli(symbol);
auto order = compute(group, vgens);
std::cout << "Order: " << order << ":" << target << std::endl;
return order != target;

22
Testing/special.cpp
View File

@@ -9,31 +9,29 @@ int main(int argc, char **argv) {
unsigned int order;
tc::Group group(0);
if (name == "E6") {
auto group = tc::group::E<6>();
order = compute<6>(group, vgens);
group = tc::group::E(6);
}
if (name == "E7") {
auto group = tc::group::E<7>();
order = compute<7>(group, vgens);
group = tc::group::E(7);
}
if (name == "E8") {
auto group = tc::group::E<8>();
order = compute<8>(group, vgens);
group = tc::group::E(8);
}
if (name == "B6") {
auto group = tc::group::B<6>();
order = compute<6>(group, vgens);
group = tc::group::B(6);
}
if (name == "B7") {
auto group = tc::group::B<7>();
order = compute<7>(group, vgens);
group = tc::group::B(7);
}
if (name == "B8") {
auto group = tc::group::B<8>();
order = compute<8>(group, vgens);
group = tc::group::B(8);
}
order = compute(group, vgens);
std::cout << "Order: " << order << ":" << target << std::endl;
return order != target;
}

16
examples/bench.cpp
View File

@@ -22,17 +22,17 @@ void test(const G &group) {
}
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;
}

15
examples/group.cpp
View File

@@ -4,20 +4,21 @@
#include <tc/groups.hpp>
int main() {
constexpr unsigned int Rank = 5;
constexpr unsigned int SRank = 3;
tc::Symbol symbol(4);
symbol << 5, 3, 2, 3;
tc::Group group = tc::schlafli(symbol);
// tc::Group group = tc::group::E(7);
tc::Group<Rank> group = tc::schlafli<Rank>({5, 3, 2, 3});
// tc::Group<Rank> group = tc::group::E<Rank>();
size_t srank = 3;
tc::SubGroups<Rank, SRank> subs = tc::subgroups<Rank, SRank>(group);
tc::SubGroups subs = tc::subgroups(group, srank);
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) {
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));

2
examples/path.cpp
View File

@@ -6,7 +6,7 @@
#include <tc/groups.hpp>
int main() {
auto cube = tc::group::B<3>();
auto cube = tc::group::B(3);
auto vars = tc::solve(cube);
std::string start;

48
include/tc/cosets.hpp
View File

@@ -3,58 +3,65 @@
#include <vector>
namespace tc {
template<unsigned int Rank>
class Path;
template<unsigned int Rank>
class Cosets {
private:
std::vector<int> data;
size_t _rank;
public:
Cosets(const Cosets<Rank> &) = default;
Cosets(const Cosets &) = default;
Cosets() = default;
explicit Cosets(size_t rank) : _rank(rank) {}
void add_row() {
data.resize(data.size() + Rank, -1);
data.resize(data.size() + rank(), -1);
}
void put(int coset, int gen, int target) {
data[coset * Rank + gen] = target;
data[target * Rank + gen] = coset;
data[coset * rank() + gen] = target;
data[target * rank() + gen] = coset;
}
[[nodiscard]] int get(int coset, int gen) const {
return data[coset * Rank + gen];
return data[coset * rank() + gen];
}
[[nodiscard]] size_t rank() const {
return _rank;
}
[[nodiscard]] size_t order() const {
return data.size() / Rank;
return data.size() / _rank;
}
Path<Rank> path() const;
Path path() const;
};
template<unsigned int Rank>
class Path {
private:
friend class Cosets<Rank>;
friend class Cosets;
std::vector<unsigned int> source;
std::vector<unsigned int> gen;
size_t _rank;
size_t _order;
explicit Path(size_t order) : _order(order), source(order), gen(order) {}
explicit Path(size_t rank, size_t order) : _rank(rank), _order(order), source(order), gen(order) {}
public:
size_t order() const {
[[nodiscard]] size_t rank() const {
return _rank;
}
[[nodiscard]] size_t order() const {
return _order;
}
template<class T, class F>
std::vector <T> walk(const T &start, const F &op) {
std::vector <T> res;
std::vector<T> walk(const T &start, const F &op) {
std::vector<T> res;
res.reserve(order());
res.push_back(start);
@@ -67,20 +74,19 @@ namespace tc {
}
template<class T, class E, class F>
std::vector <T> walk(const T &start, const E &gens, const F &op) {
std::vector<T> walk(const T &start, const E &gens, const F &op) {
return walk(start, [&](const T &s, const int g) {
return op(s, gens[g]);
});
}
};
template<unsigned int Rank>
Path<Rank> Cosets<Rank>::path() const {
Path<Rank> res(order());
Path Cosets::path() const {
Path res(rank(), order());
std::vector<bool> set(order());
for (int coset = 0; coset < order(); ++coset) {
for (int gen = 0; gen < Rank; ++gen) {
for (int gen = 0; gen < rank(); ++gen) {
int target = get(coset, gen);
if (!set[target]) {

122
include/tc/group.hpp
View File

@@ -16,38 +16,39 @@ namespace {
}
namespace tc {
template<unsigned int Rank>
using Symbol = Eigen::Vector<unsigned int, Rank>;
using Symbol = Eigen::Vector<unsigned int, Eigen::Dynamic>;
template<unsigned int Rank>
Symbol<Rank> iota() {
Symbol<Rank> res;
for (int i = 0; i < Rank; ++i) {
res(i) = i;
}
return res;
}
using MatrixXui = Eigen::Matrix<unsigned int, Eigen::Dynamic, Eigen::Dynamic>;
/// A Coxeter Matrix
template<unsigned int Rank>
class Group : public Eigen::Matrix<unsigned int, Rank, Rank> {
class Group : public MatrixXui {
public:
using Base = Eigen::Matrix<unsigned int, Rank, Rank>;
using Base = MatrixXui;
std::string name = "G";
Symbol<Rank> gens = iota<Rank>();
Symbol gens;
using Base::Base;
explicit Group(size_t rank) : Base(rank, rank), gens(rank) {
for (Eigen::Index i = 0; i < rank; ++i) {
gens(i) = i;
}
}
[[nodiscard]] size_t rank() const {
return rows();
}
};
template<unsigned int Rank, unsigned int SRank>
Group<SRank> subgroup(const Group<Rank> &group, Symbol<SRank> gens) {
Group<SRank> res;
Group subgroup(const Group &group, const Symbol &gens) {
size_t rank = group.size();
size_t srank = gens.size();
Group res(srank);
res.name = group.name + ":" + stringify(gens);
res.gens = gens;
for (int i = 0; i < SRank; ++i) {
for (int j = 0; j < SRank; ++j) {
for (Eigen::Index i = 0; i < srank; ++i) {
for (Eigen::Index j = 0; j < srank; ++j) {
res(i, j) = group(gens[i], gens[j]);
}
}
@@ -55,17 +56,17 @@ namespace tc {
return res;
}
template<unsigned int Rank, unsigned int SRank>
Symbol<Rank> inverse(Symbol<SRank> gens) {
Symbol<Rank> res;
Symbol inverse(size_t rank, const Symbol &gens) {
size_t srank = gens.size();
Symbol res(rank);
res.fill(0);
for (int i = 0; i < SRank; ++i) {
for (int i = 0; i < srank; ++i) {
res(gens(i)) = i;
}
return res;
}
constexpr unsigned int Factorial(unsigned int n) {
unsigned int factorial(unsigned int n) {
unsigned int res = 1;
for (int i = 1; i <= n; ++i) {
res *= i;
@@ -73,28 +74,28 @@ namespace tc {
return res;
}
constexpr unsigned int Choose(unsigned int n, unsigned int k) {
return Factorial(n) / Factorial(k) / Factorial(n - k);
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)>;
using SubGroups = std::vector<Group>;
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 subgroups(const Group &group, size_t srank) {
size_t rank = group.rank();
SubGroups<Rank, SRank> res;
std::vector<bool> mask(rank, false);
std::fill_n(mask.begin(), srank, true);
size_t i = 0;
Symbol<SRank> row;
SubGroups res;
res.reserve(choose(rank, srank));
Symbol row(srank);
do {
for (int j = 0, k = 0; j < Rank; ++j) {
for (int j = 0, k = 0; j < rank; ++j) {
if (mask[j])
row(k++) = j;
}
res[i++] = subgroup<Rank, SRank>(group, row);
res.push_back(subgroup(group, row));
} while (std::prev_permutation(mask.begin(), mask.end()));
return res;
@@ -103,15 +104,16 @@ namespace tc {
/**
* Create a named coxeter matrix from a simplified schlafli symbol
*/
template<unsigned int Rank>
Group<Rank> schlafli(const Symbol<Rank - 1> &mults, const std::string &name) {
Group<Rank> res;
Group schlafli(const Symbol &mults, const std::string &name) {
size_t rank = mults.size() + 1;
Group res(rank);
res.name = name;
res.fill(2);
res.diagonal().fill(1);
res.topRightCorner(Rank - 1, Rank - 1).diagonal() << mults;
res.bottomLeftCorner(Rank - 1, Rank - 1).diagonal() << mults;
res.topRightCorner(rank - 1, rank - 1).diagonal() << mults;
res.bottomLeftCorner(rank - 1, rank - 1).diagonal() << mults;
return res;
}
@@ -119,38 +121,36 @@ namespace tc {
/**
* Create a coxeter matrix from a simplified schlafli symbol.
*/
template<unsigned int Rank>
Group<Rank> schlafli(const Symbol<Rank - 1> &mults) {
return schlafli<Rank>(mults, stringify(mults));
Group schlafli(const Symbol &mults) {
return schlafli(mults, stringify(mults));
}
template<unsigned int GR, unsigned int HR>
Group<GR + HR> product(const Group<GR> &g, const Group<HR> &h) {
Group<GR + HR> res;
Group product(const Group &g, const Group &h) {
Group res(g.rank() + h.rank());
res.name = g.name + "*" + h.name;
res.fill(2);
int off = 0;
res.block(off, off, GR, GR) << g.array() + off;
off += GR;
Eigen::Index off = 0;
res.block(off, off, g.rank(), g.rank()) << g.array() + off;
off += (Eigen::Index) g.rank();
res.block(off, off, HR, HR) << h.array() + off;
off += HR;
res.block(off, off, h.rank(), h.rank()) << h.array() + off;
off += (Eigen::Index) h.rank();
return res;
}
template<unsigned int GR, unsigned int P>
Group<GR * P> power(const Group<GR> &g) {
Group<GR * P> res;
res.name = g.name + "^" + P;
Group power(const Group &g, size_t p) {
Group res(g.rank() * p);
res.name = g.name + "^" + std::to_string(p);
res.fill(2);
for (int k = 0; k < P; ++k) {
int off = k * GR;
res.block(off, off, GR, GR) << g.array() + off;
for (Eigen::Index k = 0; k < p; ++k) {
auto off = (Eigen::Index) g.rank() * k;
res.block(off, off, g.rank(), g.rank()) << g.array() + off;
}
return res;

133
include/tc/groups.hpp
View File

@@ -6,13 +6,11 @@ namespace tc::group {
/**
* Universal Coxeter Group
*/
template<unsigned int Rank>
Group<Rank> U() {
std::stringstream ss;
ss << "U(" << Rank << ")";
Group U(size_t rank) {
std::string name = "U(" + std::to_string(rank) + ")";
Group<Rank> res;
res.name = ss.str();
Group res(rank);
res.name = name;
res.fill(2);
return res;
}
@@ -20,53 +18,47 @@ namespace tc::group {
/**
* Simplex
*/
template<unsigned int Rank>
Group<Rank> A() {
std::stringstream ss;
ss << "A(" << Rank << ")";
Group A(size_t rank) {
std::string name = "A(" + std::to_string(rank) + ")";
if (Rank == 0) {
Group<Rank> res;
res.name = ss.str();
if (rank == 0) {
Group res(rank);
res.name = name;
return res;
}
tc::Symbol<Rank - 1> mults;
mults.fill(3);
tc::Symbol symbol(rank - 1);
symbol.fill(3);
return schlafli<Rank>(mults, ss.str());
return schlafli(symbol, name);
}
/**
* Cube, Orthoplex
*/
template<unsigned int Rank>
Group<Rank> B() {
std::stringstream ss;
ss << "B(" << Rank << ")";
Group B(size_t rank) {
std::string name = "B(" + std::to_string(rank) + ")";
tc::Symbol<Rank - 1> mults;
mults.fill(3);
mults(0) = 4;
tc::Symbol symbol(rank - 1);
symbol.fill(3);
symbol(0) = 4;
return schlafli<Rank>(mults, ss.str());
return schlafli(symbol, name);
}
/**
* Demicube, Orthoplex
*/
template<unsigned int Rank>
Group<Rank> D() {
std::stringstream ss;
ss << "D(" << Rank << ")";
Group D(size_t rank) {
std::string name = "D(" + std::to_string(rank) + ")";
tc::Symbol<Rank - 1> mults;
mults.fill(3);
mults(Rank - 2) = 2;
tc::Symbol symbol(rank - 1);
symbol.fill(3);
symbol((Eigen::Index) rank - 2) = 2;
Group<Rank> g = schlafli<Rank>(mults, ss.str());
g(1, Rank - 1) = 3;
g(Rank - 1, 1) = 3;
Group g = schlafli(symbol, name);
g(1, (Eigen::Index) rank - 1) = 3;
g((Eigen::Index) rank - 1, 1) = 3;
return g;
}
@@ -74,18 +66,16 @@ namespace tc::group {
/**
* E groups
*/
template<unsigned int Rank>
Group<Rank> E() {
std::stringstream ss;
ss << "E(" << Rank << ")";
Group E(size_t rank) {
std::string name = "E(" + std::to_string(rank) + ")";
tc::Symbol<Rank - 1> mults;
mults.fill(3);
mults(Rank - 2) = 2;
tc::Symbol symbol(rank - 1);
symbol.fill(3);
symbol((Eigen::Index) rank - 2) = 2;
Group<Rank> g = schlafli<Rank>(mults, ss.str());
g(2, Rank - 1) = 3;
g(Rank - 1, 2) = 3;
Group g = schlafli(symbol, name);
g(2, (Eigen::Index) rank - 1) = 3;
g((Eigen::Index) rank - 1, 2) = 3;
return g;
}
@@ -93,59 +83,66 @@ namespace tc::group {
/**
* 24 Cell
*/
Group<4> F4() {
return schlafli<4>({3, 4, 3}, "F4");
Group F4() {
tc::Symbol symbol(3);
symbol << 3, 4, 3;
return schlafli(symbol, "F4");
}
/**
* Hexagon
*/
Group<2> G2() {
return schlafli<2>(tc::Symbol<1>(6), "G2");
Group G2() {
tc::Symbol symbol(1);
symbol << 6;
return schlafli(symbol, "G2");
}
/**
* Icosahedron
*/
template<unsigned int Rank>
Group<Rank> H() {
std::stringstream ss;
ss << "H(" << Rank << ")";
Group H(size_t rank) {
std::string name = "H(" + std::to_string(rank) + ")";
tc::Symbol<Rank - 1> mults;
mults.fill(3);
mults(0) = 5;
tc::Symbol symbol(rank - 1);
symbol.fill(3);
symbol(0) = 5;
return schlafli<Rank>(mults, ss.str());
return schlafli(symbol, name);
}
/**
* Polygonal
*/
Group<2> I2(unsigned int n) {
std::stringstream ss;
ss << "I2(" << n << ")";
Group I2(unsigned int n) {
std::string name = "I2(" + std::to_string(n) + ")";
return schlafli<2>(tc::Symbol<1>(n), ss.str());
tc::Symbol symbol(1);
symbol << n;
return schlafli(symbol, name);
}
/**
* Toroidal. I2(n) * I2(m)
*/
Group<4> T(unsigned int n, unsigned int m) {
std::stringstream ss;
ss << "T(" << n << "," << m << ")";
Group T(unsigned int n, unsigned int m) {
std::string name = "T(" + std::to_string(n) + "," + std::to_string(m) + ")";
return schlafli<4>({n, 2, m}, ss.str());
tc::Symbol symbol(3);
symbol << n, 2, m;
return schlafli(symbol, name);
}
/**
* Toroidal. T(n, n)
*/
Group<4> T(unsigned int n) {
std::stringstream ss;
ss << "T(" << n << ")";
Group T(unsigned int n) {
std::string name = "T(" + std::to_string(n) + ")";
return schlafli<4>({n, 2, n}, ss.str());
tc::Symbol symbol(3);
symbol << n, 2, n;
return schlafli(symbol, name);
}
}

62
include/tc/solver.hpp
View File

@@ -20,7 +20,7 @@ namespace {
public:
int i, j, mult;
std::vector <Row> rows;
std::vector<Row> rows;
public:
explicit Table(int i, int j, int mult) :
@@ -60,31 +60,38 @@ namespace {
}
};
template<unsigned int Rank>
class Tables {
public:
static constexpr unsigned int Rels = Rank * (Rank + 1) / 2 - Rank;
private:
int *null_lst_ptr = new int;
BlockAllocator<int> alloc;
std::array <std::shared_ptr<Table>, Rels> tables;
std::array <std::vector<std::shared_ptr < Table>>, Rank>
deps;
std::vector<std::shared_ptr<Table>> tables;
std::vector<std::vector<std::shared_ptr<Table>>> deps;
size_t _rank;
size_t _rels;
public:
explicit Tables(const tc::Group<Rank> &group) {
for (int i = 0, irel = 0; i < Rank - 1; ++i) {
for (int j = i + 1; j < Rank; ++j, ++irel) {
explicit Tables(const tc::Group &group) : _rank(group.rank()), _rels(rank() * (rank() + 1) / 2 - rank()) {
deps.resize(rank());
for (int i = 0; i < rank() - 1; ++i) {
for (int j = i + 1; j < rank(); ++j) {
auto table = std::make_shared<Table>(i, j, group(i, j));
tables[irel] = table;
tables.push_back(table);
deps[i].push_back(table);
deps[j].push_back(table);
}
}
}
[[nodiscard]] size_t rank() const {
return _rank;
}
[[nodiscard]] size_t rels() const {
return _rels;
}
void add_row() {
// std::vector already does block allocation.
for (const auto &table: tables) {
@@ -92,7 +99,7 @@ namespace {
}
}
void initialize(int target, const tc::Cosets<Rank> &cosets) {
void initialize(int target, const tc::Cosets &cosets) {
for (auto &table: tables) {
Row &row = table->rows[target];
@@ -113,7 +120,7 @@ namespace {
delete null_lst_ptr;
}
void learn(int coset, int gen, int target, const tc::Cosets<Rank> &cosets, std::priority_queue<int> &facts) {
void learn(int coset, int gen, int target, const tc::Cosets &cosets, std::priority_queue<size_t> &facts) {
if (target == coset) {
for (auto &table: deps[gen]) {
Row &target_row = table->rows[target];
@@ -140,11 +147,11 @@ namespace {
// forward learn
int lst = *target_row.lst;
int gen_ = (table->i == gen) ? table->j : table->i;
facts.push(lst * Rank + gen_);
facts.push(lst * rank() + gen_);
} else if (target_row.gnr == -table->mult) {
// stationary learn
int gen_ = (table->i == gen) ? table->j : table->i;
facts.push(target * Rank + gen_);
facts.push(target * rank() + gen_);
} else if (target_row.gnr == table->mult - 1) {
// determined family
*target_row.lst = target;
@@ -159,35 +166,36 @@ namespace tc {
/**
* Assumes that g is a coxeter group - that is, self-adjoint and the diagonal is 2.
*/
template<unsigned int Rank>
tc::Cosets<Rank> solve(const Group <Rank> &group, const std::vector<unsigned int> &sub_gens = {}) {
tc::Cosets<Rank> cosets;
tc::Cosets solve(const Group &group, const std::vector<unsigned int> &sub_gens = {}) {
size_t rank = group.rank();
tc::Cosets cosets(rank);
cosets.add_row();
if (Rank == 0) {
if (rank == 0) {
return cosets;
}
for (unsigned int gen: sub_gens) {
if (gen < Rank)
if (gen < rank)
cosets.put(0, gen, 0);
}
Tables<Rank> tables(group);
Tables tables(group);
tables.add_row();
tables.initialize(0, cosets);
std::priority_queue<int> facts;
std::priority_queue<size_t> facts;
for (int coset = 0; coset < cosets.order(); coset++) {
for (int gen = 0; gen < Rank; ++gen) {
for (int gen = 0; gen < rank; ++gen) {
if (cosets.get(coset, gen) >= 0) continue; // todo vector<bool> set
int target = cosets.order();
cosets.add_row();
tables.add_row();
facts.push(coset * Rank + gen);
facts.push(coset * rank + gen);
// todo nothing before the current coset will be used.
// delete all table rows using old cosets to free memory early.
@@ -198,8 +206,8 @@ namespace tc {
int fact_idx = facts.top();
facts.pop();
int coset_ = fact_idx / Rank;
int gen_ = fact_idx % Rank;
int coset_ = fact_idx / rank;
int gen_ = fact_idx % rank;
if (cosets.get(coset_, gen_) != -1)
continue;