int p;

int cells;

bool squared;

eGeometry base;

vector<primeinfo> primes;

vector<int> dualval;

int current_prime_id;

fgeomextra(eGeometry b, int i) : base(b), current_prime_id(i) {}

bool operator == (const matrix& B) const {

for(int i=0; i<MWDIM; i++) for(int j=0; j<MWDIM; j++)

if(self[i][j] != B[i][j]) return false;

return true;

}

bool operator != (const matrix& B) const {

for(int i=0; i<MWDIM; i++) for(int j=0; j<MWDIM; j++)

if(self[i][j] != B[i][j]) return true;

return false;

}

bool operator < (const matrix& B) const {

for(int i=0; i<MWDIM; i++) for(int j=0; j<MWDIM; j++)

if(self[i][j] != B[i][j]) return self[i][j] < B[i][j];

return false;

}

unsigned force_hash;

int Prime, wsquare, Field, dual;

// we perform our computations in the field Z_Prime[w] where w^2 equals wsquare

// (or simply Z_Prime for wsquare == 0)

#define EASY

// 'easy' assumes that all elements of the field actually used

// are of form n or mw (not n+mw), and cs and ch are both of form n

// by experimentation, such cs and ch always exist

// many computations are much simpler under that assumption

#ifndef EASY

static int neasy;

int m(int x) { x %= Prime; if(x<0) x+= Prime; return x; }

#endif

int sub(int a, int b) {

#ifdef EASY

return (a + b * (Prime-1)) % Prime;

#else

return m(a%Prime-b%Prime) + Prime * m(a/Prime-b/Prime);

#endif

}

int add(int a, int b) {

#ifdef EASY

if(a == ERR || b == ERR || a*b<0) return ERR;

return (a+b)%Prime;

#else

return m(a%Prime+b%Prime) + Prime * m(a/Prime+b/Prime);

#endif

}

int mul(int tx, int ty) {

#ifdef EASY

return (tx*ty*((tx<0&&ty<0)?wsquare:1)) % Prime;

#else

if(tx >= Prime && tx % Prime) neasy++;

if(ty >= Prime && ty % Prime) neasy++;

int x[2], y[2], z[3];

for(int i=0; i<3; i++) z[i] = 0;

for(int i=0; i<2; i++)

x[i] = tx%Prime, tx /= Prime;

for(int i=0; i<2; i++)

y[i] = ty%Prime, ty /= Prime;

for(int i=0; i<2; i++)

for(int j=0; j<2; j++)

z[i+j] = (z[i+j] + x[i] * y[j]) % Prime;

z[0] += z[2] * wsquare;

return m(z[0]) + Prime * m(z[1]);

#endif

}

int sqr(int x) { return mul(x,x); }

int err;

matrix mmul(const matrix& A, const matrix& B) {

matrix res;

for(int i=0; i<MWDIM; i++) for(int k=0; k<MWDIM; k++) {

int t = 0;

#ifdef EASY

int tp = 0, tn = 0;

for(int j=0; j<MWDIM; j++) {

int val = mul(A[i][j], B[j][k]);

if(val > 0) tp += val;

else tn += val;

}

tp %= Prime; tn %= Prime;

if(tp && tn) err++;

t = tp + tn;

#else

for(int j=0; j<MWDIM; j++) t = add(t, mul(A[i][j], B[j][k]));

#endif

res[i][k] = t;

}

return res;

}

map<matrix, int> matcode;

vector<matrix> matrices;

vector<string> qpaths;

vector<matrix> qcoords;

// S7 in 2D, but e.g. 4 for a 3D cube

int rotations;

// S7 in 2D, but e.g. 24 for a 3D cube

int local_group;

// Id: Identity

// R : rotate by 1/rotations of the full circle

// P : make a step and turn backwards

// X : in 3-dim, turn by 90 degrees

matrix Id, R, P, X;

