// 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 { EX hyperpoint final_coords(hyperpoint h) { if(sn::in() || !bt::in()) return ultra_normalize(h); #if CAP_BT if(bt::in() && !prod) return bt::bt_to_minkowski(h); #endif return h; } void subcellshape::compute_common() { reg3::make_vertices_only(vertices_only, faces); faces_local = faces; for(auto& face: faces_local) for(auto& v: face) v = from_cellcenter * final_coords(v); vertices_only_local = vertices_only; for(auto& v: vertices_only_local) v = from_cellcenter * final_coords(v); int N = isize(faces); dirdist.resize(N); for(int i=0; i cface; for(auto& v: faces[i]) cface.insert(bucketer(v)); for(int j=0; j 0) next_dir[a][b] = c; } } void subcellshape::compute_hept() { cellcenter = C0; to_cellcenter = Id; from_cellcenter = Id; compute_common(); } EX namespace reg3 { EX void make_vertices_only(vector& vo, const vector>& csh) { vo.clear(); for(auto& v: csh) for(hyperpoint h: v) { bool found = false; for(hyperpoint h2: vo) if(hdist(h, h2) < 1e-6) found = true; if(!found) vo.push_back(h); } } EX } #if MAXMDIM >= 4 void subcellshape::compute_sub() { hyperpoint gres = Hypc; for(auto& face: faces) { hyperpoint res = Hypc; for(auto& vertex: face) res += vertex; face_centers.push_back(normalize(res)); gres += res; } cellcenter = normalize(gres); to_cellcenter = rgpushxto0(cellcenter); from_cellcenter = gpushxto0(cellcenter); compute_common(); } /** \brief regular three-dimensional tessellations */ EX namespace reg3 { EX int subcube_count = 1; EX flagtype coxeter_param = 0; const flagtype cox_othercell = 1; const flagtype cox_midedges = 2; const flagtype cox_vertices = 4; #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 cubes_reg3; EX bool in() { if(fake::in()) return FPIU(in()); if(geometry == gCubeTiling && (cubes_reg3 || !PURE)) return true; 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.heptshape->vertices_only) { hyperpoint nei; auto& faces = cgi.heptshape->faces; for(int i=0; ifaces[0][0][0]); auto vz = abs(cgi.heptshape->faces[0][0][3]); for(int x=1-sub; xfaces; for(auto& face: ss.faces) for(auto& v: face) { v[0] += vx * x; v[1] += vx * y; v[2] += vx * z; v[3] += vz * (sub-1); } } } EX void generate_coxeter(flagtype f) { auto& ssh = cgi.subshapes; for(auto& fac: cgi.heptshape->faces) { hyperpoint facectr = Hypc; vector ring; hyperpoint last = fac.back(); ring.push_back(last); for(hyperpoint h: fac) { if(f & cox_midedges) ring.push_back(mid(last, h)); ring.push_back(last = h); facectr += h; } facectr = normalize(facectr); hyperpoint fc2 = rspintox(facectr) * xpush0(2*hdist0(facectr)); if(f & cox_vertices) { for(int i=1; ifaces[0][0][0]); auto vz = abs(cgi.heptshape->faces[0][0][3]); auto step = hdist0(tC0(cgi.adjmoves[0])); array co; int s = bch ? 1 : 2; for(co[0]=-sub; co[0]<=sub; co[0]+=s) for(co[1]=-sub; co[1]<=sub; co[1]+=s) for(co[2]=-sub; co[2]<=sub; co[2]+=s) if(((co[0]^co[1]^1)&1) && ((co[0]^co[2]^1)&1)) { hyperpoint ctr = Hypc; ctr[3] = vz * sub; struct direction { hyperpoint dir; int limit; transmatrix mirror; void flip() { dir = -dir; limit = 200 - limit; } }; array di; int mi = 0; for(int i=0; i<3; i++) { ctr[i] += co[i] * vx; auto& dii = di[i]; if(co[i] >= 0) { dii.dir = ctangent(i, vx); dii.limit = sub - co[i]; dii.mirror = cpush(i, +step/2) * cmirror(i) * cpush(i, -step/2); } else { dii.dir = ctangent(i, -vx); dii.limit = co[i] + sub; dii.mirror = cpush(i, -step/2) * cmirror(i) * cpush(i, +step/2); } if(dii.limit == 0) mi++; } sort(di.begin(), di.end(), [] (direction& d1, direction& d2) { return d1.limit > d2.limit; }); cgi.subshapes.emplace_back(); auto &ss = cgi.subshapes.back(); auto pt0 = [&] (const array& x) { transmatrix M = Id; hyperpoint res = ctr; for(int i=0; i<3; i++) { ld xx = x[i]; if(xx > di[i].limit) xx = 2*di[i].limit-xx, M = di[i].mirror * M; res += di[i].dir * xx; } return normalize(M * res); }; auto pt = [&] (ld x, ld y, ld z) { if(sub == 1 || !bch || sphere) return pt0(make_array(x,y,z)); // Unfortunately using the rule above for bch (with sub > 1) generates faces which are not flat. // Therefore, we replace the vertices by the centers of their Voronoi cells // we do this only in the hyperbolic case -- it does not work correctly in the spherical case because of Voronoi not working as expected // the arguments for pt1 are the Voronoi centers for: (x,y,z) = (1,.5,0) // pt1 rearranges them to whatever (x,y,z) actually is array arg1 = {x, y, z}; auto pt1 = [&] (ld x1, ld y1, ld z1) { array arg0; for(int i=0; i<3; i++) { if(arg1[i] == 1) arg0[i] = x1; else if(arg1[i] == -1) arg0[i] = -x1; else if(arg1[i] == .5) arg0[i] = y1; else if(arg1[i] == -.5) arg0[i] = -y1; else if(arg1[i] == 0) arg0[i] = z1; else throw hr_exception("unknown number in pt1"); } return normalize(pt0(arg0)); }; hyperpoint a = pt1(0,0,0); hyperpoint b = pt1(2,0,0); hyperpoint c = pt1(1,1,1); hyperpoint d = pt1(1,1,-1); hyperpoint res = circumscribe(a, b, c, d); return res; }; auto add_face = [&] (std::initializer_list vh) { ss.faces.push_back(vh); }; const ld h = .5; if(mi == 0) { for(int s: {-1, 1}) { for(int i=0; i<3; i++) { if(bch) add_face({pt(0,.5,s), pt(.5,0,s), pt(0,-.5,s), pt(-.5,0,s)}); else add_face({pt(-1,-1,s), pt(-1,+1,s), pt(+1,+1,s), pt(+1,-1,s)}); tie(di[0], di[1], di[2]) = make_tuple(di[1], di[2], di[0]); } } if(bch) for(int u=0; u<8; u++) { for(int j=0; j<3; j++) if((u>>j)&1) di[j].flip(); add_face({pt(0,.5,1), pt(0,1,.5), pt(.5,1,0), pt(1,.5,0), pt(1,0,.5), pt(.5,0,1)}); for(int j=0; j<3; j++) if((u>>j)&1) di[j].flip(); } } else if(mi == 1) { auto& M = di[2].mirror; for(int s: {-1, 1}) { if(bch) add_face({pt(0,h,s), pt(h,0,s), pt(0,-h,s), pt(-h,0,s)}); else add_face({pt(-1,-1,s), pt(-1,+1,s), pt(+1,+1,s), pt(+1,-1,s)}); for(int i=0; i<2; i++) { if(bch) add_face({pt(1,0,-.5), pt(1,-.5,0), M*pt(1,0,-.5), pt(1,.5,0)}); else add_face({pt(-1,-1,-1), pt(-1,+1,-1), pt(-1,+1,+1), pt(-1,-1,+1)}); tie(di[0], di[1]) = make_tuple(di[1], di[0]); di[0].flip(); } } if(bch) for(ld s: {-1, 1}) for(int i=0; i<4; i++) { add_face({pt(0,.5,s), pt(0,1,s/2), pt(.5,1,0), pt(1,.5,0), pt(1,0,s/2), pt(.5,0,s)}); tie(di[0], di[1]) = make_tuple(di[1], di[0]); di[0].flip(); } } else { transmatrix spi = mi == 2 ? di[1].mirror * di[2].mirror : di[0].mirror * di[1].mirror; if(cgi.loop == 5) spi = spi * spi; vector spi_power = {Id}; for(int i=1; i f0, f1; for(auto P: spi_power) f0.push_back(bch ? P*pt(-1,-.5,0) : P*pt(-1,-1,-1)); for(auto P: spi_power) f1.push_back(bch ? P*pt(+1,-.5,0) : P*pt(+1,-1,-1)); ss.faces.push_back(f0); ss.faces.push_back(f1); if(bch) for(auto P: spi_power) for(int s: {-1,1}) add_face({P*pt(-.5*s,0,-1), P*pt(0,-.5,-1), P*pt(0,-1,-.5), P*pt(-.5*s,-1,0), P*pt(-1*s,-.5,0), P*pt(-1*s,0,-.5)}); } else { vector face_dirs = {Id}; for(int i=0; i f0; for(auto P: spi_power) f0.push_back(bch ? F*P*pt(-.5,0,-1) : F*P*pt(-1,-1,-1)); ss.faces.push_back(f0); } vector vertex_dirs; hyperpoint pter = normalize(pt0(make_array(-.5,-.5,-.5))); for(auto& F: face_dirs) for(auto& P: spi_power) { transmatrix T = F * P; bool fnd = false; for(auto T1: vertex_dirs) if(hdist(T * pter, T1*pter) < 1e-3) fnd = true; if(!fnd) vertex_dirs.push_back(T); } if(cgi.loop == 3) hassert(isize(vertex_dirs) == 4); if(cgi.loop == 5) hassert(isize(vertex_dirs) == 20); if(bch) for(auto& V: vertex_dirs) add_face({V*pt(-1,-.5,0), V*pt(-1,0,-.5), V*pt(-.5,0,-1), V*pt(0,-.5,-1), V*pt(0,-1,-.5), V*pt(-.5,-1,0)}); } } } } } EX void generate_bch_oct() { if(S7 != 6) throw hr_exception("generate_bch_oct but no cubes"); const int sub = subcube_count; if(1) { auto vx = abs(cgi.heptshape->faces[0][0][0]); auto vz = abs(cgi.heptshape->faces[0][0][3]); array co; // vx = 1; vz = 0; for(co[0]=-sub; co[0] sgn = {1,1,1}; if((co[1] ^ co[0]) & 1) co1[1]++, sgn[1] = -1; if((co[2] ^ co[0]) & 1) co1[2]++, sgn[2] = -1; hyperpoint ctr = Hypc; ctr[3] = vz * sub; auto pt = [&] (int m, ld x0, ld x1, ld x2) { hyperpoint res = ctr; auto x = make_array(x0, x1, x2); for(int i=0; i<3; i++) res[i] = vx * (co1[i] + x[(m+i)%3] * sgn[i]); return res; }; for(int it=0; it<2; it++) { cgi.subshapes.emplace_back(); auto &ss = cgi.subshapes.back(); for(int m=0; m<3; m++) { ss.faces.push_back({pt(m,0,0,0), pt(m,1,0,0), pt(m,1,0,.5), pt(m,.5,0,1), pt(m,0,0,1)}); ss.faces.push_back({pt(m,1,0,0), pt(m,1,0,.5), pt(m,1,.5,0) }); } ss.faces.push_back({pt(0,1,0,.5), pt(0,1,.5,0), pt(0,.5,1,0), pt(0,0,1,.5), pt(0,0,.5,1), pt(0,.5,0,1)}); for(int d=0; d<3; d++) co1[d] += sgn[d], sgn[d] *= -1; println(hlog, ss.faces); } } } } EX void generate_subcells() { cgi.subshapes.clear(); switch(variation) { case eVariation::subcubes: generate_plain_subcubes(); break; case eVariation::dual_subcubes: generate_special_subcubes(false); break; case eVariation::bch: generate_special_subcubes(true); break; case eVariation::bch_oct: generate_bch_oct(); break; case eVariation::coxeter: generate_coxeter(coxeter_param); break; case eVariation::pure: { cgi.subshapes.emplace_back(); cgi.subshapes[0].faces = cgi.heptshape->faces; break; } default: throw hr_exception("unknown variation in generate_subcells"); } for(auto& ss: cgi.subshapes) ss.compute_sub(); println(hlog, "subcells generated = ", isize(cgi.subshapes)); } void binary_rebase(heptagon *h, const transmatrix& V) { } void test(); #if HDR /** \brief vertex_adjacencies[heptagon id] is a list of other heptagons which are vertex adjacent * note: in case of ideal vertices this is just the face adjacency **/ struct vertex_adjacency_info { /** id of the adjacent heptagon */ int h_id; /** transition matrix to that heptagon */ transmatrix T; /** the sequence of moves we need to make to get there */ vector move_sequence; }; struct hrmap_closed3 : hrmap { vector allh; vector>> strafe_data; vector> tmatrices; vector> tmatrices_cell; vector acells; map > local_id; /* acells index, subshape index */ vector> acells_by_master; vector > vertex_adjacencies; vector>> move_sequences; transmatrix adj(heptagon *h, int d) override { return tmatrices[h->fieldval][d]; } transmatrix adj(cell *c, int d) override { return tmatrices_cell[local_id.at(c).first][d]; } heptagon *getOrigin() override { return allh[0]; } transmatrix relative_matrixc(cell *h2, cell *h1, const hyperpoint& hint) override; void initialize(int cell_count); vector& allcells() override { return acells; } ~hrmap_closed3() { clearfrom(getOrigin()); } subcellshape& get_cellshape(cell *c) override { if(PURE) return *cgi.heptshape ; int id = local_id.at(c).second; return cgi.subshapes[id]; } transmatrix master_relative(cell *c, bool get_inverse) override { int id = local_id.at(c).second; auto& ss = cgi.subshapes[id]; return get_inverse ? ss.from_cellcenter : ss.to_cellcenter; } void make_subconnections(); int wall_offset(cell *c) override; int shvid(cell *c) override { return local_id.at(c).second; } transmatrix ray_iadj(cell *c, int i) override; cellwalker strafe(cellwalker cw, int j) override { int id = local_id.at(cw.at).first; return cellwalker(cw.at->cmove(j), strafe_data[id][j][cw.spin]); } const vector& get_move_seq(cell *c, int i) override { int id = local_id.at(c).first; return move_sequences[id][i]; } }; struct hrmap_quotient3 : hrmap_closed3 { }; #endif transmatrix hrmap_closed3::ray_iadj(cell *c, int i) { if(PURE) return iadj(c, i); auto& v = get_face_vertices(c, i); hyperpoint h = project_on_triangle(v[0], v[1], v[2]); transmatrix T = rspintox(h); return T * xpush(-2*hdist0(h)) * spintox(h); } int hrmap_closed3::wall_offset(cell *c) { if(PURE) return 0; auto& wo = cgi.walloffsets[ local_id.at(c).second ]; if(wo.second == nullptr) wo.second = c; return wo.first; } EX const vector& get_face_vertices(cell *c, int d) { return cgi.subshapes[currentmap->shvid(c)].faces_local[d]; } EX int get_face_vertex_count(cell *c, int d) { return isize(get_face_vertices(c, d)); } void hrmap_closed3::initialize(int big_cell_count) { allh.resize(big_cell_count); tmatrices.resize(big_cell_count); acells.clear(); for(int a=0; afieldval = a; } } const static bool testing_subconnections = false; void hrmap_closed3::make_subconnections() { auto& ss = cgi.subshapes; auto& vas = vertex_adjacencies; vas.resize(isize(allh)); for(int a=0; a buckets; for(auto& v: cgi.heptshape->vertices_only) buckets.insert(bucketer(v)); if(cgflags & qIDEAL) { for(int d=0; dmove(d)->fieldval, T, {d}}); } } else for(int i=0; ivertices_only) if(buckets.count(bucketer(T*w))) found_va = true; if(!found_va) continue; va.emplace_back(vertex_adjacency_info{allh[va[i].h_id]->move(d)->fieldval, T, va[i].move_sequence}); va.back().move_sequence.push_back(d); } } } map by_sides; vector > > which_cell_0; which_cell_0.resize(isize(allh)); acells_by_master.resize(isize(allh)); for(int a=0; ac7) allh[a]->c7 = c; local_id[c] = {isize(acells), id}; acells.push_back(c); acells_by_master[a].push_back(c); which_cell_0[a][bucketer(cc)].push_back(cc); } } println(hlog, "found ", isize(acells), " cells, ", by_sides); tmatrices_cell.resize(isize(acells)); move_sequences.resize(isize(acells)); int failures = 0; vector > > > which_cell; which_cell.resize(isize(allh)); for(cell *c: acells) { int id = local_id[c].second; for(int i=0; itype; i++) which_cell[c->master->fieldval][bucketer(ss[id].face_centers[i])].emplace_back(c, i); } strafe_data.resize(isize(acells)); for(cell *c: acells) { int id = local_id[c].second; int cid = local_id[c].first; auto& tmcell = tmatrices_cell[cid]; vector foundtab; vector> foundtab_ids; strafe_data[cid].resize(c->type); for(int i=0; itype; i++) { int found = 0; hyperpoint ctr = ss[id].face_centers[i]; transmatrix T1 = Id; int h_id = c->master->fieldval; vector path; while(true) { int d = -1; ld dist = hdist0(ctr); for(int d1=0; d1move(d)->fieldval; } for(auto& va: vertex_adjacencies[h_id]) { hyperpoint ctr1 = iso_inverse(va.T) * ctr; auto bucket = bucketer(ctr1); for(auto p: which_cell[va.h_id][bucket]) { cell *c1 = p.first; int j = p.second; int id1 = local_id[c1].second; if(hdist(ctr1, ss[id1].face_centers[j]) < 1e-6) { transmatrix T2 = T1 * va.T; if(id == id1 && eqmatrix(T2, Id)) continue; c->c.connect(i, c1, j, false); if(!found) { tmcell.push_back(ss[id].from_cellcenter * T2 * ss[id1].to_cellcenter); if(elliptic) fixelliptic(tmcell.back()); auto& ms = move_sequences[local_id[c].first]; ms.push_back(path); for(auto dir: va.move_sequence) ms.back().push_back(dir); auto& sd = strafe_data[cid][i]; sd.resize(c->type, -1); for(int i1=0; i1type; i1++) { set facevertices; for(auto v: ss[id].faces[i1]) facevertices.insert(bucketer(v)); if(ss[id].dirdist[i][i1] == 1) { int found_strafe = 0; for(int j1=0; j1type; j1++) if(j1 != j) { int num = 0; for(auto v: ss[id1].faces[j1]) if(facevertices.count(bucketer(T2*v))) num++; if(num == 2) sd[i1] = j1, found_strafe++; } if(found_strafe != 1) println(hlog, "found_strafe = ", found_strafe); } } /* for bch, also provide second-order strafe */ if(variation == eVariation::bch) for(int i1=0; i1type; i1++) { if(ss[id].dirdist[i][i1] != 2) continue; if(isize(ss[id].faces[i]) == 6) continue; if(isize(ss[id].faces[i1]) == 6) continue; vector fac; for(int i2=0; i2type; i2++) if(ss[id].dirdist[i][i2] == 1 && ss[id].dirdist[i2][i1] == 1) fac.push_back(sd[i2]); if(isize(fac) != 2) { println(hlog, "fac= ", fac); throw hr_exception("fac error"); } int found_strafe = 0; for(int j1=0; j1type; j1++) if(j1 != j) if(ss[id1].dirdist[j1][fac[0]] == 1) if(ss[id1].dirdist[j1][fac[1]] == 1) { sd[i1] = j1; found_strafe++; } if(found_strafe != 1) println(hlog, "found_strafe = ", found_strafe, " (second order)"); } } foundtab_ids.emplace_back(va.h_id, id1, j); found++; } } if(found && !testing_subconnections) break; } if(testing_subconnections && !found) { c->c.connect(i, c, i, false); tmcell.push_back(Id); } foundtab.push_back(found); if(found != 1) failures++; } if(testing_subconnections) println(hlog, "foundtab = ", foundtab); } println(hlog, "total failures = ", failures); if(failures && !testing_subconnections) throw hr_exception("hrmap_closed3 failures"); } transmatrix hrmap_closed3::relative_matrixc(cell *c2, cell *c1, const hyperpoint& hint) { if(c1 == c2) return Id; int d = hr::celldistance(c2, c1); for(int a=0; atype; a++) if(hr::celldistance(c2, c1->move(a)) < d) return adj(c1, a) * relative_matrix(c2, c1->move(a), hint); for(int a=0; atype; a++) println(hlog, "d=", d, " vs ", hr::celldistance(c2, c1->move(a))); println(hlog, "error in hrmap_quotient3:::relative_matrix"); return Id; } #if CAP_CRYSTAL int encode_coord(int bits, const crystal::coord& co) { int c = 0; for(int i=0; i<(1<>1) & ((1<>= bits; return co; } struct hrmap_from_crystal : hrmap_quotient3 { int bits; hrmap_from_crystal(int b) : bits(b) { int size = 1 << (4*bits); initialize(size); if(1) { auto m = crystal::new_map(); dynamicval cm(currentmap, m); for(int a=0; acmove(d))); allh[a]->c.connect(d, allh[b], h1->c.spin(d), false); tmatrices[a].push_back(crystal::get_adj(h1, d)); } } delete m; } make_subconnections(); } }; #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 moveid(S7), movedir(lgr); for(int s=0; sfullv[s] * cgi.adjmoves[0], cgi.adjmoves[i])) moveid[i] = s; for(int s=0; sfullv[s]) * cgi.adjmoves[0]), tC0(cgi.adjmoves[i])) < 1e-4) movedir[s] = i; for(int a=0; amatcode[ f->mmul(f->mmul(f->matrices[k], f->matrices[moveid[b]]), f->P) ]; for(int l=0; lgmul(k, l) % lgr == 0) { tmatrices[a][b] = cgi.adjmoves[b] * f->fullv[l]; allh[a]->c.connect(b, allh[k/lgr], movedir[l], false); } } } make_subconnections(); create_patterns(); } set plane; void make_plane(cellwalker cw) { if(plane.count(cw)) return; plane.insert(cw); auto& ss = get_cellshape(cw.at); for(int i=0; itype; i++) if(ss.dirdist[i][cw.spin] == 1) make_plane(strafe(cw, i)); } void create_patterns() { DEBB(DF_GEOM, ("creating pattern = ", isize(allh))); if(!PURE) { println(hlog, "create_patterns not implemented"); return; } // also, strafe needs currentmap dynamicval 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 plane_indices; for(auto cw: plane) plane_indices.insert(cw.at->master->fieldval); int fN = isize(f->matrices); set nwi; for(int i=0; igmul(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 periods; int flip(int x) { return (x+6) % 12; } void build_reps() { // start_game(); auto& hsh = get_hsh(); set boundaries; for(int a=0; a<12; a++) for(int b=0; b<12; b++) if(hsh.dirdist[a][b] == 1) { 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 = hsh.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 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 = std::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; bc.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); } make_subconnections(); } }; 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); } } make_subconnections(); } }; } EX bool minimize_quotient_maps = false; EX bool subrule = false; EX bool strafe_test = false; hrmap_quotient3 *gen_quotient_map(bool minimized, fieldpattern::fpattern &fp) { #if CAP_FIELD #if CAP_CRYSTAL if(geometry == gSpace344) { return new hrmap_from_crystal(minimized ? 1 : 2); } else #endif if(geometry == gSpace535 && minimized) { return new seifert_weber::hrmap_singlecell(108*degree); } else if(geometry == gSpace535) return new seifert_weber::hrmap_seifert_cover; else if(hyperbolic) { return new hrmap_field3(&fp); } else if(geometry == gCubeTiling) return new seifert_weber::hrmap_singlecell(0); #endif return nullptr; } struct hrmap_h3_abstract : hrmap { reg3::hrmap_quotient3 *quotient_map; transmatrix adj(heptagon *h, int d) override { throw hr_exception("any adj"); } transmatrix relative_matrixc(cell *c2, cell *c1, const hyperpoint& hint) override { if(PURE) return relative_matrix(c2->master, c1->master, hint); return relative_matrix_via_masters(c2, c1, hint); } transmatrix master_relative(cell *c, bool get_inverse) override { if(PURE) return Id; int aid = cell_id.at(c); return quotient_map->master_relative(quotient_map->acells[aid], get_inverse); } int shvid(cell *c) override { if(PURE) return 0; if(!cell_id.count(c)) return quotient_map->shvid(c); int aid = cell_id.at(c); return quotient_map->shvid(quotient_map->acells[aid]); } int wall_offset(cell *c) override { if(PURE) return 0; if(!cell_id.count(c)) return quotient_map->wall_offset(c); /* necessary because ray samples are from quotient_map */ int aid = cell_id.at(c); return quotient_map->wall_offset(quotient_map->acells[aid]); } transmatrix adj(cell *c, int d) override { if(PURE) return adj(c->master, d); if(!cell_id.count(c)) return quotient_map->adj(c, d); /* necessary because ray samples are from quotient_map */ int aid = cell_id.at(c); return quotient_map->tmatrices_cell[aid][d]; } subcellshape& get_cellshape(cell *c) override { if(PURE) return *cgi.heptshape; int aid = cell_id.at(c); return quotient_map->get_cellshape(quotient_map->acells[aid]); } map cell_id; map, cell*> cell_at; cell *get_cell_at(heptagon *h, int acell_id) { pair p(h, acell_id); auto& ca = cell_at[p]; if(!ca) { ca = newCell(quotient_map->acells[acell_id]->type, h); cell_id[ca] = acell_id; if(!h->c7) h->c7 = ca; } return ca; } void find_cell_connection(cell *c, int d) override { if(PURE) { auto h = c->master->cmove(d); c->c.connect(d, h->c7, c->master->c.spin(d), false); return; } int id = cell_id.at(c); heptagon *h = c->master; for(int dir: quotient_map->move_sequences[id][d]) h = h->cmove(dir); auto ac = quotient_map->acells[id]; cell *c1 = get_cell_at(h, quotient_map->local_id[ac->move(d)].first); c->c.connect(d, c1, ac->c.spin(d), false); } transmatrix ray_iadj(cell *c, int i) override { if(PURE) return iadj(c, i); if(!cell_id.count(c)) return quotient_map->ray_iadj(c, i); /* necessary because ray samples are from quotient_map */ int aid = cell_id.at(c); return quotient_map->ray_iadj(quotient_map->acells[aid], i); } const vector& get_move_seq(cell *c, int i) override { int aid = cell_id.at(c); return quotient_map->get_move_seq(quotient_map->acells[aid], i); } cellwalker strafe(cellwalker cw, int j) override { int aid = PURE ? cw.at->master->fieldval : cell_id.at(cw.at); auto ress = quotient_map->strafe(cellwalker(quotient_map->acells[aid], cw.spin), j); cellwalker res = cellwalker(cw.at->cmove(j), ress.spin); if(PURE && strafe_test) { hyperpoint hfront = tC0(cgi.adjmoves[cw.spin]); cw.at->cmove(j); transmatrix T = currentmap->adj(cw.at, j); for(int i=0; ic.spin(j)) if(hdist(hfront, T * tC0(cgi.adjmoves[i])) < cgi.strafedist + .01) { auto resx = cellwalker(cw.at->cmove(j), i); if(res == resx) return res; if(PURE && res != resx) println(hlog, "h3: ", res, " vs ", resx); } throw hr_exception("incorrect strafe"); } return res; } }; struct hrmap_h3 : hrmap_h3_abstract { heptagon *origin; hrmap *binary_map; vector extra_origins; map> reg_gmatrix; map > > altmap; vector& allcells() override { return hrmap::allcells(); } hrmap_h3() { origin = init_heptagon(S7); heptagon& h = *origin; h.s = hsOrigin; if(PURE) h.c7 = newCell(S7, origin); worst_error1 = 0, worst_error2 = 0; dynamicval cr(currentmap, this); heptagon *alt = NULL; transmatrix T = Id; binary_map = nullptr; quotient_map = gen_quotient_map(minimize_quotient_maps, currfp); h.zebraval = quotient_map ? quotient_map->allh[0]->zebraval : 0; #if CAP_BT if(hyperbolic) { dynamicval 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); if(PURE) { 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); } } if(!PURE) get_cell_at(origin, 0); } ld worst_error1, worst_error2; heptagon *getOrigin() override { return origin; } void fix_distances(heptagon *h, heptagon *h2) { vector 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; imove(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 g(geometry, gBinary3); dynamicval cm(currentmap, binary_map); for(int i=0; itype; 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 #if CAP_BT if(hyperbolic) { dynamicval g(geometry, gBinary3); dynamicval cm(currentmap, binary_map); binary_map->virtualRebase(alt, T); } #endif fixmatrix(T); auto hT = tC0(T); 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 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 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; d2c.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 heptagon *created = init_heptagon(S7); if(PURE && parent->c7) created->c7 = newCell(S7, created); #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_h3() { #if CAP_BT if(binary_map) { dynamicval g(geometry, gBinary3); delete binary_map; } #endif if(quotient_map) delete quotient_map; clearfrom(origin); for(auto e: extra_origins) clearfrom(e); } map reducers; bool link_alt(heptagon *h, heptagon *alt, hstate firststate, int dir) override { altdist(h) = 0; if(firststate != hsOrigin) reducers[h] = dir; return true; } void extend_altmap(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; icmove(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_matrixh(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 g(geometry, gBinary3); dynamicval 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; } cell* gen_extra_origin(int fv) override { auto orig = isize(extra_origins) ? extra_origins.back() : origin; auto& p1 = reg_gmatrix[orig]; heptagon *alt = p1.first; transmatrix T = p1.second; for(int a=0; a<10; a++) { T = T * xpush(10); #if CAP_BT if(hyperbolic) { dynamicval g(geometry, gBinary3); dynamicval cm(currentmap, binary_map); binary_map->virtualRebase(alt, T); fixmatrix(T); } #endif } heptagon *created = init_heptagon(S7); created->s = hsOrigin; created->fieldval = quotient_map->acells[fv]->master->fieldval; reg_gmatrix[created] = make_pair(alt, T); altmap[alt].emplace_back(created, T); extra_origins.push_back(created); return get_cell_at(created, fv); } }; EX int get_aid(cell *c) { auto m = dynamic_cast (currentmap); if(!m) throw hr_exception("get_aid incorrect"); if(PURE) return c->master->fieldval; return m->cell_id[c]; } EX int get_size_of_aid(int aid) { auto m = dynamic_cast (currentmap); if(!m) throw hr_exception("get_size_of_fv incorrect"); return m->quotient_map->acells[aid]->type; }; struct hrmap_sphere3 : hrmap_closed3 { vector locations; hrmap_sphere3() { heptagon *h = init_heptagon(S7); h->s = hsOrigin; locations.push_back(Id); allh.push_back(h); for(int i=0; ic.connect(d, allh[i1], d2, false); fb++; tmi.push_back(inverse(T1) * locations[i1]); } } if(fb != 1) throw hr_exception("friend not found"); goto next_d; } if(1) { int d2 = 0; if(hopf) { for(d2=0; d2zebraval = hrand(10); h->fieldval = isize(allh); h->fiftyval = 9999; allh.push_back(h); locations.push_back(T); if(isnan(T[0][0])) exit(1); allh[i]->c.connect(d, h, d2, false); tmi.push_back(inverse(T1) * T); if(elliptic) fixelliptic(tmi.back()); } next_d: ; } } make_subconnections(); } transmatrix relative_matrixh(heptagon *h2, heptagon *h1, const hyperpoint& hint) override { return iso_inverse(locations[h1->fieldval]) * locations[h2->fieldval]; } }; EX const transmatrix& get_sphere_loc(int v) { return ((hrmap_sphere3*)currentmap)->locations[v]; } struct ruleset { fieldpattern::fpattern fp; vector root; string other; vector children; vector childpos; vector otherpos; virtual hrmap_quotient3 *qmap() = 0; ruleset() : fp(0) {} void load_ruleset(string fname) { shstream ins(decompress_string(read_file_as_string(fname))); dynamicval q(fieldpattern::use_quotient_fp, true); hread_fpattern(ins, fp); hread(ins, root); hread(ins, children); hread(ins, other); // hread(ins, childpos); int t = S7; int qty = isize(children) / t; for(int i=0; i<=qty; i++) childpos.push_back(i * t); } void load_ruleset_new(string fname) { shstream ins(decompress_string(read_file_as_string(fname))); ins.read(ins.vernum); mapstream::load_geometry(ins); ins.read(fieldpattern::use_rule_fp); ins.read(fieldpattern::use_quotient_fp); ins.read(reg3::minimize_quotient_maps); hread_fpattern(ins, fp); hread(ins, root); hread(ins, children); hread(ins, other); hread(ins, childpos); } /** \brief address = (fieldvalue, state) */ typedef pair address; /** nles[x] lists the addresses from which we can reach address x * without ever ending in the starting point */ map> nonlooping_earlier_states; vector> possible_states; virtual int connection(int fv, int d) = 0; void find_mappings() { 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++; } } /* find all back paths */ auto &nles = nonlooping_earlier_states; nles.clear(); vector
bfs; int qty = isize(root); for(int i=0; ifieldval = h->fieldval; if(firststate == hsOrigin) { alt->fiftyval = root[alt->fieldval % isize(root)]; return true; } vector& choices = possible_states[alt->fieldval]; vector choices2; for(auto c: choices) { bool ok = true; for(int d=0; dcmove(d)->distance < h->distance) if(children[childpos[c]+d] == -1) ok = false; if(ok) choices2.push_back(c); } alt->fiftyval = hrand_elt(choices2, -1); return alt->fiftyval != -1; } }; struct hrmap_h3_rule : hrmap_h3_abstract, ruleset { heptagon *origin; reg3::hrmap_quotient3 *emerald_map; hrmap_quotient3 *qmap() override { return quotient_map; } hrmap_h3_rule() { load_ruleset(get_rule_filename()); quotient_map = gen_quotient_map(is_minimized(), fp); find_mappings(); origin = init_heptagon(S7); heptagon& h = *origin; h.s = hsOrigin; h.fiftyval = root[0]; if(PURE) h.c7 = newCell(S7, origin); emerald_map = gen_quotient_map(false, currfp); h.emeraldval = 0; if(!PURE) get_cell_at(origin, 0); } int connection(int fv, int d) override { return qmap()->allh[fv]->move(d)->fieldval; } heptagon *getOrigin() override { return origin; } #define DEB 0 heptagon *counterpart(heptagon *h) { return quotient_map->allh[h->fieldval]; } vector evmemo; void find_emeraldval(heptagon *target, heptagon *parent, int d) { generate_cellrotations(); auto& cr = cgi.cellrotations; if(evmemo.empty()) { println(hlog, "starting"); map 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; iemeraldval; 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(PURE && 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 possible; int pfv = parent->fieldval; 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_h3_rule() { if(quotient_map) delete quotient_map; clearfrom(origin); } transmatrix adj(heptagon *h, int d) override { return quotient_map->adj(h, d); } transmatrix relative_matrixh(heptagon *h2, heptagon *h1, const hyperpoint& hint) override { return relative_matrix_recursive(h2, h1); } virtual bool link_alt(heptagon *h, heptagon *alt, hstate firststate, int dir) override { return ruleset_link_alt(h, alt, firststate, dir); } }; vector starts = {nullptr}; struct hrmap_h3_subrule : hrmap, ruleset { heptagon *origin; hrmap_quotient3 *quotient_map; hrmap_quotient3 *qmap() override { return quotient_map; } int connection(int fv, int d) override { return quotient_map->local_id[quotient_map->acells[fv]->move(d)].first; } void explain_conflict(vector hs) { vector v; for(auto h: hs) v.push_back(h.at); int N = isize(v); vector> paths(N); while(true) { bool eq = true; for(auto c: v) if(c != v[0]) eq = false; if(eq) break; int mindist = 999999; int maxdist = -999999; for(auto c: v) { if(c->distance < mindist) mindist = c->distance; if(c->distance > maxdist) maxdist = c->distance; } int goal = min(mindist, maxdist-1); println(hlog, "mindist = ", mindist, " maxdist = ", maxdist, " goal = ", goal); int id = 0; for(auto& c: v) { while(c->distance > goal) { println(hlog, c, " distance is ", c); int d = find_parent(c); paths[id].push_back(c->c.spin(d)); c = c->move(d); } id++; } } for(auto& p: paths) reverse(p.begin(), p.end()); hs.push_back(heptspin(v[0], find_parent(v[0]))); for(auto h: hs) { println(hlog, h, " : dist = ", h.at->distance, " id = ", h.at->fiftyval, " qid = ", h.at->fieldval); } println(hlog, "paths = ", paths); } hrmap_h3_subrule() { println(hlog, "loading a subrule ruleset"); load_ruleset_new(get_rule_filename()); quotient_map = gen_quotient_map(is_minimized(), fp); int t = quotient_map->acells[0]->type; find_mappings(); origin = init_heptagon(t); heptagon& h = *origin; h.s = hsOrigin; h.fieldval = 0; h.fiftyval = root[0]; h.c7 = newCell(t, origin); } heptagon *getOrigin() override { return origin; } heptagon *create_step(heptagon *parent, int d) override { heptspin parentd(parent, d); if(starts[isize(starts)/2] == parentd) { int i = 0; vector cut; for(auto s: starts) if(i++ >= isize(starts)/2) cut.push_back(s); println(hlog, "cycle detected is ", cut); explain_conflict(cut); throw hr_exception("create_step cycle detected"); } starts.push_back(parentd); finalizer f([] { starts.pop_back(); }); int id = parent->fiftyval; if(id < 0) id += (1<<16); int qid = parent->fieldval; int d2 = quotient_map->acells[qid]->c.spin(d); int qid2 = quotient_map->local_id[quotient_map->acells[qid]->move(d)].first; heptagon *res = nullptr; int id1 = children[childpos[id]+d]; int pos = otherpos[childpos[id]+d]; if(id1 < -1) id1 += (1<<16); if(id1 != -1) { int t = childpos[id1+1] - childpos[id1]; res = init_heptagon(t); res->c7 = newCell(t, res); res->fieldval = qid2; res->distance = parent->distance + 1; res->fiftyval = id1; } else if(other[pos] == ('A' + d) && other[pos+1] == ',') { vector possible; for(auto s: nonlooping_earlier_states[address{qid, id}]) possible.push_back(s.second); if(possible.empty()) throw hr_exception("impossible"); id1 = hrand_elt(possible, 0); int t = childpos[id1+1] - childpos[id1]; res = init_heptagon(t); res->alt = parent->alt; res->fieldval = qid2; res->distance = parent->distance - 1; res->fiftyval = id1; } else { heptagon *at = parent; while(other[pos] != ',') { int dir = other[pos++] - 'a'; at = at->cmove(dir); } res = at; } if(!res) throw hr_exception("res missing"); if(res->move(d2)) { heptspin a(res, d2); heptspin b = a + wstep; heptspin c(parent, d); println(hlog, "res conflict: ", a, " already connected to ", b, " and should be connected to ", c); explain_conflict({a, b, c}); } res->c.connect(d2, parent, d, false); return res; } int find_parent(heptagon *h) { int id = h->fiftyval; int pos = childpos[id]; for(int i=0; iadj(quotient_map->acells[h->fieldval], d); } transmatrix relative_matrixh(heptagon *h2, heptagon *h1, const hyperpoint& hint) override { return relative_matrix_recursive(h2, h1); } int shvid(cell *c) override { return quotient_map->shvid(quotient_map->acells[c->master->fieldval]); } int wall_offset(cell *c) override { return quotient_map->wall_offset(quotient_map->acells[c->master->fieldval]); } subcellshape& get_cellshape(cell *c) override { return quotient_map->get_cellshape(quotient_map->acells[c->master->fieldval]); } transmatrix ray_iadj(cell *c, int i) override { return quotient_map->ray_iadj(quotient_map->acells[c->master->fieldval], i); } cellwalker strafe(cellwalker cw, int j) override { int aid = cw.at->master->fieldval; auto ress = quotient_map->strafe(cellwalker(quotient_map->acells[aid], cw.spin), j); return cellwalker(cw.at->cmove(j), ress.spin); } bool link_alt(heptagon *h, heptagon *alt, hstate firststate, int dir) override { return ruleset_link_alt(h, alt, firststate, dir); } }; struct hrmap_h3_rule_alt : hrmap { heptagon *origin; hrmap_h3_rule_alt(heptagon *o) { origin = o; } }; EX hrmap *new_alt_map(heptagon *o) { println(hlog, "new_alt_map called"); return new hrmap_h3_rule_alt(o); } EX bool reg3_rule_available = true; EX string other_rule = ""; EX bool is_minimized() { return geometry == gSpace535; } 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() != ""; } ruleset& get_ruleset() { if(subrule) return *((hrmap_h3_subrule*)currentmap); if(in_rule()) return *((hrmap_h3_rule*)currentmap); throw hr_exception("get_ruleset called but not in rule"); } EX int rule_get_root(int i) { return get_ruleset().root[i]; } EX const vector& rule_get_children() { return get_ruleset().children; } EX const vector& rule_get_childpos() { return get_ruleset().childpos; } 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(subrule) return new hrmap_h3_subrule; if(in_rule()) return new hrmap_h3_rule; if(sphere) return new hrmap_sphere3; return new hrmap_h3; } hrmap_h3* hypmap() { return ((hrmap_h3*) currentmap); } EX int quotient_count() { return isize(hypmap()->quotient_map->allh); } EX int quotient_count_sub() { return isize(hypmap()->quotient_map->acells); } /** 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 s1 = {c1}; vector s2 = {c2}; int best = 99999999; int d0 = 0; auto go_nearer = [&] (vector& v, int& d) { vector 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 && PURE) return celldistance_534(c1, c2); auto r = hypmap(); 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 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) { if(sphere) { auto m = currentmap; hyperpoint h = tC0(m->relative_matrix(c->master, m->getOrigin(), 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; } auto m = hypmap(); if(cgflags & qSINGLE) return true; if(fake::in()) return FPIU(reg3::pseudohept(c)); // 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 (currentmap); if(mr) { if(geometry == gSpace535) return c->master->fieldval % 31 == 0; return c->master->fieldval == 0; } auto ms = dynamic_cast (currentmap); if(ms) 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 0.01) continue; vector perm(S7); for(int x=0; xreg_gmatrix[c->master].second); h = bt::deparabolic3(h); return regmap()->reg_gmatrix[c->master].first->distance * log(2) - h[0]; } map, 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 v[2]; v[0].push_back(c1); v[1].push_back(c2); int steps = 0; map 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 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<> 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 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 all; vector formulas; formulas.push_back(""); all.push_back(Id); auto& faces = cgi.heptshape->faces; hyperpoint v = faces[0][0]; auto add = [&] (transmatrix T) { for(int i=0; i 5) return; for(auto& hv: faces) 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