// Hyperbolic Rogue -- regular honeycombs // Copyright (C) 2011-2019 Zeno Rogue, see 'hyper.cpp' for details /** \file reg3.cpp * \brief regular honeycombs * * works with spherical and hyperbolic ones -- Euclidean cubic tiling implemented in euclid.cpp * includes non-quotient spaces as well as field quotient and elliptic spaces * hyperbolic honeycombs rely on bt:: to deal with floating point errors (just like archimedean) */ #include "hyper.h" namespace hr { #if MAXMDIM >= 4 /** \brief regular three-dimensional tessellations */ EX namespace reg3 { #if HDR inline short& altdist(heptagon *h) { return h->emeraldval; } #endif EX int extra_verification; EX bool ultra_mirror_on; EX bool ultra_mirror_in() { return (cgflags & qULTRA) && ultra_mirror_on; } EX bool in() { if(fake::in()) return FPIU(in()); return WDIM == 3 && !euclid && !bt::in() && !nonisotropic && !hybri && !kite::in(); } EX void compute_ultra() { cgi.ultra_mirror_part = .99; cgi.ultra_material_part = .99; cgi.ultra_mirrors.clear(); if(cgflags & qULTRA) { for(auto& v: cgi.vertices_only) { hyperpoint nei; for(int i=0; i<isize(cgi.cellshape); i++) for(int j=0; j<isize(cgi.cellshape[i]); j++) if(sqhypot_d(WDIM, cgi.cellshape[i][j]-v) < 1e-6) nei = cgi.cellshape[i][j?j-1:j+1]; transmatrix T = spintox(v); hyperpoint a = T * v; hyperpoint b = T * nei; ld f0 = 0.5; ld f1 = binsearch(0.5, 1, [&] (ld d) { hyperpoint c = lerp(b, a, d); if(debugflags & DF_GEOM) println(hlog, "d=", d, " c= ", c, " material = ", material(c)); return material(c) <= 0; }); cgi.ultra_material_part = f1; auto f = [&] (ld d) { hyperpoint c = lerp(b, a, d); c = normalize(c); return c[1] * c[1] + c[2] * c[2]; }; for(int it=0; it<100; it++) { ld fa = (f0*2+f1) / 3; ld fb = (f0*1+f1*2) / 3; if(debugflags & DF_GEOM) println(hlog, "f(", fa, ") = ", f(fa), " f(", fb, ") = ", f(fb)); if(f(fa) > f(fb)) f0 = fa; else f1 = fb; } cgi.ultra_mirror_part = f0; hyperpoint c = lerp(b, a, f0); c = normalize(c); c[1] = c[2] = 0; c = normalize(c); cgi.ultra_mirror_dist = hdist0(c); if(cgi.ultra_mirror_part >= 1-1e-6) continue; cgi.ultra_mirrors.push_back(rspintox(v) * xpush(cgi.ultra_mirror_dist*2) * MirrorX * spintox(v)); } } } EX void make_vertices_only() { auto& vertices_only = cgi.vertices_only; vertices_only.clear(); for(auto& v: cgi.cellshape) for(hyperpoint h: v) { bool found = false; for(hyperpoint h2: vertices_only) if(hdist(h, h2) < 1e-6) found = true; if(!found) vertices_only.push_back(h); } } EX void generate() { if(fake::in()) { fake::generate(); return; } int& loop = cgi.loop; int& face = cgi.face; auto& spins = cgi.spins; auto& cellshape = cgi.cellshape; auto& adjcheck = cgi.adjcheck; auto& dirs_adjacent = cgi.dirs_adjacent; int& mid = cgi.schmid; mid = 3; face = 3; if(S7 == 6) face = 4; if(S7 == 8) mid = 4; if(S7 == 12) face = 5; if(S7 == 20) mid = 5; /* icosahedron not implemented */ loop = ginf[geometry].tiling_name[5] - '0'; DEBB(DF_GEOM, ("face = ", face, " loop = ", loop, " S7 = ", S7)); ld angle_between_faces, hcrossf; /* frontal face direction */ hyperpoint h0, h1, h2, h3, h012, h013; if(1) { dynamicval<eGeometry> dg(geometry, gSphere); angle_between_faces = edge_of_triangle_with_angles(2*M_PI/mid, M_PI/face, M_PI/face); h0 = xtangent(1); h1 = cspin(0, 1, angle_between_faces) * h0; h2 = cspin(1, 2, 2*M_PI/face) * h1; h3 = cspin(1, 2, -2*M_PI/face) * h1; hcrossf = edge_of_triangle_with_angles(M_PI/2, M_PI/mid, M_PI/face); h012 = cspin(1, 2, M_PI/face) * cspin(0, 1, hcrossf) * h0; h013 = cspin(1, 2, -M_PI/face) * cspin(0, 1, hcrossf) * h0; } for(auto hx: {&h0, &h1, &h2, &h3, &h012, &h013}) (*hx)[3] = 0; ld klein_scale = binsearch(0, 10, [&] (ld d) { dynamicval<eGeometry> g(geometry, elliptic ? gCell120 : geometry); /* center of an edge */ hyperpoint u = C0 + (h012 + h013) * d / 2; if(material(u) <= 0) { println(hlog, "klein_scale = ", d, " bad"); return true; } u = normalize(u); hyperpoint h = C0 * face; for(int i=0; i<face; i++) h += d * (cspin(1, 2, M_PI*2*i/face) * h012); h = normalize(h); hyperpoint h2 = rspintox(h) * xpush0(2 * hdist0(h)); h2 = spintox(u) * h2; u = spintox(u) * u; h2 = gpushxto0(u) * h2; u = gpushxto0(u) * u; ld x = hypot(h2[1], h2[2]); ld y = h2[0]; ld loop2 = 360 / (90 + atan(y/x) / degree); println(hlog, "d=", d, " loop2= ", loop2); if(sphere) return loop2 < loop; return loop2 > loop; }); /* precise ideal vertex */ if(klein_scale > 1-1e-5 && klein_scale < 1+1e-5) klein_scale = 1; /* actual vertex */ hyperpoint v2 = C0 + klein_scale * h012; hyperpoint midface = Hypc; for(int i=0; i<face; i++) midface += cspin(1, 2, 2*i*M_PI/face) * v2; midface = normalize(midface); ld between_centers = 2 * hdist0(midface); DEBB(DF_GEOM, ("between_centers = ", between_centers)); if(S7 == 20) { spins[0] = Id; spins[1] = cspin(0, 1, angle_between_faces) * cspin(1, 2, M_PI); spins[2] = spins[1] * cspin(1, 2, -2 * M_PI/face) * spins[1]; spins[3] = spins[1] * cspin(1, 2, +2 * M_PI/face) * spins[1]; for(int a=4; a<10; a++) spins[a] = cspin(1, 2, 2*M_PI/face) * spins[a-3]; for(int a=S7/2; a<S7; a++) spins[a] = spins[a-S7/2] * cspin(0, 1, M_PI); } if(S7 == 12 || S7 == 8) { spins[0] = Id; spins[1] = cspin(0, 1, angle_between_faces) * cspin(1, 2, M_PI); for(int a=2; a<face+1; a++) spins[a] = cspin(1, 2, 2*M_PI*(a-1)/face) * spins[1]; for(int a=S7/2; a<S7; a++) spins[a] = cspin(0, 1, M_PI) * spins[a-S7/2]; if(S7 == 8) swap(spins[6], spins[7]); if(S7 == 12) swap(spins[8], spins[11]); if(S7 == 12) swap(spins[9], spins[10]); } if(S7 == 6) { spins[0] = Id; spins[1] = cspin(0, 1, angle_between_faces) * cspin(1, 2, M_PI); spins[2] = cspin(1, 2, M_PI/2) * spins[1]; for(int a=S7/2; a<S7; a++) spins[a] = spins[a-S7/2] * cspin(0, 1, M_PI); } if(S7 == 4) { spins[0] = Id; spins[1] = cspin(0, 1, angle_between_faces) * cspin(1, 2, M_PI); for(int a=2; a<face+1; a++) spins[a] = cspin(1, 2, 2*M_PI*(a-1)/face) * spins[1]; } cellshape.clear(); cellshape.resize(S7); for(int a=0; a<S7; a++) { for(int b=0; b<face; b++) cellshape[a].push_back(spins[a] * cspin(1, 2, 2*M_PI*b/face) * v2); } cgi.adjmoves[0] = cpush(0, between_centers) * cspin(0, 2, M_PI); for(int i=1; i<S7; i++) cgi.adjmoves[i] = spins[i] * cgi.adjmoves[0]; for(int a=0; a<S7; a++) DEBB(DF_GEOM, ("center of ", a, " is ", tC0(cgi.adjmoves[a]))); DEBB(DF_GEOM, ("doublemove = ", tC0(cgi.adjmoves[0] * cgi.adjmoves[0]))); adjcheck = hdist(tC0(cgi.adjmoves[0]), tC0(cgi.adjmoves[1])) * 1.0001; int numedges = 0; for(int a=0; a<S7; a++) for(int b=0; b<S7; b++) { dirs_adjacent[a][b] = a != b && hdist(tC0(cgi.adjmoves[a]), tC0(cgi.adjmoves[b])) < adjcheck; if(dirs_adjacent[a][b]) numedges++; } DEBB(DF_GEOM, ("numedges = ", numedges)); if(loop == 4) cgi.strafedist = adjcheck; else cgi.strafedist = hdist(cgi.adjmoves[0] * C0, cgi.adjmoves[1] * C0); if(stretch::applicable()) { transmatrix T = cspin(0, 2, 90 * degree); transmatrix iT = inverse(T); for(auto& v: cgi.adjmoves) v = T * v * iT; for(auto& vv: cellshape) for(auto& v: vv) v = T * v; } make_vertices_only(); compute_ultra(); for(int a=0; a<S7; a++) for(int b=0; b<S7; b++) if(cgi.dirs_adjacent[a][b]) for(int c=0; c<S7; c++) if(cgi.dirs_adjacent[a][c] && cgi.dirs_adjacent[b][c]) { transmatrix t = build_matrix(tC0(cgi.adjmoves[a]), tC0(cgi.adjmoves[b]), tC0(cgi.adjmoves[c]), C0); if(det(t) > 1e-3) cgi.next_dir[a][b] = c; } } void binary_rebase(heptagon *h, const transmatrix& V) { } void test(); #if HDR struct hrmap_quotient3 : hrmap { vector<heptagon*> allh; vector<vector<transmatrix>> tmatrices; vector<cell*> acells; transmatrix adj(heptagon *h, int d) override { return tmatrices[h->fieldval][d]; } heptagon *getOrigin() override { return allh[0]; } transmatrix relative_matrix(heptagon *h2, heptagon *h1, const hyperpoint& hint) override; void initialize(int cell_count); vector<cell*>& allcells() override { return acells; } vector<hyperpoint> get_vertices(cell* c) override { return cgi.vertices_only; } }; #endif void hrmap_quotient3::initialize(int cell_count) { allh.resize(cell_count); acells.clear(); tmatrices.resize(cell_count); for(int a=0; a<cell_count; a++) { allh[a] = init_heptagon(S7); allh[a]->c7 = newCell(S7, allh[a]); allh[a]->fieldval = a; acells.push_back(allh[a]->c7); } } transmatrix hrmap_quotient3::relative_matrix(heptagon *h2, heptagon *h1, const hyperpoint& hint) { if(h1 == h2) return Id; int d = hr::celldistance(h2->c7, h1->c7); for(int a=0; a<S7; a++) if(hr::celldistance(h1->move(a)->c7, h2->c7) < d) return adj(h1, a) * relative_matrix(h2, h1->move(a), hint); for(int a=0; a<S7; a++) println(hlog, "d=", d, " vs ", hr::celldistance(h1->move(a)->c7, h2->c7)); println(hlog, "error in hrmap_quotient3:::relative_matrix"); return Id; } #if CAP_CRYSTAL int encode_coord(const crystal::coord& co) { int c = 0; for(int i=0; i<4; i++) c |= ((co[i]>>1) & 3) << (2*i); return c; } EX crystal::coord decode_coord(int a) { crystal::coord co; for(int i=0; i<4; i++) co[i] = (a & 3) * 2, a >>= 2; return co; } struct hrmap_from_crystal : hrmap_quotient3 { hrmap_from_crystal() { initialize(256); if(1) { auto m = crystal::new_map(); dynamicval<hrmap*> cm(currentmap, m); for(int a=0; a<256; a++) { auto co = decode_coord(a); heptagon *h1 = get_heptagon_at(co); for(int d=0; d<8; d++) { int b = encode_coord(crystal::get_coord(h1->cmove(d))); allh[a]->c.connect(d, allh[b], h1->c.spin(d), false); tmatrices[a].push_back(crystal::get_adj(h1, d)); } } delete m; } } }; #endif struct hrmap_field3 : reg3::hrmap_quotient3 { fieldpattern::fpattern *f; hrmap_field3(fieldpattern::fpattern *ptr) { f = ptr; auto lgr = f->local_group; int N = isize(f->matrices) / lgr; initialize(N); vector<int> moveid(S7), movedir(lgr); for(int s=0; s<lgr; s++) for(int i=0; i<S7; i++) if(eqmatrix(f->fullv[s] * cgi.adjmoves[0], cgi.adjmoves[i])) moveid[i] = s; for(int s=0; s<lgr; s++) for(int i=0; i<S7; i++) if(hdist(tC0(inverse(f->fullv[s]) * cgi.adjmoves[0]), tC0(cgi.adjmoves[i])) < 1e-4) movedir[s] = i; for(int a=0; a<N; a++) { tmatrices[a].resize(S7); for(int b=0; b<S7; b++) { int k = lgr*a; k = f->matcode[ f->mmul(f->mmul(f->matrices[k], f->matrices[moveid[b]]), f->P) ]; for(int l=0; l<lgr; l++) if(f->gmul(k, l) % lgr == 0) { tmatrices[a][b] = cgi.adjmoves[b] * f->fullv[l]; allh[a]->c.connect(b, allh[k/lgr], movedir[l], false); } } } create_patterns(); } set<cellwalker> plane; void make_plane(cellwalker cw) { if(plane.count(cw)) return; plane.insert(cw); for(int i=0; i<S7; i++) if(cgi.dirs_adjacent[i][cw.spin]) make_plane(reg3::strafe(cw, i)); } void create_patterns() { DEBB(DF_GEOM, ("creating pattern = ", isize(allh))); // also, strafe needs currentmap dynamicval<hrmap*> c(currentmap, this); if(S7 == 12) { // Emerald in 534 cell *a = gamestart(); cell *b = a; for(cell *c: allcells()) if(bounded_celldistance(a, c) == 5) { b = c; break; } for(cell *c: allcells()) if(bounded_celldistance(a, c) > bounded_celldistance(b, c)) c->master->zebraval |= 1; // Vineyard in 534 b = (cellwalker(a, 0) + wstep + rev + wstep).at; for(cell *c: allcells()) if(bounded_celldistance(a, c) == bounded_celldistance(b, c)) c->master->zebraval |= 2; } if(S7 == 6 && ginf[geometry].vertex == 5) { // Emerald in 534 cell *a = gamestart(); for(cell *c: allcells()) if(bounded_celldistance(a, c) > 3) c->master->zebraval |= 1; // Vineyard in 435 make_plane(cellwalker(gamestart(), 0)); DEBB(DF_GEOM, ("plane size = ", isize(plane))); set<int> plane_indices; for(auto cw: plane) plane_indices.insert(cw.at->master->fieldval); int fN = isize(f->matrices); set<int> nwi; for(int i=0; i<fN; i++) { bool ok = true; for(auto o: plane_indices) { int j = f->gmul(i, o * f->local_group) / f->local_group; if(plane_indices.count(j)) ok = false; forCellEx(c1, allcells()[j]) if(plane_indices.count(c1->master->fieldval)) ok = false; } if(ok) nwi.insert(i); } int gpow = 0; for(int i: nwi) { int pw = 1; int at = i; while(true) { at = f->gmul(at, i); if(!nwi.count(at)) break; pw++; } if(pw == 4) gpow = i; } int u = 0; for(int a=0; a<5; a++) { for(int o: plane_indices) { int j = f->gmul(u, o * f->local_group) / f->local_group; allcells()[j]->master->zebraval |= 2; } u = f->gmul(u, gpow); } } } }; /** \brief homology cover of the Seifert-Weber space */ namespace seifert_weber { using crystal::coord; vector<coord> periods; int flip(int x) { return (x+6) % 12; } void build_reps() { // start_game(); for(int a=0; a<12; a++) for(int b=0; b<12; b++) if(cgi.dirs_adjacent[a][b]) for(int c=0; c<12; c++) if(cgi.dirs_adjacent[a][c] && cgi.dirs_adjacent[b][c]) { transmatrix t = build_matrix(tC0(cgi.adjmoves[a]), tC0(cgi.adjmoves[b]), tC0(cgi.adjmoves[c]), C0); if(det(t) > 0) cgi.next_dir[a][b] = c; } set<coord> boundaries; for(int a=0; a<12; a++) for(int b=0; b<12; b++) if(cgi.dirs_adjacent[a][b]) { coord res = crystal::c0; int sa = a, sb = b; do { // printf("%d ", sa); if(sa < 6) res[sa]++; else res[sa-6]--; sa = flip(sa); sb = flip(sb); swap(sa, sb); sb = cgi.next_dir[sa][sb]; // sb = next_dirsa][sb]; } while(a != sa || b != sb); // printf("\n"); if(res > crystal::c0) boundaries.insert(res); } periods.clear(); for(int index = 5; index >= 0; index--) { for(auto k: boundaries) println(hlog, k); DEBB(DF_GEOM, ("simplifying...")); for(auto by: boundaries) if(among(by[index], 1, -1)) { DEBB(DF_GEOM, ("simplifying by ", by)); periods.push_back(by); set<coord> nb; for(auto v: boundaries) if(v == by) ; else if(v[index] % by[index] == 0) nb.insert(v - by * (v[index] / by[index])); else println(hlog, "error"); boundaries = move(nb); break; } } } int get_rep(coord a) { a = a - periods[0] * (a[5] / periods[0][5]); a = a - periods[1] * (a[4] / periods[1][4]); a = a - periods[2] * (a[3] / periods[2][3]); for(int i=0; i<3; i++) a[i] = gmod(a[i], 5); return a[2] * 25 + a[1] * 5 + a[0]; } coord decode(int id) { coord res = crystal::c0; for(int a=0; a<3; a++) res[a] = id % 5, id /= 5; return res; } struct hrmap_singlecell : hrmap_quotient3 { hrmap_singlecell(ld angle) { initialize(1); tmatrices[0].resize(S7); for(int b=0; b<S7; b++) { allh[0]->c.connect(b, allh[0], (b+S7/2) % S7, false); transmatrix T = cgi.adjmoves[b]; hyperpoint p = tC0(T); tmatrices[0][b] = rspintox(p) * xpush(hdist0(p)) * cspin(2, 1, angle) * spintox(p); } } }; struct hrmap_seifert_cover : hrmap_quotient3 { hrmap_seifert_cover() { if(periods.empty()) build_reps(); initialize(125); for(int a=0; a<125; a++) { tmatrices[a].resize(12); for(int b=0; b<12; b++) { coord x = decode(a); if(b < 6) x[b]++; else x[b-6]--; int a1 = get_rep(x); allh[a]->c.connect(b, allh[a1], flip(b), false); transmatrix T = cgi.adjmoves[b]; hyperpoint p = tC0(T); tmatrices[a][b] = rspintox(p) * xpush(hdist0(p)) * cspin(2, 1, 108 * degree) * spintox(p); } } } }; } struct hrmap_reg3 : hrmap { heptagon *origin; hrmap *binary_map; hrmap_quotient3 *quotient_map; map<heptagon*, pair<heptagon*, transmatrix>> reg_gmatrix; map<heptagon*, vector<pair<heptagon*, transmatrix> > > altmap; vector<cell*> spherecells; vector<cell*>& allcells() override { if(sphere) return spherecells; return hrmap::allcells(); } hrmap_reg3() { origin = init_heptagon(S7); heptagon& h = *origin; h.s = hsOrigin; h.c7 = newCell(S7, origin); if(sphere) spherecells.push_back(h.c7); worst_error1 = 0, worst_error2 = 0; dynamicval<hrmap*> cr(currentmap, this); heptagon *alt = NULL; transmatrix T = Id; binary_map = nullptr; quotient_map = nullptr; #if CAP_FIELD #if CAP_CRYSTAL if(geometry == gSpace344) { quotient_map = new hrmap_from_crystal; } else #endif if(geometry == gSpace535) { quotient_map = new seifert_weber::hrmap_seifert_cover; } else if(hyperbolic) { quotient_map = new hrmap_field3(&currfp); } #endif h.zebraval = quotient_map ? quotient_map->allh[0]->zebraval : 0; #if CAP_BT if(hyperbolic) { dynamicval<eGeometry> g(geometry, gBinary3); bt::build_tmatrix(); alt = init_heptagon(S7); alt->s = hsOrigin; alt->alt = alt; binary_map = bt::new_alt_map(alt); T = xpush(.01241) * spin(1.4117) * xpush(0.1241) * cspin(0, 2, 1.1249) * xpush(0.07) * Id; } #endif reg_gmatrix[origin] = make_pair(alt, T); altmap[alt].emplace_back(origin, T); celllister cl(origin->c7, 4, 100000, NULL); for(cell *c: cl.lst) { hyperpoint h = tC0(relative_matrix(c->master, origin, C0)); cgi.close_distances[bucketer(h)] = cl.getdist(c); } } ld worst_error1, worst_error2; heptagon *getOrigin() override { return origin; } void fix_distances(heptagon *h, heptagon *h2) { vector<heptagon*> to_fix; auto fix_pair = [&] (heptagon *h, heptagon *h2) { if(!h2) return; if(h->distance > h2->distance+1) { h->distance = h2->distance + 1; to_fix.push_back(h); } else if(h2->distance > h->distance+1) { h2->distance = h->distance + 1; to_fix.push_back(h2); } if(h->alt && h->alt == h2->alt) { if(altdist(h) > altdist(h2) + 1) { altdist(h) = altdist(h2) + 1; to_fix.push_back(h); } else if (altdist(h2) > altdist(h) + 1) { altdist(h2) = altdist(h) + 1; to_fix.push_back(h2); } } }; if(!h2) to_fix = {h}; else fix_pair(h, h2); for(int i=0; i<isize(to_fix); i++) { h = to_fix[i]; for(int j=0; j<S7; j++) fix_pair(h, h->move(j)); } } #define DEB 0 heptagon *counterpart(heptagon *h) { return quotient_map->allh[h->fieldval]; } void verify_neighbors(heptagon *alt, int steps, const hyperpoint& hT) { ld err; for(auto& p2: altmap[alt]) if((err = intval(tC0(p2.second), hT)) < 1e-3) { println(hlog, "FAIL"); exit(3); } #if CAP_BT if(steps) { dynamicval<eGeometry> g(geometry, gBinary3); dynamicval<hrmap*> cm(currentmap, binary_map); for(int i=0; i<alt->type; i++) verify_neighbors(alt->cmove(i), steps-1, currentmap->iadj(alt, i) * hT); } #endif } heptagon *create_step(heptagon *parent, int d) override { auto& p1 = reg_gmatrix[parent]; if(DEB) println(hlog, "creating step ", parent, ":", d, ", at ", p1.first, tC0(p1.second)); heptagon *alt = p1.first; #if CAP_FIELD transmatrix T = p1.second * (quotient_map ? quotient_map->tmatrices[parent->fieldval][d] : cgi.adjmoves[d]); #else transmatrix T = p1.second * cgi.adjmoves[d]; #endif transmatrix T1 = T; #if CAP_BT if(hyperbolic) { dynamicval<eGeometry> g(geometry, gBinary3); dynamicval<hrmap*> cm(currentmap, binary_map); binary_map->virtualRebase(alt, T); } #endif fixmatrix(T); auto hT = tC0(T); bool hopf = stretch::applicable(); if(hopf) T = stretch::translate(hT); if(DEB) println(hlog, "searching at ", alt, ":", hT); if(DEB) for(auto& p2: altmap[alt]) println(hlog, "for ", tC0(p2.second), " intval is ", intval(tC0(p2.second), hT)); ld err; for(auto& p2: altmap[alt]) if((err = intval(tC0(p2.second), hT)) < 1e-3) { if(err > worst_error1) println(hlog, format("worst_error1 = %lg", double(worst_error1 = err))); // println(hlog, "YES found in ", isize(altmap[alt])); if(DEB) println(hlog, "-> found ", p2.first); int fb = 0; hyperpoint old = tC0(p1.second);; if(!hopf) T * (inverse(T1) * old); #if CAP_FIELD if(quotient_map) { p2.first->c.connect(counterpart(parent)->c.spin(d), parent, d, false); fix_distances(p2.first, parent); return p2.first; } #endif for(int d2=0; d2<S7; d2++) { hyperpoint back = p2.second * tC0(cgi.adjmoves[d2]); if((err = intval(back, old)) < 1e-3) { if(err > worst_error2) println(hlog, format("worst_error2 = %lg", double(worst_error2 = err))); if(p2.first->move(d2)) println(hlog, "error: repeated edge"); p2.first->c.connect(d2, parent, d, false); fix_distances(p2.first, parent); fb++; } } if(fb != 1) { println(hlog, "found fb = ", fb); println(hlog, old); for(int d2=0; d2<S7; d2++) { println(hlog, p2.second * tC0(cgi.adjmoves[d2]), " in distance ", intval(p2.second * tC0(cgi.adjmoves[d2]), old)); } parent->c.connect(d, parent, d, false); return parent; } return p2.first; } if(extra_verification) verify_neighbors(alt, extra_verification, hT); if(DEB) println(hlog, "-> not found"); int d2 = 0, fv = isize(reg_gmatrix); #if CAP_FIELD if(quotient_map) { auto cp = counterpart(parent); d2 = cp->c.spin(d); fv = cp->c.move(d)->fieldval; } #endif if(hopf) { hyperpoint old = tC0(p1.second); for(d2=0; d2<S7; d2++) { hyperpoint back = T * tC0(cgi.adjmoves[d2]); if((err = intval(back, old)) < 1e-3) break; } if(d2 == S7) { d2 = 0; println(hlog, "Hopf connection failed"); } println(hlog, "found d2 = ", d2); } heptagon *created = init_heptagon(S7); created->c7 = newCell(S7, created); if(sphere) spherecells.push_back(created->c7); #if CAP_FIELD if(quotient_map) { created->emeraldval = fv; created->zebraval = quotient_map->allh[fv]->zebraval; } else #endif created->zebraval = hrand(10); created->fieldval = fv; created->distance = parent->distance + 1; created->fiftyval = 9999; fixmatrix(T); reg_gmatrix[created] = make_pair(alt, T); altmap[alt].emplace_back(created, T); created->c.connect(d2, parent, d, false); return created; } ~hrmap_reg3() { #if CAP_BT if(binary_map) { dynamicval<eGeometry> g(geometry, gBinary3); delete binary_map; } #endif if(quotient_map) delete quotient_map; clearfrom(origin); } map<heptagon*, int> reducers; void link_alt(const cellwalker& hs) override { auto h = hs.at->master; altdist(h) = 0; if(h->alt->s != hsOrigin) reducers[h] = hs.spin; } void generateAlts(heptagon* h, int levs, bool link_cdata) override { if(reducers.count(h)) { heptspin hs(h, reducers[h]); reducers.erase(h); hs += wstep; hs += rev; altdist(hs.at) = altdist(h) - 1; hs.at->alt = h->alt; reducers[hs.at] = hs.spin; fix_distances(hs.at, NULL); } for(int i=0; i<S7; i++) { auto h2 = h->cmove(i); if(h2->alt == NULL) { h2->alt = h->alt; altdist(h2) = altdist(h) + 1; fix_distances(h2, NULL); } } } transmatrix adj(heptagon *h, int d) override { #if CAP_FIELD if(quotient_map) return quotient_map->adj(h, d); else #endif return relative_matrix(h->cmove(d), h, C0); } transmatrix relative_matrix(heptagon *h2, heptagon *h1, const hyperpoint& hint) override { auto p1 = reg_gmatrix[h1]; auto p2 = reg_gmatrix[h2]; transmatrix T = Id; #if CAP_BT if(hyperbolic) { dynamicval<eGeometry> g(geometry, gBinary3); dynamicval<hrmap*> cm(currentmap, binary_map); T = binary_map->relative_matrix(p2.first, p1.first, hint); } #endif T = inverse(p1.second) * T * p2.second; if(elliptic && T[LDIM][LDIM] < 0) T = centralsym * T; return T; } vector<hyperpoint> get_vertices(cell* c) override { return cgi.vertices_only; } }; struct hrmap_reg3_rule : hrmap { heptagon *origin; reg3::hrmap_quotient3 *quotient_map; reg3::hrmap_quotient3 *emerald_map; fieldpattern::fpattern fp; vector<int> root; string other; vector<short> children; vector<int> otherpos; void load_ruleset(string fname) { string buf; #if ISANDROID || ISIOS buf = get_asset(fname); #else FILE *f = fopen(fname.c_str(), "rb"); if(!f) f = fopen((rsrcdir + fname).c_str(), "rb"); buf.resize(1000000); int qty = fread(&buf[0], 1, 1000000, f); buf.resize(qty); fclose(f); #endif shstream ins(decompress_string(buf)); dynamicval<bool> q(fieldpattern::use_quotient_fp, true); hread_fpattern(ins, fp); hread(ins, root); hread(ins, children); hread(ins, other); } /** \brief address = (fieldvalue, state) */ typedef pair<int, int> address; /** nles[x] lists the addresses from which we can reach address x * without ever ending in the starting point */ map<address, set<address>> nonlooping_earlier_states; vector<vector<int>> possible_states; void find_mappings() { auto &nles = nonlooping_earlier_states; nles.clear(); vector<address> bfs; int qty = isize(quotient_map->allh); if(geometry == gSpace535) qty = 1; for(int i=0; i<qty; i++) bfs.emplace_back(i, root[i]); auto mov = [&] (int fv, int d) { if(geometry == gSpace535) return 0; return quotient_map->allh[fv]->move(d)->fieldval; }; int qstate = isize(children) / S7; DEBB(DF_GEOM, ("qstate = ", qstate)); for(int i=0; i<isize(bfs); i++) { address last = bfs[i]; int state = last.second; int fv = last.first; for(int d=0; d<S7; d++) { int nstate = children[state*S7+d]; if(nstate < -1) nstate += (1<<16); if(nstate >= 0) { address next = {mov(fv, d), nstate}; if(!nles.count(next)) bfs.push_back(next); nles[next].insert(last); } } } vector<int> q(qstate, 0); for(auto p: bfs) q[p.second]++; vector<int> q2(isize(quotient_map->allh)+1, 0); for(auto p: q) q2[p]++; DEBB(DF_GEOM, ("q2 = ", q2)); bfs = {}; for(int i=0; i<qty; i++) bfs.emplace_back(i, root[i]); for(int i=0; i<isize(bfs); i++) { address last = bfs[i]; int state = last.second; int fv = last.first; for(int d=0; d<S7; d++) { int nstate = children[state*S7+d]; if(nstate < -1) nstate += (1<<16); if(nstate >= 0) { address next = {mov(fv, d), nstate}; if(!nles.count(next)) continue; int c = isize(nles[next]); nles[next].erase(last); if(nles[next].empty() && c) { nles.erase(next); bfs.push_back(next); } } } } DEBB(DF_GEOM, ("removed cases = ", isize(bfs))); // just the number of FV's int pstable = 0; for(auto& p: nonlooping_earlier_states) pstable = max(pstable, p.first.first+1); println(hlog, "pstable size = ", pstable, " (states: ", qstate, ")"); possible_states.resize(pstable); for(auto& p: nonlooping_earlier_states) possible_states[p.first.first].push_back(p.first.second); } hrmap_reg3_rule() : fp(0) { load_ruleset(get_rule_filename()); origin = init_heptagon(S7); heptagon& h = *origin; h.s = hsOrigin; h.fiftyval = root[0]; h.c7 = newCell(S7, origin); int opos = 0; for(int c: children) { if(c < -1) c += (1<<16); if(c >= 0) otherpos.push_back(-1); else { otherpos.push_back(opos); while(other[opos] != ',') opos++; opos++; } } quotient_map = nullptr; if(geometry == gSpace535) quotient_map = new seifert_weber::hrmap_seifert_cover(); #if CAP_CRYSTAL else if(geometry == gSpace344) quotient_map = new hrmap_from_crystal; #endif else quotient_map = new hrmap_field3(&fp); if(geometry == gSpace535) emerald_map = new seifert_weber::hrmap_seifert_cover(); #if CAP_CRYSTAL else if(geometry == gSpace344) emerald_map = new hrmap_from_crystal; #endif else emerald_map = new hrmap_field3(&currfp); h.emeraldval = 0; find_mappings(); } heptagon *getOrigin() override { return origin; } #define DEB 0 heptagon *counterpart(heptagon *h) { return quotient_map->allh[h->fieldval]; } vector<short> evmemo; void find_emeraldval(heptagon *target, heptagon *parent, int d) { if(geometry == gSpace535) { target->emeraldval = target->fieldval; target->zebraval = 0; return; } generate_cellrotations(); auto& cr = cgi.cellrotations; if(evmemo.empty()) { println(hlog, "starting"); map<int, int> matrix_hashtable; auto matrix_hash = [] (const transmatrix& M) { return bucketer(M[0][0]) + bucketer(M[0][1]) * 71 + bucketer(M[0][2]) * 113 + bucketer(M[1][0]) * 1301 + bucketer(M[1][1]) * 1703 + bucketer(M[1][2]) * 17031 + bucketer(M[2][2]) * 2307 + bucketer(M[2][0]) * 2311 + bucketer(M[2][1]) * 10311; }; for(int i=0; i<isize(cr); i++) matrix_hashtable[matrix_hash(cr[i].M)] = cr[i].inverse_id; println(hlog, "ids size = ", isize(matrix_hashtable)); for(int eid=0; eid<isize(emerald_map->allh); eid++) for(int k0=0; k0<isize(cr); k0++) for(int fv=0; fv<isize(quotient_map->allh); fv++) { for(int d=0; d<S7; d++) { int ed = cr[k0].mapping[d]; auto cpart = emerald_map->allh[eid]; int eid1 = emerald_map->allh[eid]->move(ed)->fieldval; const transmatrix& X = cr[cr[k0].inverse_id].M; transmatrix U = quotient_map->iadj(quotient_map->allh[fv], d) * X * emerald_map->adj(cpart, ed); int k1 = matrix_hashtable[matrix_hash(U)]; /* for(int ik1=0; ik1<isize(cr); ik1++) { auto& mX1 = cr[ik1].M; if(eqmatrix(mX1, U)) k1 = cr[ik1].inverse_id; } */ evmemo.push_back(eid1 * isize(cr) + k1); } } println(hlog, "generated ", isize(evmemo)); } int memo_id = parent->emeraldval; memo_id = memo_id * isize(quotient_map->allh) + parent->fieldval; memo_id = memo_id * S7 + d; target->emeraldval = evmemo[memo_id]; target->zebraval = emerald_map->allh[target->emeraldval / isize(cr)]->zebraval; } heptagon *create_step(heptagon *parent, int d) override { int id = parent->fiftyval; if(id < 0) id += (1<<16); auto cp = counterpart(parent); int d2 = cp->c.spin(d); int fv = cp->c.move(d)->fieldval; // indenter ind(2); heptagon *res = nullptr; int id1 = children[S7*id+d]; int pos = otherpos[S7*id+d]; if(id1 < -1) id1 += (1<<16); if(id1 == -1 && false) { int kk = pos; string s; while(other[kk] != ',') s += other[kk++]; println(hlog, "id=", id, " d=", d, " d2=", d2, " id1=", id1, " pos=", pos, " s = ", s); } if(id1 != -1) { res = init_heptagon(S7); if(parent->c7) res->c7 = newCell(S7, res); res->fieldval = fv; res->distance = parent->distance + 1; res->fiftyval = id1; find_emeraldval(res, parent, d); // res->c.connect(d2, parent, d, false); } else if(other[pos] == ('A' + d) && other[pos+1] == ',') { res = init_heptagon(S7); res->alt = parent->alt; res->fieldval = fv; res->distance = parent->distance - 1; vector<int> possible; int pfv = parent->fieldval; if(geometry == gSpace535) pfv = 0; for(auto s: nonlooping_earlier_states[address{pfv, id}]) possible.push_back(s.second); id1 = hrand_elt(possible, 0); res->fiftyval = id1; find_emeraldval(res, parent, d); } else { heptagon *at = parent; while(other[pos] != ',') { int dir = (other[pos++] & 31) - 1; // println(hlog, "from ", at, " go dir ", dir); at = at->cmove(dir); } res = at; } if(!res) throw hr_exception("res missing"); if(res->move(d2)) println(hlog, "res conflict"); res->c.connect(d2, parent, d, false); return res; } ~hrmap_reg3_rule() { if(quotient_map) delete quotient_map; clearfrom(origin); } transmatrix adj(heptagon *h, int d) override { return quotient_map->adj(h, d); } transmatrix relative_matrix(heptagon *h2, heptagon *h1, const hyperpoint& hint) override { return relative_matrix_recursive(h2, h1); } vector<hyperpoint> get_vertices(cell* c) override { return cgi.vertices_only; } }; struct hrmap_reg3_rule_alt : hrmap { heptagon *origin; hrmap_reg3_rule_alt(heptagon *o) { origin = o; } }; EX hrmap *new_alt_map(heptagon *o) { return new hrmap_reg3_rule_alt(o); } EX void link_structures(heptagon *h, heptagon *alt, hstate firststate) { auto cm = (hrmap_reg3_rule*) currentmap; alt->fieldval = h->fieldval; if(geometry == gSpace535) alt->fieldval = 0; if(firststate == hsOrigin) { alt->fiftyval = cm->root[alt->fieldval]; return; } vector<int>& choices = cm->possible_states[alt->fieldval]; vector<int> choices2; for(auto c: choices) { bool ok = true; for(int d=0; d<12; d++) if(h->cmove(d)->distance < h->distance) if(cm->children[S7*c+d] == -1) ok = false; if(ok) choices2.push_back(c); } alt->fiftyval = hrand_elt(choices2, -1); } EX bool reg3_rule_available = true; EX string other_rule = ""; EX string get_rule_filename() { if(other_rule != "") return other_rule; switch(geometry) { case gSpace336: return "honeycomb-rules-336.dat"; case gSpace344: return "honeycomb-rules-344.dat"; // case gSpace345: return "honeycomb-rules-345.dat"; case gSpace353: return "honeycomb-rules-353.dat"; case gSpace354: return "honeycomb-rules-354.dat"; // case gSpace355: return "honeycomb-rules-355.dat"; case gSpace435: return "honeycomb-rules-435.dat"; case gSpace436: return "honeycomb-rules-436.dat"; case gSpace534: return "honeycomb-rules-534.dat"; case gSpace535: return "honeycomb-rules-535.dat"; case gSpace536: return "honeycomb-rules-536.dat"; default: return ""; } } EX bool in_rule() { return reg3_rule_available && get_rule_filename() != ""; } EX int rule_get_root(int i) { return ((hrmap_reg3_rule*)currentmap)->root[i]; } EX const vector<short>& rule_get_children() { return ((hrmap_reg3_rule*)currentmap)->children; } EX hrmap* new_map() { if(geometry == gSeifertCover) return new seifert_weber::hrmap_seifert_cover; if(geometry == gSeifertWeber) return new seifert_weber::hrmap_singlecell(108*degree); if(geometry == gHomologySphere) return new seifert_weber::hrmap_singlecell(36*degree); if(quotient && !sphere) return new hrmap_field3(&currfp); if(in_rule()) return new hrmap_reg3_rule; return new hrmap_reg3; } hrmap_reg3* regmap() { return ((hrmap_reg3*) currentmap); } EX int quotient_count() { return isize(regmap()->quotient_map->allh); } /** This is a generalization of hyperbolic_celldistance in expansion.cpp to three dimensions. It still assumes that there are at most 4 cells around every edge, and that distances from the origin are known, so it works only in {5,3,4}. */ int celldistance_534(cell *c1, cell *c2) { int d1 = celldist(c1); int d2 = celldist(c2); vector<cell*> s1 = {c1}; vector<cell*> s2 = {c2}; int best = 99999999; int d0 = 0; auto go_nearer = [&] (vector<cell*>& v, int& d) { vector<cell*> w; for(cell *c: v) forCellEx(c1, c) if(celldist(c1) < d) w.push_back(c1); sort(w.begin(), w.end()); d--; d0++; auto last = std::unique(w.begin(), w.end()); w.erase(last, w.end()); v = w; }; while(d0 < best) { for(cell *a1: s1) for(cell *a2: s2) { if(a1 == a2) best = min(best, d0); else if(isNeighbor(a1, a2)) best = min(best, d0+1); } if(d1 == 0 && d2 == 0) break; if(d1 >= d2) go_nearer(s1, d1); if(d1 < d2) go_nearer(s2, d2); } return best; } EX int celldistance(cell *c1, cell *c2) { if(c1 == c2) return 0; if(c1 == currentmap->gamestart()) return c2->master->distance; if(c2 == currentmap->gamestart()) return c1->master->distance; if(geometry == gSpace534) return celldistance_534(c1, c2); auto r = regmap(); hyperpoint h = tC0(r->relative_matrix(c1->master, c2->master, C0)); int b = bucketer(h); if(cgi.close_distances.count(b)) return cgi.close_distances[b]; if(in_rule()) return clueless_celldistance(c1, c2); dynamicval<eGeometry> g(geometry, gBinary3); #if CAP_BT return 20 + bt::celldistance3(r->reg_gmatrix[c1->master].first, r->reg_gmatrix[c2->master].first); #else return 20; #endif } EX bool pseudohept(cell *c) { auto m = regmap(); if(cgflags & qSINGLE) return true; if(fake::in()) return FPIU(reg3::pseudohept(c)); if(sphere) { hyperpoint h = tC0(m->relative_matrix(c->master, regmap()->origin, C0)); if(S7 == 12) { hyperpoint h1 = cspin(0, 1, atan2(16, 69) + M_PI/4) * h; for(int i=0; i<4; i++) if(abs(abs(h1[i]) - .5) > .01) return false; return true; } if(S7 == 8) return h[3] >= .99 || h[3] <= -.99 || abs(h[3]) < .01; if(cgi.loop == 3 && cgi.face == 3 && S7 == 4) return c == m->gamestart(); if(cgi.loop == 4 && cgi.face == 3) return abs(h[3]) > .9; if(cgi.loop == 3 && cgi.face == 4) return abs(h[3]) > .9; if(cgi.loop == 5 && cgi.face == 3) return abs(h[3]) > .99 || abs(h[0]) > .99 || abs(h[1]) > .99 || abs(h[2]) > .99; } // chessboard pattern in 534 if(geometry == gField534) return hr::celldistance(c, currentmap->gamestart()) & 1; if(geometry == gCrystal344 || geometry == gCrystal534 || geometry == gSeifertCover) return false; if(quotient) return false; /* added */ auto mr = dynamic_cast<hrmap_reg3_rule*> (currentmap); if(mr) { if(geometry == gSpace535) return c->master->fieldval % 31 == 0; return c->master->fieldval == 0; } if(m && hyperbolic) { heptagon *h = m->reg_gmatrix[c->master].first; return (h->zebraval == 1) && (h->distance & 1); } return false; } EX void generate_cellrotations() { auto &cr = cgi.cellrotations; if(isize(cr)) return; for(int a=0; a<S7; a++) for(int b=0; b<S7; b++) for(int c=0; c<S7; c++) { transmatrix T = build_matrix(cgi.adjmoves[a]*C0, cgi.adjmoves[b]*C0, cgi.adjmoves[c]*C0, C0); if(abs(det(T)) < 0.001) continue; transmatrix U = build_matrix(cgi.adjmoves[0]*C0, cgi.adjmoves[1]*C0, cgi.adjmoves[2]*C0, C0); transmatrix S = U * inverse(T); if(abs(det(S) - 1) > 0.01) continue; vector<int> perm(S7); for(int x=0; x<S7; x++) perm[x] = -1; for(int x=0; x<S7; x++) for(int y=0; y<S7; y++) if(hdist(S * cgi.adjmoves[x] * C0, cgi.adjmoves[y] * C0) < .1) perm[x] = y; bool bad = false; for(int x=0; x<S7; x++) if(perm[x] == -1) bad = true; if(bad) continue; cr.emplace_back(geometry_information::cellrotation_t{S, perm, 0}); } int rots = isize(cr); for(int i=0; i<rots; i++) for(int j=0; j<rots; j++) if(cr[i].mapping[cr[j].mapping[0]] == 0 && cr[i].mapping[cr[j].mapping[1]] == 1 && cr[i].mapping[cr[j].mapping[2]] == 2) cr[i].inverse_id = j; } #endif #if 0 /* More precise, but very slow distance. Not used/optimized for now */ ld adistance(cell *c) { hyperpoint h = tC0(regmap()->reg_gmatrix[c->master].second); h = bt::deparabolic3(h); return regmap()->reg_gmatrix[c->master].first->distance * log(2) - h[0]; } map<pair<cell*, cell*>, int> memo; bool cdd; int celldistance(cell *c1, cell *c2) { if(memo.count(make_pair(c1, c2))) return memo[make_pair(c1, c2)]; if(c1 == c2) return 0; vector<cell*> v[2]; v[0].push_back(c1); v[1].push_back(c2); int steps = 0; map<cell*, int> visited; visited[c1] = 1; visited[c2] = 2; while(true) { if(cdd) { println(hlog, "state ", steps, "/",isize(v[0]), "/", isize(v[1])); println(hlog, " A: ", v[0]); println(hlog, " B: ", v[1]); } for(int i: {0,1}) { vector<cell*> new_v; for(cell *c: v[i]) forCellCM(cn, c) if(adistance(cn) < adistance(c)) { auto &vi = visited[cn]; if((vi&3) == 0) { vi = 4 * (steps+1); vi |= (1<<i); new_v.push_back(cn); } else if((vi&3) == 2-i) { vector<pair<cell*, int>> ca1, ca2; int b1 = 4*steps-4; int b2 = ((vi>>2)<<2) - 4; for(auto p: visited) { if(cdd) println(hlog, p); int ps = p.second & 3; if(ps == 1+i && p.second >= b1) ca1.emplace_back(p.first, p.second/4); if(ps == 2-i && p.second >= b2 && p.second <= b2+8) ca2.emplace_back(p.first, p.second/4); } int bound = 1<<16; for(auto p1: ca1) for(auto p2: ca2) { hyperpoint h = tC0(relative_matrix(p1.first->master, p2.first->master)); int b = bucketer(h); if(close_distances.count(b)) { int d = close_distances[b] + p1.second + p2.second; if(cdd) println(hlog, "candidate: close=", close_distances[b], p1, p2, "; h = ", h); if(d < bound) bound = d; } else if(cdd) println(hlog, "bucket missing"); } return memo[make_pair(c1, c2)] = bound; return bound; } } v[i] = std::move(new_v); } steps++; } } cellwalker target; int tsteps; int dist_alt(cell *c) { if(!target.at) { target = cellwalker(currentmap->gamestart(), 0); tsteps = 0; for(int i=0; i<30; i++) target += wstep, target += rev, tsteps++; } if(specialland == laCamelot) return reg3::celldistance(c, target.at); else { int d = reg3::celldistance(c, target.at) - tsteps; if(d < 10) target += wstep, target += rev, tsteps++; return d; } } #endif // Construct a cellwalker in direction j from cw.at, such that its direction is as close // as possible to cw.spin. Assume that j and cw.spin are adjacent #if MAXMDIM >= 4 EX cellwalker strafe(cellwalker cw, int j) { hyperpoint hfront = tC0(cgi.adjmoves[cw.spin]); cw.at->cmove(j); transmatrix T = currentmap->adj(cw.at, j); for(int i=0; i<S7; i++) if(i != cw.at->c.spin(j)) if(hdist(hfront, T * tC0(cgi.adjmoves[i])) < cgi.strafedist + .01) return cellwalker(cw.at->cmove(j), i); println(hlog, "incorrect strafe"); exit(1); } EX int matrix_order(const transmatrix A) { transmatrix T = A; int res = 1; while(!eqmatrix(T, Id)) { res++; T = T * A; } return res; } EX void generate_fulls() { reg3::generate_cellrotations(); auto cons = [&] (int i0, int i1, int i2) { transmatrix T = build_matrix(cgi.adjmoves[ 0]*C0, cgi.adjmoves[ 1]*C0, cgi.adjmoves[ 2]*C0, C0); transmatrix U = build_matrix(cgi.adjmoves[i0]*C0, cgi.adjmoves[i1]*C0, cgi.adjmoves[i2]*C0, C0); return U * inverse(T); }; cgi.full_P = cgi.adjmoves[0]; cgi.full_R = S7 == 8 ? cons(1, 7, 0) : S7 == 20 ? cons(1,2,6) : cons(1, 2, 0); cgi.full_X = S7 == 8 ? cons(1, 0, 6) : S7 == 6 ? cons(1, 0, 5) : S7 == 20 ? cons(1,0,7) : cons(1, 0, cgi.face); cgi.xp_order = matrix_order(cgi.full_X * cgi.full_P); cgi.r_order = matrix_order(cgi.full_R); cgi.rx_order = matrix_order(cgi.full_R * cgi.full_X); println(hlog, "orders = ", tie(cgi.rx_order, cgi.r_order, cgi.xp_order)); } EX void construct_relations() { auto& rels = cgi.rels; if(!rels.empty()) return; rels.clear(); reg3::generate_cellrotations(); reg3::generate_fulls(); vector<transmatrix> all; vector<string> formulas; formulas.push_back(""); all.push_back(Id); hyperpoint v = cgi.cellshape[0][0]; auto add = [&] (transmatrix T) { for(int i=0; i<isize(all); i++) if(eqmatrix(all[i], T)) return i; int S = isize(all); all.push_back(T); return S; }; println(hlog, cgi.cellshape); println(hlog, "cellshape = ", isize(cgi.cellshape)); bool ok = true; int last_i = -1; for(auto& v: cgi.cellshape) for(hyperpoint h: v) { int i = 0, j = 0; for(auto& uv: cgi.cellshape) for(hyperpoint u: uv) { if(hdist(h, cgi.full_X*u) < 5e-2) i++; if(hdist(h, cgi.full_R*u) < 5e-2) j++; } if(last_i == -1) last_i = i; if(i != j || i != last_i) ok = false; } if(!ok) { println(hlog, "something wrong"); exit(1); } add(Id); auto work = [&] (transmatrix T, int p, char c) { if(hdist0(tC0(T)) > 5) return; for(auto& hv: cgi.cellshape) for(hyperpoint h: hv) if(hdist(T * h, v) < 1e-4) goto ok; return; ok: int id = add(T); // println(hlog, p, " x ", (s0+c), " = ", id); if(id >= isize(formulas)) formulas.push_back(formulas[p] + c); else if(id == 0) println(hlog, "reached identity: ", formulas[p]+c); else if(formulas[p][0] != formulas[id][0]) rels.emplace_back(formulas[p] + c, formulas[id]); }; for(int i=0; i<isize(all); i++) { transmatrix T = all[i]; work(T * cgi.full_R, i, 'R'); work(T * cgi.full_X, i, 'X'); work(T * cgi.full_P, i, 'P'); } } EX } #endif #if MAXMDIM == 3 EX namespace reg3 { EX bool in() { return false; } EX bool in_rule() { return false; } EX } #endif }