matrix strtomatrix(string s) {

matrix res = Id;

matrix m = Id;

for(int i=isize(s)-1; i>=0; i--)

if(s[i] == 'R') res = mmul(R, res);

else if (s[i] == 'P') res = mmul(P, res);

else if (s[i] == 'x') { m[0][0] = -1; res = mmul(m, res); m[0][0] = +1; }

else if (s[i] == 'y') { m[1][1] = -1; res = mmul(m, res); m[1][1] = +1; }

else if (s[i] == 'z') { m[2][2] = -1; res = mmul(m, res); m[2][2] = +1; }

return res;

}

void addas(const matrix& M, int i) {

if(!matcode.count(M)) {

matcode[M] = i;

for(int j=0; j<isize(qcoords); j++)

addas(mmul(M, qcoords[j]), i);

}

}

void add(const matrix& M) {

if(!matcode.count(M)) {

int i = isize(matrices);

matcode[M] = i, matrices.push_back(M);

for(int j=0; j<isize(qcoords); j++)

addas(mmul(M, qcoords[j]), i);

if(WDIM == 3) add(mmul(X, M));

add(mmul(R, M));

}

}

#define MXF 1000000

vector<int> connections;

vector<int> inverses; // NYI in 3D

// 2D only

vector<int> rrf; // rrf[i] equals gmul(i, rotations-1)

vector<int> rpf; // rpf[i] equals gmul(i, rotations)

matrix mpow(matrix M, int N) {

while((N&1) == 0) N >>= 1, M = mmul(M, M);

matrix res = M;

N >>= 1;

while(N) {

M = mmul(M,M); if(N&1) res = mmul(res, M);

N >>= 1;

}

return res;

}

int gmul(int a, int b) { return matcode[mmul(matrices[a], matrices[b])]; }

int gpow(int a, int N) { return matcode[mpow(matrices[a], N)]; }

int gorder(int a) {

int b = a;

int q = 1;

while(b) b = gmul(b, a), q++;

return q;

}

pair<int,bool> gmul(pair<int, bool> a, int b) {

return make_pair(gmul(a.first,b), a.second);

}

int order(const matrix& M);

string decodepath(int i) {

string s;

while(i) {

if(i % S7) i--, s += 'R';

else i = connections[i], s += 'P';

}

return s;

}

int orderstats();

int cs, sn, ch, sh;

int solve();

void build();

static constexpr int MAXDIST = 120;

vector<char> disthep;

vector<char> disthex;

vector<char> distwall, distriver, distwall2, distriverleft, distriverright, distflower;

int distflower0;

vector<eItem> markers;

eVariation current_variation;

int getdist(pair<int,bool> a, vector<char>& dists);

int getdist(pair<int,bool> a, pair<int,bool> b);

int dijkstra(vector<char>& dists, vector<int> indist[MAXDIST]);

void analyze();

int maxdist, otherpole, circrad, wallid, wallorder, riverid;

bool easy(int i) {

return i < Prime || !(i % Prime);

}

// 11 * 25

// (1+z+z^3) * (1+z^3+z^4) ==

// 1+z+z^7 == 1+z+z^2(z^5) == 1+z+z^2(1+z^2) = 1+z+z^2+z^4

void init(int p) {

indenter_finish iinit(debug_geometry, "initializing a fieldpattern");

Prime = p;

if(solve()) {

printf("error: could not solve the fieldpattern\n");

exit(1);

}

build();

analyze();

}

fpattern(int p) {

force_hash = 0;

#if CAP_THREAD && MAXMDIM >= 4

dis = nullptr;

#endif

if(!p) return;

init(p);

}

void findsubpath();

vector<matrix> generate_isometries();

bool check_order(matrix M, int req);

unsigned compute_hash();

void set_field(int p, int sq);

unsigned hashv;

#if MAXMDIM >= 4

// general 4D

vector<transmatrix> fullv;

void add1(const matrix& M);

void add1(const matrix& M, const transmatrix& Full);

vector<matrix> generate_isometries3();

int solve3();

bool generate_all3();

#if CAP_THREAD

struct discovery *dis;

#endif

#endif

vector<triplet_info> find_triplets();

void generate_quotientgroup();

fpattern experiment;

std::shared_ptr<std::thread> discoverer;

std::mutex lock;

std::condition_variable cv;

bool is_suspended;

bool stop_it;

map<unsigned, tuple<int, int, matrix, matrix, matrix, int> > hashes_found;

discovery() : experiment(0) { is_suspended = false; stop_it = false; experiment.dis = this; experiment.Prime = experiment.Field = experiment.wsquare = 0; }

void activate();

void suspend();

void check_suspend();

void schedule_destruction();

void discovered();

~discovery();

int err = 0;

matrix P = M;

for(int i=1; i<req; i++) {

if(P == Id) return false;

P = mmul(P, M);

}

return P == Id && !err;

matrix T = Id;

int low = wsquare ? 1-Prime : 0;

vector<matrix> res;

auto colprod = [&] (int a, int b) {

return add(add(mul(T[0][a], T[0][b]), mul(T[1][a], T[1][b])), mul(T[2][a], T[2][b]));

};

for(T[0][0]=low; T[0][0]<Prime; T[0][0]++)

for(T[1][0]=low; T[1][0]<Prime; T[1][0]++)

for(T[2][0]=low; T[2][0]<Prime; T[2][0]++)

if(colprod(0, 0) == 1)

for(T[0][1]=low; T[0][1]<Prime; T[0][1]++)

for(T[1][1]=low; T[1][1]<Prime; T[1][1]++)

for(T[2][1]=low; T[2][1]<Prime; T[2][1]++)

if(colprod(1, 1) == 1)

if(colprod(1, 0) == 0)

for(T[0][2]=low; T[0][2]<Prime; T[0][2]++)

for(T[1][2]=low; T[1][2]<Prime; T[1][2]++)

for(T[2][2]=low; T[2][2]<Prime; T[2][2]++)

if(colprod(2, 2) == 1)

if(colprod(2, 0) == 0)

if(colprod(2, 1) == 0)

res.push_back(T);

return res;

matrix T = Id;

int low = wsquare ? 1-Prime : 0;

vector<matrix> res;

auto colprod = [&] (int a, int b) {

return add(add(mul(T[0][a], T[0][b]), mul(T[1][a], T[1][b])), sub(mul(T[2][a], T[2][b]), mul(T[3][a], T[3][b])));

};

auto rowcol = [&] (int a, int b) {

return add(add(mul(T[a][0], T[0][b]), mul(T[a][1], T[1][b])), add(mul(T[a][2], T[2][b]), mul(T[a][3], T[3][b])));

};

for(T[0][0]=low; T[0][0]<Prime; T[0][0]++)

for(T[1][0]=low; T[1][0]<Prime; T[1][0]++)

for(T[2][0]=low; T[2][0]<Prime; T[2][0]++)

for(T[3][0]=low; T[3][0]<Prime; T[3][0]++)

if(colprod(0, 0) == 1)

for(T[0][1]=low; T[0][1]<Prime; T[0][1]++)

for(T[1][1]=low; T[1][1]<Prime; T[1][1]++)

for(T[2][1]=low; T[2][1]<Prime; T[2][1]++)

for(T[3][1]=low; T[3][1]<Prime; T[3][1]++)

if(colprod(1, 1) == 1)

if(colprod(1, 0) == 0) {

#if CAP_THREAD && MAXMDIM >= 4

if(dis) dis->check_suspend();

if(dis && dis->stop_it) return res;

#endif

for(T[0][2]=low; T[0][2]<Prime; T[0][2]++)

for(T[0][3]=low; T[0][3]<Prime; T[0][3]++)

if(rowcol(0, 0) == 1)

if(rowcol(0, 1) == 0)

for(T[1][2]=low; T[1][2]<Prime; T[1][2]++)

for(T[1][3]=low; T[1][3]<Prime; T[1][3]++)

if(rowcol(1, 0) == 0)

if(rowcol(1, 1) == 1)

for(T[2][2]=low; T[2][2]<Prime; T[2][2]++)

for(T[3][2]=low; T[3][2]<Prime; T[3][2]++)

if(colprod(2, 2) == 1)

if(colprod(2, 0) == 0)

if(colprod(2, 1) == 0)

for(T[2][3]=low; T[2][3]<Prime; T[2][3]++)

for(T[3][3]=low; T[3][3]<Prime; T[3][3]++)

if(rowcol(2, 0) == 0)

if(rowcol(2, 1) == 0)

if(rowcol(2, 2) == 1)

// if(colprod(3, 3) == 1)

if(add(colprod(3, 3), 1) == 0)

if(colprod(3, 0) == 0)

if(colprod(3, 1) == 0)

if(colprod(3, 2) == 0)

if(rowcol(3, 3) == 1)

if(rowcol(3, 0) == 0)

if(rowcol(3, 1) == 0)

if(rowcol(3, 2) == 0)

res.push_back(T);

if(isize(res) > limitp) return res;

}

return res;

if(!matcode.count(M)) {

int i = isize(matrices);

matcode[M] = i, matrices.push_back(M);

}

if(!matcode.count(M)) {

int i = isize(matrices);

matcode[M] = i, matrices.push_back(M), fullv.push_back(Full);

}

unsigned hashv = 0;

int iR = matcode[R];

int iP = matcode[P];

int iX = matcode[X];

for(int i=0; i<isize(matrices); i++) {

hashv = 3 * hashv + gmul(i, iP) + 7 * gmul(i, iR);

if(MWDIM == 4) hashv += 11 * gmul(i, iX);

}

return hashv;

reg3::generate_fulls();

err = 0;

matrices.clear();

matcode.clear();

add1(Id);

fullv = {hr::Id};

for(int i=0; i<isize(matrices); i++) {

add1(mmul(matrices[i], R), fullv[i] * cgi.full_R);

add1(mmul(matrices[i], X), fullv[i] * cgi.full_X);

if(err) return false;

}

local_group = isize(matrices);

if(local_group != isize(cgi.cellrotations)) return false;

for(int i=0; i<(int)matrices.size(); i++) {

matrix E = mmul(matrices[i], P);

if(!matcode.count(E))

for(int j=0; j<local_group; j++) add1(mmul(E, matrices[j]));

if(err) return false;

if(isize(matrices) >= limitv) { println(hlog, "limitv exceeded"); return false; }

}

hashv = compute_hash();

if(debug_field)

println(hlog, "all = ", isize(matrices), "/", local_group, " = ", isize(matrices) / local_group, " hash = ", hashv, " count = ", ++hash_found[hashv]);

if(use_quotient_fp)

generate_quotientgroup();

return true;

int MS = isize(matrices);

int best_p = 0, best_i = 0;

for(int i=0; i<MS; i++) {

int j = i, p = 1;

while(j >= local_group)

j = gmul(j, i), p++;

if(j == 0 && p > best_p) {

bool okay = true;

vector<bool> visited(MS, false);

for(int ii=0; ii<MS; ii++) if(!visited[ii]) {

int jj = ii;

for(int k=0; k<p; k++) {

if(k && jj/local_group == ii/local_group) okay = false;

visited[jj] = true;

jj = gmul(i, jj);

}

}

if(okay) {

bool chk = (MS/p) % local_group;

println(hlog, "quotient by ", i, " : ", p, " times less, ", (MS/p/local_group), " tiles, check ", chk);

best_p = p; best_i = i;

if(chk) {

exit(1);

}

}

}

}

if(best_p > 1) {

vector<int> new_id(MS, -1);

vector<int> orig_id(MS, -1);

vector<matrix> new_matrices;

int nv = 0;

for(int i=0; i<MS; i++) if(new_id[i] == -1) {

int prode = i;

for(int l=0; l<local_group; l++) {

new_matrices.push_back(matrices[i+l]);

}

for(int k=0; k<best_p; k++) {

for(int l=0; l<local_group; l++) {

new_id[gmul(prode, l)] = nv + l;

}

prode = gmul(best_i, prode);

}

nv += local_group;

}

println(hlog, "got nv = ", nv, " / ", local_group);

for(int i=0; i<MS; i++)

matcode[matrices[i]] = new_id[i];

matrices = std::move(new_matrices);

println(hlog, "size matrices = ", isize(matrices), " size matcode = ", isize(matcode));

println(hlog, tie(P, R, X));

/*println(hlog, "TRY AGAIN");

}

indenter_finish isolve3(debug_field, "fpattern::solve3");

reg3::generate_fulls();

if(debug_field) println(hlog, "generating isometries for ", Field);

auto iso3 = generate_isometries();

auto iso4 = generate_isometries3();

int cmb = 0;

vector<matrix> possible_P, possible_X, possible_R;

for(auto& M: iso3) {

if(check_order(M, 2))

possible_X.push_back(M);

if(check_order(M, cgi.r_order))

possible_R.push_back(M);

}

for(auto& M: iso4)

if(check_order(M, 2))

possible_P.push_back(M);

if(debug_field)

println(hlog, "field = ", Field, " #P = ", isize(possible_P), " #X = ", isize(possible_X), " #R = ", isize(possible_R), " r_order = ", cgi.r_order, " xp_order = ", cgi.xp_order);

for(auto& xX: possible_X)

for(auto& xP: possible_P) if(check_order(mmul(xP, xX), cgi.xp_order))

for(auto& xR: possible_R) if(check_order(mmul(xR, xX), cgi.rx_order)) {

err = 0;

if(mmul(xX, xP) != mmul(xR, mmul(mmul(xP, xX), xR))) continue;

if(err) continue;

#if CAP_THREAD && MAXMDIM >= 4

if(dis) dis->check_suspend();

if(dis && dis->stop_it) return 0;

#endif

P = xP; R = xR; X = xX;

if(!generate_all3()) continue;

callhooks(hooks_solve3);

#if CAP_THREAD && MAXMDIM >= 4

if(dis) { dis->discovered(); continue; }

#endif

if(force_hash && hashv != force_hash) continue;

cmb++;

goto ok;

}

ok:

if(debug_field)

println(hlog, "cmb = ", cmb, " for field = ", Field);

return cmb;

Prime = p;

Field = sq ? Prime*Prime : Prime;

wsquare = sq;

for(int a=0; a<MWDIM; a++) for(int b=0; b<MWDIM; b++) Id[a][b] = a==b?1:0;

indenter_finish isolve(debug_field, "fpattern::solve");

for(int a=0; a<MWDIM; a++) for(int b=0; b<MWDIM; b++) Id[a][b] = a==b?1:0;

if(!isprime(Prime)) {

return 1;

}

rotations = WDIM == 2 ? S7 : 4;

local_group = WDIM == 2 ? S7 : 24;

for(dual=0; dual<3; dual++) {

for(int pw=1; pw<3; pw++) {

if(pw>3) break;

Field = pw==1? Prime : Prime*Prime;

if(pw == 2) {

for(wsquare=1; wsquare<Prime; wsquare++) {

int roots = 0;

for(int a=0; a<Prime; a++) if((a*a)%Prime == wsquare) roots++;

if(!roots) break;

}

} else wsquare = 0;

#if MAXMDIM >= 4

if(WDIM == 3) {

if(dual == 0 && (Prime <= limitsq || pw == 1)) {

int s = solve3();

if(s) return 0;

}

continue;

}

#endif

if(dual == 2) {

if(Field <= 10) {

vector<matrix> all_isometries = generate_isometries();

for(auto& X: all_isometries)

if(check_order(X, rotations))

for(auto& Y: all_isometries)

if(check_order(Y, 2) && check_order(mmul(X, Y), S3)) {

R = X; P = Y;

return 0;

}

}

continue;

}

#ifdef EASY

std::vector<int> sqrts(Prime, 0);

for(int k=1-Prime; k<Prime; k++) sqrts[sqr(k)] = k;

int fmax = Prime;

#else

std::vector<int> sqrts(Field);

for(int k=0; k<Field; k++) sqrts[sqr(k)] = k;

int fmax = Field;

#endif

R = P = X = Id;

X[1][1] = 0; X[2][2] = 0;

X[1][2] = 1; X[2][1] = Prime-1;

for(cs=0; cs<fmax; cs++) {

int sb = sub(1, sqr(cs));

sn = sqrts[sb];

R[0][0] = cs; R[1][1] = cs;

R[0][1] = sn; R[1][0] = sub(0, sn);

if(!check_order(R, dual ? S3 : rotations)) continue;

if(R[0][0] == 1) continue;

for(ch=2; ch<fmax; ch++) {

int chx = sub(mul(ch,ch), 1);

sh = sqrts[chx];

P[0][0] = sub(0, ch);

P[0][WDIM] = sub(0, sh);

P[1][1] = Prime-1;

P[WDIM][0] = sh;

P[WDIM][WDIM] = ch;

if(!check_order(mmul(P, R), dual ? rotations : S3)) continue;

if(dual) R = mmul(P, R);

return 0;

}

}

}

}

return 2;

int cnt = 1;

matrix Po = M;

while(Po != Id) Po = mmul(Po, M), cnt++;

return cnt;

int N = isize(matrices);

auto compute_transcript = [&] (int i, int j) {

vector<int> indices(N, -1);

vector<int> transcript;

vector<int> q;

int qty = 0;

auto visit = [&] (int id) {

transcript.push_back(indices[id]);

if(indices[id] == -1) {

indices[id] = isize(q);

q.push_back(id);

qty++;

}

};

visit(0);

for(int x=0; x<isize(q); x++) {

int at = q[x];

visit(gmul(at, i));

visit(gmul(at, j));

}

transcript.push_back(qty);

return transcript;

};

if(debug_field) println(hlog, "looking for alternate solutions");

auto orig_transcript = compute_transcript(1, S7);

set<vector<int>> transcripts_seen;

transcripts_seen.insert(orig_transcript);

set<int> conjugacy_classes;

vector<int> cc;

for(int i=0; i<N; i++) conjugacy_classes.insert(i);

for(int i=0; i<N; i++) {

if(!conjugacy_classes.count(i)) continue;

vector<int> removals;

for(int j=0; j<N; j++) {

int c = gmul(inverses[j], gmul(i, j));

if(c > i) removals.push_back(c);

}

for(auto r: removals) conjugacy_classes.erase(r);

cc.push_back(i);

}

if(debug_field) println(hlog, "conjugacy_classes = ", cc);

vector<triplet_info> tinf;

triplet_info ti;

ti.i = 1; ti.j = S7; ti.size = orig_transcript.back();

tinf.push_back(ti);

for(int i: conjugacy_classes) if(gorder(i) == S7) {

if(debug_field) println(hlog, "checking i=", i);

for(int j=1; j<N; j++) if(gorder(j) == 2 && gorder(gmul(i, j)) == S3) {

auto t = compute_transcript(i, j);

if(!transcripts_seen.count(t)) {

transcripts_seen.insert(t);

triplet_info ti;

ti.i = i; ti.j = j; ti.size = t.back();

tinf.push_back(ti);

}

}

}

if(debug_field) println(hlog, "solutions found = ", isize(transcripts_seen));

return tinf;

if(WDIM == 3) return;

indenter_finish ibuild(debug_geometry, "fpattern::build");

for(int i=0; i<isize(qpaths); i++) {

matrix M = strtomatrix(qpaths[i]);

qcoords.push_back(M);

printf("Solved %s as matrix of order %d\n", qpaths[i].c_str(), order(M));

}

matcode.clear(); matrices.clear();

add(Id);

if(isize(matrices) != local_group) { printf("Error: rotation crash #1 (%d)\n", isize(matrices)); exit(1); }

connections.clear();

for(int i=0; i<(int)matrices.size(); i++) {

matrix M = matrices[i];

matrix PM = mmul(P, M);

add(PM);

if(isize(matrices) % local_group) { printf("Error: rotation crash (%d)\n", isize(matrices)); exit(1); }

if(!matcode.count(PM)) { printf("Error: not marked\n"); exit(1); }

connections.push_back(matcode[PM]);

}

if(debug_field) println(hlog, "Computing inverses...");

int N = isize(matrices);

if(debug_field) println(hlog, "Number of heptagons: ", N);

if(WDIM == 3) return;

rrf.resize(N); rrf[0] = S7-1;

for(int i=0; i<N; i++)

rrf[btspin(i,1)] = btspin(rrf[i], 1),

rrf[connections[i]] = connections[rrf[i]];

rpf.resize(N); rpf[0] = S7;

for(int i=0; i<N; i++)

rpf[btspin(i,1)] = btspin(rpf[i], 1),

rpf[connections[i]] = connections[rpf[i]];

inverses.resize(N);

inverses[0] = 0;

for(int i=0; i<N; i++) // inverses[i] = gpow(i, N-1);

inverses[btspin(i,1)] = rrf[inverses[i]], // btspin(inverses[i],6),

inverses[connections[i]] = rpf[inverses[i]];

int errs = 0;

for(int i=0; i<N; i++) if(gmul(i, inverses[i])) errs++;

if(errs) printf("errs = %d\n", errs);

if(0) for(int i=0; i<isize(matrices); i++) {

printf("%5d/%4d", connections[i], inverses[i]);

if(i%S7 == S7-1) printf("\n");

}

if(debug_field) println(hlog, "triplet_id = ", triplet_id, " N = ", N);

if(triplet_id) {

auto triplets = find_triplets();

if(triplet_id >= 0 && triplet_id < isize(triplets)) {

auto ti = triplets[triplet_id];

R = matrices[ti.i];

P = matrices[ti.j];

dynamicval<int> t(triplet_id, 0);

build();

if(debug_field) println(hlog, "triplet built successfully");

return;

}

}

if(debug_field) println(hlog, "Built.");

if(!a.second) return dists[a.first];

int m = MAXDIST;

int ma = dists[a.first];

int mb = dists[connections[btspin(a.first, 3)]];

int mc = dists[connections[btspin(a.first, 4)]];

m = min(m, 1 + ma);

m = min(m, 1 + mb);

m = min(m, 1 + mc);

if(m <= 2 && ma+mb+mc <= m*3-2) return m-1; // special case

m = min(m, 2 + dists[connections[btspin(a.first, 2)]]);

m = min(m, 2 + dists[connections[btspin(a.first, 5)]]);

m = min(m, 2 + dists[connections[btspin(connections[btspin(a.first, 3)], 5)]]);

return m;

if(a.first == b.first) return a.second == b.second ? 0 : 1;

if(b.first) a.first = gmul(a.first, inverses[b.first]), b.first = 0;

return getdist(a, b.second ? disthex : disthep);

int N = isize(matrices);

dists.resize(N);

for(int i=0; i<N; i++) dists[i] = MAXDIST-1;

int maxd = 0;

for(int i=0; i<MAXDIST; i++) while(!indist[i].empty()) {

int at = indist[i].back();

indist[i].pop_back();

if(dists[at] <= i) continue;

maxd = i;

dists[at] = i;

int lg = MWDIM == 4 ? local_group : S7;

for(int q=0; q<lg; q++) {

dists[at] = i;

if(WDIM == 3)

indist[i+1].push_back(gmul(at, local_group));

else if(PURE) // todo-variation: PURE here?

indist[i+1].push_back(connections[at]);

else {

indist[i+2].push_back(connections[at]);

indist[i+3].push_back(connections[btspin(connections[at], 2)]);

}

at = groupspin(at, 1, lg);

}

}

return maxd;

indenter_finish dif(debug_geometry, "fpattern::analyze");

if(MWDIM == 4) {

/* we need to compute inverses */

int N = isize(matrices);