1 // Hyperbolic Rogue -- regular honeycombs
2 // Copyright (C) 2011-2019 Zeno Rogue, see 'hyper.cpp' for details
3
4 /** \file reg3.cpp
5 * \brief regular honeycombs
6 *
7 * works with spherical and hyperbolic ones -- Euclidean cubic tiling implemented in euclid.cpp
8 * includes non-quotient spaces as well as field quotient and elliptic spaces
9 * hyperbolic honeycombs rely on bt:: to deal with floating point errors (just like archimedean)
10 */
11
12 #include "hyper.h"
13 namespace hr {
14 #if MAXMDIM >= 4
15
final_coords(hyperpoint h)16 EX hyperpoint final_coords(hyperpoint h) {
17 if(sn::in() || !bt::in())
18 return ultra_normalize(h);
19 #if CAP_BT
20 if(bt::in()) {
21 ld yy = log(2) / 2;
22 return bt::parabolic3(h[0], h[1]) * xpush0(yy*h[2]);
23 }
24 #endif
25 return h;
26 }
27
compute_common()28 void subcellshape::compute_common() {
29 reg3::make_vertices_only(vertices_only, faces);
30
31 faces_local = faces;
32 for(auto& face: faces_local) for(auto& v: face) v = from_cellcenter * final_coords(v);
33
34 vertices_only_local = vertices_only;
35 for(auto& v: vertices_only_local) v = from_cellcenter * final_coords(v);
36
37 int N = isize(faces);
38
39 dirdist.resize(N);
40 for(int i=0; i<N; i++) {
41 auto& da = dirdist[i];
42 da.resize(N, false);
43 set<unsigned> cface;
44 for(auto& v: faces[i]) cface.insert(bucketer(v));
45 for(int j=0; j<N; j++) {
46 int mutual = 0;
47 for(auto& w: faces[j]) if(cface.count(bucketer(w))) mutual++;
48 da[j] = i == j ? 0 : mutual == 2 ? 1 : INFD;
49 }
50 }
51 floyd_warshall(dirdist);
52
53 next_dir.resize(N);
54 for(int a=0; a<N; a++) next_dir[a].resize(N);
55
56 for(int a=0; a<N; a++)
57 for(int b=0; b<N; b++)
58 if(dirdist[a][b] == 1)
59 for(int c=0; c<N; c++)
60 if(dirdist[a][c] == 1 && dirdist[b][c] == 1) {
61 transmatrix t = build_matrix(tC0(cgi.adjmoves[a]), tC0(cgi.adjmoves[b]), tC0(cgi.adjmoves[c]), C0);
62 if(det(t) > 0) next_dir[a][b] = c;
63 }
64 }
65
compute_hept()66 void subcellshape::compute_hept() {
67 cellcenter = C0;
68 to_cellcenter = Id;
69 from_cellcenter = Id;
70 compute_common();
71 }
72
compute_sub()73 void subcellshape::compute_sub() {
74 hyperpoint gres = Hypc;
75 for(auto& face: faces) {
76 hyperpoint res = Hypc;
77 for(auto& vertex: face)
78 res += vertex;
79 face_centers.push_back(normalize(res));
80 gres += res;
81 }
82 cellcenter = normalize(gres);
83 to_cellcenter = rgpushxto0(cellcenter);
84 from_cellcenter = gpushxto0(cellcenter);
85 compute_common();
86 }
87
88 /** \brief regular three-dimensional tessellations */
89 EX namespace reg3 {
90
91 EX int subcube_count = 1;
92
93 EX flagtype coxeter_param = 0;
94
95 const flagtype cox_othercell = 1;
96 const flagtype cox_midedges = 2;
97 const flagtype cox_vertices = 4;
98
99 #if HDR
altdist(heptagon * h)100 inline short& altdist(heptagon *h) { return h->emeraldval; }
101 #endif
102
103 EX int extra_verification;
104
105 EX bool ultra_mirror_on;
106
ultra_mirror_in()107 EX bool ultra_mirror_in() { return (cgflags & qULTRA) && ultra_mirror_on; }
108
in()109 EX bool in() {
110 if(fake::in()) return FPIU(in());
111 return WDIM == 3 && !euclid && !bt::in() && !nonisotropic && !hybri && !kite::in();
112 }
113
compute_ultra()114 EX void compute_ultra() {
115 cgi.ultra_mirror_part = .99;
116 cgi.ultra_material_part = .99;
117
118 cgi.ultra_mirrors.clear();
119
120 if(cgflags & qULTRA) {
121
122 for(auto& v: cgi.heptshape->vertices_only) {
123
124 hyperpoint nei;
125 auto& faces = cgi.heptshape->faces;
126
127 for(int i=0; i<isize(faces); i++)
128 for(int j=0; j<isize(faces[i]); j++)
129 if(sqhypot_d(WDIM, faces[i][j]-v) < 1e-6)
130 nei = faces[i][j?j-1:j+1];
131
132 transmatrix T = spintox(v);
133 hyperpoint a = T * v;
134 hyperpoint b = T * nei;
135 ld f0 = 0.5;
136 ld f1 = binsearch(0.5, 1, [&] (ld d) {
137 hyperpoint c = lerp(b, a, d);
138 if(debugflags & DF_GEOM)
139 println(hlog, "d=", d, " c= ", c, " material = ", material(c));
140 return material(c) <= 0;
141 });
142 cgi.ultra_material_part = f1;
143 auto f = [&] (ld d) {
144 hyperpoint c = lerp(b, a, d);
145 c = normalize(c);
146 return c[1] * c[1] + c[2] * c[2];
147 };
148 for(int it=0; it<100; it++) {
149 ld fa = (f0*2+f1) / 3;
150 ld fb = (f0*1+f1*2) / 3;
151 if(debugflags & DF_GEOM)
152 println(hlog, "f(", fa, ") = ", f(fa), " f(", fb, ") = ", f(fb));
153 if(f(fa) > f(fb)) f0 = fa;
154 else f1 = fb;
155 }
156
157 cgi.ultra_mirror_part = f0;
158
159 hyperpoint c = lerp(b, a, f0);
160 c = normalize(c);
161 c[1] = c[2] = 0;
162 c = normalize(c);
163 cgi.ultra_mirror_dist = hdist0(c);
164
165 if(cgi.ultra_mirror_part >= 1-1e-6) continue;
166
167 cgi.ultra_mirrors.push_back(rspintox(v) * xpush(cgi.ultra_mirror_dist*2) * MirrorX * spintox(v));
168 }
169 }
170 }
171
make_vertices_only(vector<hyperpoint> & vo,const vector<vector<hyperpoint>> & csh)172 EX void make_vertices_only(vector<hyperpoint>& vo, const vector<vector<hyperpoint>>& csh) {
173 vo.clear();
174 for(auto& v: csh)
175 for(hyperpoint h: v) {
176 bool found = false;
177 for(hyperpoint h2: vo) if(hdist(h, h2) < 1e-6) found = true;
178 if(!found) vo.push_back(h);
179 }
180 }
181
generate()182 EX void generate() {
183
184 if(fake::in()) {
185 fake::generate();
186 return;
187 }
188
189 auto& hsh = get_hsh();
190
191 int& loop = cgi.loop;
192 int& face = cgi.face;
193 auto& spins = cgi.spins;
194 auto& cellshape = hsh.faces;
195 auto& adjcheck = cgi.adjcheck;
196
197 int& mid = cgi.schmid;
198 mid = 3;
199
200 face = 3;
201 if(S7 == 6) face = 4;
202 if(S7 == 8) mid = 4;
203 if(S7 == 12) face = 5;
204 if(S7 == 20) mid = 5;
205 /* icosahedron not implemented */
206 loop = ginf[geometry].tiling_name[5] - '0';
207 DEBB(DF_GEOM, ("face = ", face, " loop = ", loop, " S7 = ", S7));
208
209 ld angle_between_faces, hcrossf;
210
211 /* frontal face direction */
212 hyperpoint h0, h1, h2, h3, h012, h013;
213
214 if(1) {
215 dynamicval<eGeometry> dg(geometry, gSphere);
216 angle_between_faces = edge_of_triangle_with_angles(2*M_PI/mid, M_PI/face, M_PI/face);
217
218 h0 = xtangent(1);
219 h1 = cspin(0, 1, angle_between_faces) * h0;
220 h2 = cspin(1, 2, 2*M_PI/face) * h1;
221 h3 = cspin(1, 2, -2*M_PI/face) * h1;
222
223 hcrossf = edge_of_triangle_with_angles(M_PI/2, M_PI/mid, M_PI/face);
224
225 h012 = cspin(1, 2, M_PI/face) * cspin(0, 1, hcrossf) * h0;
226 h013 = cspin(1, 2, -M_PI/face) * cspin(0, 1, hcrossf) * h0;
227 }
228
229 for(auto hx: {&h0, &h1, &h2, &h3, &h012, &h013}) (*hx)[3] = 0;
230
231 ld klein_scale = binsearch(0, 10, [&] (ld d) {
232 dynamicval<eGeometry> g(geometry, elliptic ? gCell120 : geometry);
233
234 /* center of an edge */
235 hyperpoint u = C0 + (h012 + h013) * d / 2;
236
237 if(material(u) <= 0) {
238 println(hlog, "klein_scale = ", d, " bad");
239 return true;
240 }
241
242 u = normalize(u);
243
244 hyperpoint h = C0 * face;
245 for(int i=0; i<face; i++) h += d * (cspin(1, 2, M_PI*2*i/face) * h012);
246 h = normalize(h);
247
248 hyperpoint h2 = rspintox(h) * xpush0(2 * hdist0(h));
249
250 h2 = spintox(u) * h2;
251 u = spintox(u) * u;
252
253 h2 = gpushxto0(u) * h2;
254 u = gpushxto0(u) * u;
255
256 ld x = hypot(h2[1], h2[2]);
257 ld y = h2[0];
258
259 ld loop2 = 360 / (90 + atan(y/x) / degree);
260
261 println(hlog, "d=", d, " loop2= ", loop2);
262
263 if(sphere) return loop2 < loop;
264 return loop2 > loop;
265 });
266
267 /* precise ideal vertex */
268 if(klein_scale > 1-1e-5 && klein_scale < 1+1e-5) klein_scale = 1;
269
270 /* actual vertex */
271 hyperpoint v2 = C0 + klein_scale * h012;
272
273 hyperpoint midface = Hypc;
274 for(int i=0; i<face; i++) midface += cspin(1, 2, 2*i*M_PI/face) * v2;
275 midface = normalize(midface);
276 ld between_centers = 2 * hdist0(midface);
277 DEBB(DF_GEOM, ("between_centers = ", between_centers));
278
279 if(S7 == 20) {
280 spins[0] = Id;
281 spins[1] = cspin(0, 1, angle_between_faces) * cspin(1, 2, M_PI);
282 spins[2] = spins[1] * cspin(1, 2, -2 * M_PI/face) * spins[1];
283 spins[3] = spins[1] * cspin(1, 2, +2 * M_PI/face) * spins[1];
284 for(int a=4; a<10; a++) spins[a] = cspin(1, 2, 2*M_PI/face) * spins[a-3];
285 for(int a=S7/2; a<S7; a++) spins[a] = spins[a-S7/2] * cspin(0, 1, M_PI);
286 }
287
288 if(S7 == 12 || S7 == 8) {
289 spins[0] = Id;
290 spins[1] = cspin(0, 1, angle_between_faces) * cspin(1, 2, M_PI);
291 for(int a=2; a<face+1; a++) spins[a] = cspin(1, 2, 2*M_PI*(a-1)/face) * spins[1];
292 for(int a=S7/2; a<S7; a++) spins[a] = cspin(0, 1, M_PI) * spins[a-S7/2];
293 if(S7 == 8) swap(spins[6], spins[7]);
294 if(S7 == 12) swap(spins[8], spins[11]);
295 if(S7 == 12) swap(spins[9], spins[10]);
296 }
297
298 if(S7 == 6) {
299 spins[0] = Id;
300 spins[1] = cspin(0, 1, angle_between_faces) * cspin(1, 2, M_PI);
301 spins[2] = cspin(1, 2, M_PI/2) * spins[1];
302 for(int a=S7/2; a<S7; a++) spins[a] = spins[a-S7/2] * cspin(0, 1, M_PI);
303 }
304
305 if(S7 == 4) {
306 spins[0] = Id;
307 spins[1] = cspin(0, 1, angle_between_faces) * cspin(1, 2, M_PI);
308 for(int a=2; a<face+1; a++) spins[a] = cspin(1, 2, 2*M_PI*(a-1)/face) * spins[1];
309 }
310
311 cellshape.clear();
312 cellshape.resize(S7);
313 for(int a=0; a<S7; a++) {
314 for(int b=0; b<face; b++)
315 cellshape[a].push_back(spins[a] * cspin(1, 2, 2*M_PI*b/face) * v2);
316 }
317
318 cgi.adjmoves[0] = cpush(0, between_centers) * cspin(0, 2, M_PI);
319 for(int i=1; i<S7; i++) cgi.adjmoves[i] = spins[i] * cgi.adjmoves[0];
320
321 for(int a=0; a<S7; a++)
322 DEBB(DF_GEOM, ("center of ", a, " is ", tC0(cgi.adjmoves[a])));
323
324 DEBB(DF_GEOM, ("doublemove = ", tC0(cgi.adjmoves[0] * cgi.adjmoves[0])));
325
326 adjcheck = hdist(tC0(cgi.adjmoves[0]), tC0(cgi.adjmoves[1])) * 1.0001;
327
328 if(loop == 4) cgi.strafedist = adjcheck;
329 else cgi.strafedist = hdist(cgi.adjmoves[0] * C0, cgi.adjmoves[1] * C0);
330
331 if(stretch::applicable()) {
332 transmatrix T = cspin(0, 2, 90 * degree);
333 transmatrix iT = inverse(T);
334 for(auto& v: cgi.adjmoves) v = T * v * iT;
335 for(auto& vv: cellshape) for(auto& v: vv) v = T * v;
336 }
337
338 hsh.compute_hept();
339 compute_ultra();
340
341 generate_subcells();
342 }
343
generate_plain_subcubes()344 EX void generate_plain_subcubes() {
345 if(S7 != 6) throw hr_exception("generate_plain_subcubes but no cubes");
346 auto& ssh = cgi.subshapes;
347 const int sub = subcube_count;
348 auto vx = abs(cgi.heptshape->faces[0][0][0]);
349 auto vz = abs(cgi.heptshape->faces[0][0][3]);
350 for(int x=1-sub; x<sub; x+=2)
351 for(int y=1-sub; y<sub; y+=2)
352 for(int z=1-sub; z<sub; z+=2) {
353 ssh.emplace_back();
354 auto &ss = ssh.back();
355 ss.faces = cgi.heptshape->faces;
356 for(auto& face: ss.faces) for(auto& v: face) {
357 v[0] += vx * x;
358 v[1] += vx * y;
359 v[2] += vx * z;
360 v[3] += vz * (sub-1);
361 }
362 }
363 }
364
generate_coxeter(flagtype f)365 EX void generate_coxeter(flagtype f) {
366 auto& ssh = cgi.subshapes;
367 for(auto& fac: cgi.heptshape->faces) {
368 hyperpoint facectr = Hypc;
369 vector<hyperpoint> ring;
370 hyperpoint last = fac.back();
371 ring.push_back(last);
372 for(hyperpoint h: fac) {
373 if(f & cox_midedges)
374 ring.push_back(mid(last, h));
375 ring.push_back(last = h);
376 facectr += h;
377 }
378 facectr = normalize(facectr);
379
380 hyperpoint fc2 = rspintox(facectr) * xpush0(2*hdist0(facectr));
381
382 if(f & cox_vertices) {
383 for(int i=1; i<isize(ring); i++) {
384 ssh.emplace_back();
385 auto &ss = ssh.back();
386 auto h = (f & cox_othercell) ? facectr : fc2;
387 ss.faces.push_back({C0, h, ring[i-1]});
388 ss.faces.push_back({C0, h, ring[i]});
389 ss.faces.push_back({C0, ring[i-1], ring[i]});
390 ss.faces.push_back({h, ring[i-1], ring[i]});
391 }
392 }
393 else if(f & cox_midedges) {
394 ring.push_back(ring[1]);
395 for(int i=3; i<isize(ring); i+=2) {
396 ssh.emplace_back();
397 auto &ss = ssh.back();
398 auto h = (f & cox_othercell) ? facectr : fc2;
399 ss.faces.push_back({C0, ring[i-2], ring[i-1]});
400 ss.faces.push_back({C0, ring[i-1], ring[i-0]});
401 ss.faces.push_back({C0, h, ring[i-2]});
402 ss.faces.push_back({C0, h, ring[i-0]});
403 if(f & cox_othercell) {
404 ss.faces.push_back({facectr, ring[i-2], ring[i-1], ring[i-0]});
405 }
406 else {
407 ss.faces.push_back({fc2, ring[i-1], ring[i-0]});
408 ss.faces.push_back({fc2, ring[i-2], ring[i-1]});
409 }
410 }
411 }
412 else {
413 ssh.emplace_back();
414 auto &ss = ssh.back();
415 for(int i=1; i<isize(ring); i++)
416 ss.faces.push_back({C0, ring[i-1], ring[i]});
417 if(f & cox_othercell) {
418 ring.pop_back();
419 ss.faces.push_back(ring);
420 }
421 else {
422 for(int i=1; i<isize(ring); i++)
423 ss.faces.push_back({fc2, ring[i-1], ring[i]});
424 }
425 }
426 }
427 }
428
generate_special_subcubes(bool bch)429 EX void generate_special_subcubes(bool bch) {
430 if(S7 != 6) throw hr_exception("generate_plain_subcubes but no cubes");
431 const int sub = subcube_count;
432 if(1) {
433 auto vx = abs(cgi.heptshape->faces[0][0][0]);
434 auto vz = abs(cgi.heptshape->faces[0][0][3]);
435 auto step = hdist0(tC0(cgi.adjmoves[0]));
436 array<int, 3> co;
437 int s = bch ? 1 : 2;
438 for(co[0]=-sub; co[0]<=sub; co[0]+=s)
439 for(co[1]=-sub; co[1]<=sub; co[1]+=s)
440 for(co[2]=-sub; co[2]<=sub; co[2]+=s)
441 if(((co[0]^co[1]^1)&1) && ((co[0]^co[2]^1)&1)) {
442 hyperpoint ctr = Hypc;
443 ctr[3] = vz * sub;
444
445 struct direction {
446 hyperpoint dir;
447 int limit;
448 transmatrix mirror;
449 void flip() { dir = -dir; limit = 200 - limit; }
450 };
451
452 array<direction, 3> di;
453
454 int mi = 0;
455
456 for(int i=0; i<3; i++) {
457 ctr[i] += co[i] * vx;
458 auto& dii = di[i];
459 if(co[i] >= 0) {
460 dii.dir = ctangent(i, vx);
461 dii.limit = sub - co[i];
462 dii.mirror = cpush(i, +step/2) * cmirror(i) * cpush(i, -step/2);
463 }
464 else {
465 dii.dir = ctangent(i, -vx);
466 dii.limit = co[i] + sub;
467 dii.mirror = cpush(i, -step/2) * cmirror(i) * cpush(i, +step/2);
468 }
469 if(dii.limit == 0) mi++;
470 }
471
472 sort(di.begin(), di.end(), [] (direction& d1, direction& d2) { return d1.limit > d2.limit; });
473
474 cgi.subshapes.emplace_back();
475 auto &ss = cgi.subshapes.back();
476
477 auto pt0 = [&] (const array<ld, 3>& x) {
478 transmatrix M = Id;
479 hyperpoint res = ctr;
480 for(int i=0; i<3; i++) {
481 ld xx = x[i];
482 if(xx > di[i].limit) xx = 2*di[i].limit-xx, M = di[i].mirror * M;
483 res += di[i].dir * xx;
484 }
485 return normalize(M * res);
486 };
487
488 auto pt = [&] (ld x, ld y, ld z) {
489 if(sub == 1 || !bch || sphere) return pt0(make_array(x,y,z));
490
491 // Unfortunately using the rule above for bch (with sub > 1) generates faces which are not flat.
492 // Therefore, we replace the vertices by the centers of their Voronoi cells
493 // we do this only in the hyperbolic case -- it does not work correctly in the spherical case because of Voronoi not working as expected
494
495 // the arguments for pt1 are the Voronoi centers for: (x,y,z) = (1,.5,0)
496 // pt1 rearranges them to whatever (x,y,z) actually is
497
498 array<ld, 3> arg1 = {x, y, z};
499
500 auto pt1 = [&] (ld x1, ld y1, ld z1) {
501 array<ld, 3> arg0;
502 for(int i=0; i<3; i++) {
503 if(arg1[i] == 1) arg0[i] = x1;
504 else if(arg1[i] == -1) arg0[i] = -x1;
505 else if(arg1[i] == .5) arg0[i] = y1;
506 else if(arg1[i] == -.5) arg0[i] = -y1;
507 else if(arg1[i] == 0) arg0[i] = z1;
508 else throw hr_exception("unknown number in pt1");
509 }
510 return normalize(pt0(arg0));
511 };
512 hyperpoint a = pt1(0,0,0);
513 hyperpoint b = pt1(2,0,0);
514 hyperpoint c = pt1(1,1,1);
515 hyperpoint d = pt1(1,1,-1);
516 hyperpoint res = circumscribe(a, b, c, d);
517 return res;
518 };
519
520 auto add_face = [&] (std::initializer_list<hyperpoint> vh) {
521 ss.faces.push_back(vh);
522 };
523
524 const ld h = .5;
525
526 if(mi == 0) {
527 for(int s: {-1, 1}) {
528 for(int i=0; i<3; i++) {
529 if(bch)
530 add_face({pt(0,.5,s), pt(.5,0,s), pt(0,-.5,s), pt(-.5,0,s)});
531 else
532 add_face({pt(-1,-1,s), pt(-1,+1,s), pt(+1,+1,s), pt(+1,-1,s)});
533 tie(di[0], di[1], di[2]) = make_tuple(di[1], di[2], di[0]);
534 }
535 }
536 if(bch) for(int u=0; u<8; u++) {
537 for(int j=0; j<3; j++) if((u>>j)&1) di[j].flip();
538 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)});
539 for(int j=0; j<3; j++) if((u>>j)&1) di[j].flip();
540 }
541 }
542 else if(mi == 1) {
543 auto& M = di[2].mirror;
544 for(int s: {-1, 1}) {
545 if(bch)
546 add_face({pt(0,h,s), pt(h,0,s), pt(0,-h,s), pt(-h,0,s)});
547 else
548 add_face({pt(-1,-1,s), pt(-1,+1,s), pt(+1,+1,s), pt(+1,-1,s)});
549 for(int i=0; i<2; i++) {
550 if(bch)
551 add_face({pt(1,0,-.5), pt(1,-.5,0), M*pt(1,0,-.5), pt(1,.5,0)});
552 else
553 add_face({pt(-1,-1,-1), pt(-1,+1,-1), pt(-1,+1,+1), pt(-1,-1,+1)});
554 tie(di[0], di[1]) = make_tuple(di[1], di[0]); di[0].flip();
555 }
556 }
557 if(bch) for(ld s: {-1, 1}) for(int i=0; i<4; i++) {
558 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)});
559 tie(di[0], di[1]) = make_tuple(di[1], di[0]); di[0].flip();
560 }
561 }
562 else {
563 transmatrix spi = mi == 2 ? di[1].mirror * di[2].mirror : di[0].mirror * di[1].mirror;
564 if(cgi.loop == 5) spi = spi * spi;
565 vector<transmatrix> spi_power = {Id};
566 for(int i=1; i<cgi.loop; i++) spi_power.push_back(spi_power.back() * spi);
567 if(mi == 2) {
568 for(auto P: spi_power) {
569 if(bch)
570 add_face({P*pt(.5,0,-1), P*pt(0,-.5,-1), P*pt(-.5,0,-1), P*pt(0,.5,-1)});
571 else
572 add_face({P*pt(-1,-1,-1), P*pt(1,-1,-1), P*spi*pt(1,-1,-1), P*spi*pt(-1,-1,-1)});
573 }
574 vector<hyperpoint> f0, f1;
575 for(auto P: spi_power) f0.push_back(bch ? P*pt(-1,-.5,0) : P*pt(-1,-1,-1));
576 for(auto P: spi_power) f1.push_back(bch ? P*pt(+1,-.5,0) : P*pt(+1,-1,-1));
577 ss.faces.push_back(f0);
578 ss.faces.push_back(f1);
579
580 if(bch) for(auto P: spi_power) for(int s: {-1,1})
581 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)});
582 }
583 else {
584 vector<transmatrix> face_dirs = {Id};
585 for(int i=0; i<isize(face_dirs); i++)
586 for(int j=0; j<2; j++)
587 for(auto P1: spi_power) {
588 auto T = face_dirs[i];
589 if(j == 0) T = T * P1 * di[1].mirror * di[2].mirror;
590 if(j == 1) T = T * P1 * di[2].mirror * di[0].mirror;
591 bool fnd = false;
592 for(auto& F: face_dirs)
593 for(auto P: spi_power)
594 if(eqmatrix(T, F*P)) fnd = true;
595 if(!fnd) face_dirs.push_back(T);
596 }
597 // tetrahedron in {4,3,3}; dodecahedron in {4,3,5}
598 if(cgi.loop == 3) hassert(isize(face_dirs) == 4);
599 if(cgi.loop == 5) hassert(isize(face_dirs) == 12);
600 for(auto F: face_dirs) {
601 vector<hyperpoint> f0;
602 for(auto P: spi_power) f0.push_back(bch ? F*P*pt(-.5,0,-1) : F*P*pt(-1,-1,-1));
603 ss.faces.push_back(f0);
604 }
605
606 vector<transmatrix> vertex_dirs;
607 hyperpoint pter = normalize(pt0(make_array(-.5,-.5,-.5)));
608 for(auto& F: face_dirs) for(auto& P: spi_power) {
609 transmatrix T = F * P;
610 bool fnd = false;
611 for(auto T1: vertex_dirs) if(hdist(T * pter, T1*pter) < 1e-3) fnd = true;
612 if(!fnd) vertex_dirs.push_back(T);
613 }
614 if(cgi.loop == 3) hassert(isize(vertex_dirs) == 4);
615 if(cgi.loop == 5) hassert(isize(vertex_dirs) == 20);
616 if(bch) for(auto& V: vertex_dirs)
617 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)});
618 }
619 }
620 }
621 }
622 }
623
generate_bch_oct()624 EX void generate_bch_oct() {
625 if(S7 != 6) throw hr_exception("generate_bch_oct but no cubes");
626 const int sub = subcube_count;
627 if(1) {
628 auto vx = abs(cgi.heptshape->faces[0][0][0]);
629 auto vz = abs(cgi.heptshape->faces[0][0][3]);
630 array<int, 3> co;
631 // vx = 1; vz = 0;
632 for(co[0]=-sub; co[0]<sub; co[0]++)
633 for(co[1]=-sub; co[1]<sub; co[1]++)
634 for(co[2]=-sub; co[2]<sub; co[2]++) {
635 auto co1 = co;
636 array<ld, 3> sgn = {1,1,1};
637 if((co[1] ^ co[0]) & 1) co1[1]++, sgn[1] = -1;
638 if((co[2] ^ co[0]) & 1) co1[2]++, sgn[2] = -1;
639
640 hyperpoint ctr = Hypc;
641 ctr[3] = vz * sub;
642
643 auto pt = [&] (int m, ld x0, ld x1, ld x2) {
644 hyperpoint res = ctr;
645 auto x = make_array(x0, x1, x2);
646 for(int i=0; i<3; i++)
647 res[i] = vx * (co1[i] + x[(m+i)%3] * sgn[i]);
648 return res;
649 };
650
651 for(int it=0; it<2; it++) {
652 cgi.subshapes.emplace_back();
653 auto &ss = cgi.subshapes.back();
654 for(int m=0; m<3; m++) {
655 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)});
656 ss.faces.push_back({pt(m,1,0,0), pt(m,1,0,.5), pt(m,1,.5,0) });
657 }
658 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)});
659 for(int d=0; d<3; d++)
660 co1[d] += sgn[d], sgn[d] *= -1;
661 println(hlog, ss.faces);
662 }
663 }
664 }
665 }
666
generate_subcells()667 EX void generate_subcells() {
668
669 switch(variation) {
670 case eVariation::subcubes:
671 generate_plain_subcubes();
672 break;
673
674 case eVariation::dual_subcubes:
675 generate_special_subcubes(false);
676 break;
677
678 case eVariation::bch:
679 generate_special_subcubes(true);
680 break;
681
682 case eVariation::bch_oct:
683 generate_bch_oct();
684 break;
685
686 case eVariation::coxeter:
687 generate_coxeter(coxeter_param);
688 break;
689
690 case eVariation::pure: {
691 cgi.subshapes.emplace_back();
692 cgi.subshapes[0].faces = cgi.heptshape->faces;
693 break;
694 }
695
696 default:
697 throw hr_exception("unknown variation in generate_subcells");
698 }
699
700 for(auto& ss: cgi.subshapes) ss.compute_sub();
701
702 println(hlog, "subcells generated = ", isize(cgi.subshapes));
703 }
704
binary_rebase(heptagon * h,const transmatrix & V)705 void binary_rebase(heptagon *h, const transmatrix& V) {
706 }
707
708 void test();
709
710 #if HDR
711 /** \brief vertex_adjacencies[heptagon id] is a list of other heptagons which are vertex adjacent
712 * note: in case of ideal vertices this is just the face adjacency
713 **/
714 struct vertex_adjacency_info {
715 /** id of the adjacent heptagon */
716 int h_id;
717 /** transition matrix to that heptagon */
718 transmatrix T;
719 /** the sequence of moves we need to make to get there */
720 vector<int> move_sequence;
721 };
722
723 struct hrmap_closed3 : hrmap {
724 vector<heptagon*> allh;
725 vector<vector<vector<int>>> strafe_data;
726 vector<vector<transmatrix>> tmatrices;
727 vector<vector<transmatrix>> tmatrices_cell;
728 vector<cell*> acells;
729 map<cell*, pair<int, int> > local_id; /* acells index, subshape index */
730 vector<vector<cell*>> acells_by_master;
731 vector<vector<vertex_adjacency_info> > vertex_adjacencies;
732 vector<vector<vector<int>>> move_sequences;
733
adjhr::reg3::hrmap_closed3734 transmatrix adj(heptagon *h, int d) override { return tmatrices[h->fieldval][d]; }
adjhr::reg3::hrmap_closed3735 transmatrix adj(cell *c, int d) override { return tmatrices_cell[local_id.at(c).first][d]; }
736
getOriginhr::reg3::hrmap_closed3737 heptagon *getOrigin() override { return allh[0]; }
738
739 transmatrix relative_matrixc(cell *h2, cell *h1, const hyperpoint& hint) override;
740
741 void initialize(int cell_count);
allcellshr::reg3::hrmap_closed3742 vector<cell*>& allcells() override { return acells; }
743
get_cellshapehr::reg3::hrmap_closed3744 subcellshape& get_cellshape(cell *c) override {
745 if(PURE) return *cgi.heptshape ;
746 int id = local_id.at(c).second;
747 return cgi.subshapes[id];
748 }
749
master_relativehr::reg3::hrmap_closed3750 transmatrix master_relative(cell *c, bool get_inverse) override {
751 int id = local_id.at(c).second;
752 auto& ss = cgi.subshapes[id];
753 return get_inverse ? ss.from_cellcenter : ss.to_cellcenter;
754 }
755
756 void make_subconnections();
757
758 int wall_offset(cell *c) override;
shvidhr::reg3::hrmap_closed3759 int shvid(cell *c) override { return local_id.at(c).second; }
760
761 transmatrix ray_iadj(cell *c, int i) override;
762
strafehr::reg3::hrmap_closed3763 cellwalker strafe(cellwalker cw, int j) override {
764 int id = local_id.at(cw.at).first;
765 return cellwalker(cw.at->cmove(j), strafe_data[id][j][cw.spin]);
766 }
767
get_move_seqhr::reg3::hrmap_closed3768 const vector<int>& get_move_seq(cell *c, int i) override {
769 int id = local_id.at(c).first;
770 return move_sequences[id][i];
771 }
772 };
773
774 struct hrmap_quotient3 : hrmap_closed3 { };
775 #endif
776
ray_iadj(cell * c,int i)777 transmatrix hrmap_closed3::ray_iadj(cell *c, int i) {
778 if(PURE) return iadj(c, i);
779 auto& v = get_face_vertices(c, i);
780 hyperpoint h =
781 project_on_triangle(v[0], v[1], v[2]);
782 transmatrix T = rspintox(h);
783 return T * xpush(-2*hdist0(h)) * spintox(h);
784 }
785
wall_offset(cell * c)786 int hrmap_closed3::wall_offset(cell *c) {
787 if(PURE) return 0;
788 auto& wo = cgi.walloffsets[ local_id.at(c).second ];
789 if(wo.second == nullptr)
790 wo.second = c;
791 return wo.first;
792 }
793
get_face_vertices(cell * c,int d)794 EX const vector<hyperpoint>& get_face_vertices(cell *c, int d) {
795 return cgi.subshapes[currentmap->shvid(c)].faces_local[d];
796 }
797
get_face_vertex_count(cell * c,int d)798 EX int get_face_vertex_count(cell *c, int d) {
799 return isize(get_face_vertices(c, d));
800 }
801
initialize(int big_cell_count)802 void hrmap_closed3::initialize(int big_cell_count) {
803 allh.resize(big_cell_count);
804 tmatrices.resize(big_cell_count);
805 acells.clear();
806 for(int a=0; a<big_cell_count; a++) {
807 allh[a] = init_heptagon(S7);
808 allh[a]->fieldval = a;
809 }
810 }
811
812 const static bool testing_subconnections = false;
813
make_subconnections()814 void hrmap_closed3::make_subconnections() {
815 auto& ss = cgi.subshapes;
816
817 auto& vas = vertex_adjacencies;
818 vas.resize(isize(allh));
819 for(int a=0; a<isize(allh); a++) {
820 auto& va = vas[a];
821 va.emplace_back(vertex_adjacency_info{a, Id, {}});
822
823 set<unsigned> buckets;
824 for(auto& v: cgi.heptshape->vertices_only) buckets.insert(bucketer(v));
825
826 if(cgflags & qIDEAL) {
827 for(int d=0; d<S7; d++) {
828 transmatrix T = adj(allh[a], d);
829 va.emplace_back(vertex_adjacency_info{allh[a]->move(d)->fieldval, T, {d}});
830 }
831 }
832 else
833 for(int i=0; i<isize(va); i++) {
834 for(int d=0; d<S7; d++) {
835 transmatrix T = va[i].T * adj(allh[va[i].h_id], d);
836 bool found = false;
837 for(auto& va2: va) if(eqmatrix(va2.T, T)) found = true;
838 if(found) continue;
839
840 bool found_va = false;
841 for(auto& w: cgi.heptshape->vertices_only)
842 if(buckets.count(bucketer(T*w)))
843 found_va = true;
844 if(!found_va) continue;
845 va.emplace_back(vertex_adjacency_info{allh[va[i].h_id]->move(d)->fieldval, T, va[i].move_sequence});
846 va.back().move_sequence.push_back(d);
847 }
848 }
849 }
850
851 map<int, int> by_sides;
852
853 vector<map<unsigned, vector<hyperpoint> > > which_cell_0;
854 which_cell_0.resize(isize(allh));
855
856 acells_by_master.resize(isize(allh));
857 for(int a=0; a<isize(allh); a++) {
858 for(int id=0; id<isize(ss); id++) {
859 bool exists = false;
860 auto& cc = ss[id].cellcenter;
861 for(auto& va: vertex_adjacencies[a]) {
862 hyperpoint h = iso_inverse(va.T) * cc;
863 for(auto h1: which_cell_0[va.h_id][bucketer(h)])
864 if(hdist(h1, h) < 1e-6)
865 exists = true;
866 }
867 if(exists) continue;
868 cell *c = newCell(isize(ss[id].faces), allh[a]);
869 by_sides[isize(ss[id].faces)]++;
870 if(!allh[a]->c7)
871 allh[a]->c7 = c;
872 local_id[c] = {isize(acells), id};
873 acells.push_back(c);
874 acells_by_master[a].push_back(c);
875 which_cell_0[a][bucketer(cc)].push_back(cc);
876 }
877 }
878
879 println(hlog, "found ", isize(acells), " cells, ", by_sides);
880
881 tmatrices_cell.resize(isize(acells));
882 move_sequences.resize(isize(acells));
883 int failures = 0;
884
885 vector<map<unsigned, vector<pair<cell*, int> > > > which_cell;
886 which_cell.resize(isize(allh));
887
888 for(cell *c: acells) {
889 int id = local_id[c].second;
890 for(int i=0; i<c->type; i++)
891 which_cell[c->master->fieldval][bucketer(ss[id].face_centers[i])].emplace_back(c, i);
892 }
893
894 strafe_data.resize(isize(acells));
895
896 for(cell *c: acells) {
897 int id = local_id[c].second;
898 int cid = local_id[c].first;
899 auto& tmcell = tmatrices_cell[cid];
900 vector<int> foundtab;
901 vector<tuple<int, int, int>> foundtab_ids;
902 strafe_data[cid].resize(c->type);
903 for(int i=0; i<c->type; i++) {
904 int found = 0;
905 hyperpoint ctr = ss[id].face_centers[i];
906 transmatrix T1 = Id;
907 int h_id = c->master->fieldval;
908 vector<int> path;
909 while(true) {
910 int d = -1;
911 ld dist = hdist0(ctr);
912 for(int d1=0; d1<S7; d1++) {
913 auto ctr1 = iso_inverse(tmatrices[h_id][d1]) * ctr;
914 ld dist1 = hdist0(ctr1);
915 if(dist1 < dist - 1e-6) d = d1, dist = dist1;
916 }
917 if(d == -1) break;
918 path.push_back(d);
919 T1 = T1 * tmatrices[h_id][d];
920 ctr = iso_inverse(tmatrices[h_id][d]) * ctr;
921 h_id = allh[h_id]->move(d)->fieldval;
922 }
923
924 for(auto& va: vertex_adjacencies[h_id]) {
925 hyperpoint ctr1 = iso_inverse(va.T) * ctr;
926 auto bucket = bucketer(ctr1);
927 for(auto p: which_cell[va.h_id][bucket]) {
928 cell *c1 = p.first;
929 int j = p.second;
930 int id1 = local_id[c1].second;
931 if(hdist(ctr1, ss[id1].face_centers[j]) < 1e-6) {
932 transmatrix T2 = T1 * va.T;
933 if(id == id1 && eqmatrix(T2, Id)) continue;
934 c->c.connect(i, c1, j, false);
935 if(!found) {
936 tmcell.push_back(ss[id].from_cellcenter * T2 * ss[id1].to_cellcenter);
937 if(elliptic) fixelliptic(tmcell.back());
938
939 auto& ms = move_sequences[local_id[c].first];
940 ms.push_back(path);
941 for(auto dir: va.move_sequence) ms.back().push_back(dir);
942
943 auto& sd = strafe_data[cid][i];
944 sd.resize(c->type, -1);
945
946 for(int i1=0; i1<c->type; i1++) {
947 set<unsigned> facevertices;
948 for(auto v: ss[id].faces[i1]) facevertices.insert(bucketer(v));
949 if(ss[id].dirdist[i][i1] == 1) {
950 int found_strafe = 0;
951 for(int j1=0; j1<c1->type; j1++) if(j1 != j) {
952 int num = 0;
953 for(auto v: ss[id1].faces[j1])
954 if(facevertices.count(bucketer(T2*v)))
955 num++;
956 if(num == 2) sd[i1] = j1, found_strafe++;
957 }
958 if(found_strafe != 1) println(hlog, "found_strafe = ", found_strafe);
959 }
960 }
961
962 /* for bch, also provide second-order strafe */
963 if(variation == eVariation::bch) for(int i1=0; i1<c->type; i1++) {
964 if(ss[id].dirdist[i][i1] != 2) continue;
965 if(isize(ss[id].faces[i]) == 6) continue;
966 if(isize(ss[id].faces[i1]) == 6) continue;
967 vector<int> fac;
968 for(int i2=0; i2<c->type; i2++) if(ss[id].dirdist[i][i2] == 1 && ss[id].dirdist[i2][i1] == 1)
969 fac.push_back(sd[i2]);
970 if(isize(fac) != 2) {
971 println(hlog, "fac= ", fac);
972 throw hr_exception("fac error");
973 }
974 int found_strafe = 0;
975 for(int j1=0; j1<c1->type; j1++) if(j1 != j)
976 if(ss[id1].dirdist[j1][fac[0]] == 1)
977 if(ss[id1].dirdist[j1][fac[1]] == 1) {
978 sd[i1] = j1;
979 found_strafe++;
980 }
981 if(found_strafe != 1) println(hlog, "found_strafe = ", found_strafe, " (second order)");
982 }
983 }
984 foundtab_ids.emplace_back(va.h_id, id1, j);
985 found++;
986 }
987 }
988 if(found && !testing_subconnections) break;
989 }
990 if(testing_subconnections && !found) {
991 c->c.connect(i, c, i, false);
992 tmcell.push_back(Id);
993 }
994 foundtab.push_back(found);
995 if(found != 1) failures++;
996 }
997 if(testing_subconnections) println(hlog, "foundtab = ", foundtab);
998 }
999 println(hlog, "total failures = ", failures);
1000 if(failures && !testing_subconnections) throw hr_exception("hrmap_closed3 failures");
1001 }
1002
relative_matrixc(cell * c2,cell * c1,const hyperpoint & hint)1003 transmatrix hrmap_closed3::relative_matrixc(cell *c2, cell *c1, const hyperpoint& hint) {
1004 if(c1 == c2) return Id;
1005 int d = hr::celldistance(c2, c1);
1006
1007 for(int a=0; a<c1->type; a++) if(hr::celldistance(c2, c1->move(a)) < d)
1008 return adj(c1, a) * relative_matrix(c2, c1->move(a), hint);
1009
1010 for(int a=0; a<c1->type; a++) println(hlog, "d=", d, " vs ", hr::celldistance(c2, c1->move(a)));
1011
1012 println(hlog, "error in hrmap_quotient3:::relative_matrix");
1013 return Id;
1014 }
1015
1016 #if CAP_CRYSTAL
encode_coord(const crystal::coord & co)1017 int encode_coord(const crystal::coord& co) {
1018 int c = 0;
1019 for(int i=0; i<4; i++) c |= ((co[i]>>1) & 3) << (2*i);
1020 return c;
1021 }
1022
decode_coord(int a)1023 EX crystal::coord decode_coord(int a) {
1024 crystal::coord co;
1025 for(int i=0; i<4; i++) co[i] = (a & 3) * 2, a >>= 2;
1026 return co;
1027 }
1028
1029 struct hrmap_from_crystal : hrmap_quotient3 {
1030
hrmap_from_crystalhr::reg3::hrmap_from_crystal1031 hrmap_from_crystal() {
1032 initialize(256);
1033 if(1) {
1034 auto m = crystal::new_map();
1035 dynamicval<hrmap*> cm(currentmap, m);
1036 for(int a=0; a<256; a++) {
1037 auto co = decode_coord(a);
1038 heptagon *h1 = get_heptagon_at(co);
1039 for(int d=0; d<8; d++) {
1040 int b = encode_coord(crystal::get_coord(h1->cmove(d)));
1041 allh[a]->c.connect(d, allh[b], h1->c.spin(d), false);
1042 tmatrices[a].push_back(crystal::get_adj(h1, d));
1043 }
1044 }
1045 delete m;
1046 }
1047 }
1048 };
1049 #endif
1050
1051 struct hrmap_field3 : reg3::hrmap_quotient3 {
1052
1053 fieldpattern::fpattern *f;
1054
hrmap_field3hr::reg3::hrmap_field31055 hrmap_field3(fieldpattern::fpattern *ptr) {
1056
1057 f = ptr;
1058
1059 auto lgr = f->local_group;
1060
1061 int N = isize(f->matrices) / lgr;
1062 initialize(N);
1063
1064 vector<int> moveid(S7), movedir(lgr);
1065 for(int s=0; s<lgr; s++)
1066 for(int i=0; i<S7; i++) if(eqmatrix(f->fullv[s] * cgi.adjmoves[0], cgi.adjmoves[i]))
1067 moveid[i] = s;
1068
1069 for(int s=0; s<lgr; s++)
1070 for(int i=0; i<S7; i++) if(hdist(tC0(inverse(f->fullv[s]) * cgi.adjmoves[0]), tC0(cgi.adjmoves[i])) < 1e-4)
1071 movedir[s] = i;
1072
1073 for(int a=0; a<N; a++) {
1074 tmatrices[a].resize(S7);
1075 for(int b=0; b<S7; b++) {
1076 int k = lgr*a;
1077 k = f->matcode[ f->mmul(f->mmul(f->matrices[k], f->matrices[moveid[b]]), f->P) ];
1078 for(int l=0; l<lgr; l++) if(f->gmul(k, l) % lgr == 0) {
1079 tmatrices[a][b] = cgi.adjmoves[b] * f->fullv[l];
1080 allh[a]->c.connect(b, allh[k/lgr], movedir[l], false);
1081 }
1082 }
1083 }
1084 make_subconnections();
1085 create_patterns();
1086 }
1087
1088 set<cellwalker> plane;
1089
make_planehr::reg3::hrmap_field31090 void make_plane(cellwalker cw) {
1091 if(plane.count(cw)) return;
1092 plane.insert(cw);
1093 auto& ss = get_cellshape(cw.at);
1094 for(int i=0; i<cw.at->type; i++)
1095 if(ss.dirdist[i][cw.spin] == 1)
1096 make_plane(strafe(cw, i));
1097 }
1098
1099
create_patternshr::reg3::hrmap_field31100 void create_patterns() {
1101 DEBB(DF_GEOM, ("creating pattern = ", isize(allh)));
1102
1103 if(!PURE) {
1104 println(hlog, "create_patterns not implemented");
1105 return;
1106 }
1107
1108 // also, strafe needs currentmap
1109 dynamicval<hrmap*> c(currentmap, this);
1110
1111 if(S7 == 12) {
1112 // Emerald in 534
1113 cell *a = gamestart();
1114 cell *b = a;
1115 for(cell *c: allcells())
1116 if(bounded_celldistance(a, c) == 5) {
1117 b = c;
1118 break;
1119 }
1120 for(cell *c: allcells())
1121 if(bounded_celldistance(a, c) > bounded_celldistance(b, c))
1122 c->master->zebraval |= 1;
1123
1124 // Vineyard in 534
1125 b = (cellwalker(a, 0) + wstep + rev + wstep).at;
1126 for(cell *c: allcells())
1127 if(bounded_celldistance(a, c) == bounded_celldistance(b, c))
1128 c->master->zebraval |= 2;
1129 }
1130
1131 if(S7 == 6 && ginf[geometry].vertex == 5) {
1132 // Emerald in 534
1133 cell *a = gamestart();
1134 for(cell *c: allcells())
1135 if(bounded_celldistance(a, c) > 3)
1136 c->master->zebraval |= 1;
1137
1138 // Vineyard in 435
1139 make_plane(cellwalker(gamestart(), 0));
1140 DEBB(DF_GEOM, ("plane size = ", isize(plane)));
1141
1142 set<int> plane_indices;
1143 for(auto cw: plane) plane_indices.insert(cw.at->master->fieldval);
1144
1145 int fN = isize(f->matrices);
1146
1147 set<int> nwi;
1148 for(int i=0; i<fN; i++) {
1149 bool ok = true;
1150 for(auto o: plane_indices) {
1151 int j = f->gmul(i, o * f->local_group) / f->local_group;
1152 if(plane_indices.count(j)) ok = false;
1153 forCellEx(c1, allcells()[j]) if(plane_indices.count(c1->master->fieldval)) ok = false;
1154 }
1155 if(ok) nwi.insert(i);
1156 }
1157
1158 int gpow = 0;
1159
1160 for(int i: nwi) {
1161 int pw = 1;
1162 int at = i;
1163 while(true) {
1164 at = f->gmul(at, i);
1165 if(!nwi.count(at)) break;
1166 pw++;
1167 }
1168 if(pw == 4) gpow = i;
1169 }
1170
1171 int u = 0;
1172 for(int a=0; a<5; a++) {
1173 for(int o: plane_indices) {
1174 int j = f->gmul(u, o * f->local_group) / f->local_group;
1175 allcells()[j]->master->zebraval |= 2;
1176 }
1177 u = f->gmul(u, gpow);
1178 }
1179 }
1180 }
1181 };
1182
1183 /** \brief homology cover of the Seifert-Weber space */
1184 namespace seifert_weber {
1185
1186 using crystal::coord;
1187
1188 vector<coord> periods;
1189
flip(int x)1190 int flip(int x) { return (x+6) % 12; }
1191
build_reps()1192 void build_reps() {
1193 // start_game();
1194 auto& hsh = get_hsh();
1195
1196 set<coord> boundaries;
1197
1198 for(int a=0; a<12; a++)
1199 for(int b=0; b<12; b++) if(hsh.dirdist[a][b] == 1) {
1200 coord res = crystal::c0;
1201 int sa = a, sb = b;
1202 do {
1203 // printf("%d ", sa);
1204 if(sa < 6) res[sa]++; else res[sa-6]--;
1205 sa = flip(sa);
1206 sb = flip(sb);
1207 swap(sa, sb);
1208 sb = hsh.next_dir[sa][sb];
1209 // sb = next_dirsa][sb];
1210 }
1211 while(a != sa || b != sb);
1212 // printf("\n");
1213 if(res > crystal::c0)
1214 boundaries.insert(res);
1215 }
1216
1217 periods.clear();
1218
1219 for(int index = 5; index >= 0; index--) {
1220 for(auto k: boundaries) println(hlog, k);
1221 DEBB(DF_GEOM, ("simplifying..."));
1222
1223 for(auto by: boundaries) if(among(by[index], 1, -1)) {
1224 DEBB(DF_GEOM, ("simplifying by ", by));
1225 periods.push_back(by);
1226 set<coord> nb;
1227
1228 for(auto v: boundaries)
1229 if(v == by) ;
1230 else if(v[index] % by[index] == 0)
1231 nb.insert(v - by * (v[index] / by[index]));
1232 else println(hlog, "error");
1233
1234 boundaries = move(nb);
1235 break;
1236 }
1237 }
1238 }
1239
get_rep(coord a)1240 int get_rep(coord a) {
1241 a = a - periods[0] * (a[5] / periods[0][5]);
1242 a = a - periods[1] * (a[4] / periods[1][4]);
1243 a = a - periods[2] * (a[3] / periods[2][3]);
1244 for(int i=0; i<3; i++) a[i] = gmod(a[i], 5);
1245 return a[2] * 25 + a[1] * 5 + a[0];
1246 }
1247
decode(int id)1248 coord decode(int id) {
1249 coord res = crystal::c0;
1250 for(int a=0; a<3; a++) res[a] = id % 5, id /= 5;
1251 return res;
1252 }
1253
1254 struct hrmap_singlecell : hrmap_quotient3 {
hrmap_singlecellhr::reg3::seifert_weber::hrmap_singlecell1255 hrmap_singlecell(ld angle) {
1256 initialize(1);
1257 tmatrices[0].resize(S7);
1258 for(int b=0; b<S7; b++) {
1259 allh[0]->c.connect(b, allh[0], (b+S7/2) % S7, false);
1260 transmatrix T = cgi.adjmoves[b];
1261 hyperpoint p = tC0(T);
1262 tmatrices[0][b] = rspintox(p) * xpush(hdist0(p)) * cspin(2, 1, angle) * spintox(p);
1263 }
1264 }
1265 };
1266
1267 struct hrmap_seifert_cover : hrmap_quotient3 {
1268
hrmap_seifert_coverhr::reg3::seifert_weber::hrmap_seifert_cover1269 hrmap_seifert_cover() {
1270 if(periods.empty()) build_reps();
1271 initialize(125);
1272 for(int a=0; a<125; a++) {
1273 tmatrices[a].resize(12);
1274 for(int b=0; b<12; b++) {
1275 coord x = decode(a);
1276 if(b < 6) x[b]++;
1277 else x[b-6]--;
1278 int a1 = get_rep(x);
1279 allh[a]->c.connect(b, allh[a1], flip(b), false);
1280 transmatrix T = cgi.adjmoves[b];
1281 hyperpoint p = tC0(T);
1282 tmatrices[a][b] = rspintox(p) * xpush(hdist0(p)) * cspin(2, 1, 108 * degree) * spintox(p);
1283 }
1284 }
1285 }
1286 };
1287
1288 }
1289
1290 struct hrmap_h3 : hrmap {
1291
1292 heptagon *origin;
1293 hrmap *binary_map;
1294 hrmap_quotient3 *quotient_map;
1295
1296 map<heptagon*, pair<heptagon*, transmatrix>> reg_gmatrix;
1297 map<heptagon*, vector<pair<heptagon*, transmatrix> > > altmap;
1298
allcellshr::reg3::hrmap_h31299 vector<cell*>& allcells() override {
1300 return hrmap::allcells();
1301 }
1302
hrmap_h3hr::reg3::hrmap_h31303 hrmap_h3() {
1304 origin = init_heptagon(S7);
1305 heptagon& h = *origin;
1306 h.s = hsOrigin;
1307 h.c7 = newCell(S7, origin);
1308 worst_error1 = 0, worst_error2 = 0;
1309
1310 dynamicval<hrmap*> cr(currentmap, this);
1311
1312 heptagon *alt = NULL;
1313 transmatrix T = Id;
1314
1315 binary_map = nullptr;
1316 quotient_map = nullptr;
1317
1318 #if CAP_FIELD
1319 #if CAP_CRYSTAL
1320 if(geometry == gSpace344) {
1321 quotient_map = new hrmap_from_crystal;
1322 }
1323 else
1324 #endif
1325 if(geometry == gSpace535) {
1326 quotient_map = new seifert_weber::hrmap_seifert_cover;
1327 }
1328 else if(hyperbolic) {
1329 quotient_map = new hrmap_field3(&currfp);
1330 }
1331 #endif
1332 h.zebraval = quotient_map ? quotient_map->allh[0]->zebraval : 0;
1333
1334 #if CAP_BT
1335 if(hyperbolic) {
1336 dynamicval<eGeometry> g(geometry, gBinary3);
1337 bt::build_tmatrix();
1338 alt = init_heptagon(S7);
1339 alt->s = hsOrigin;
1340 alt->alt = alt;
1341 binary_map = bt::new_alt_map(alt);
1342 T = xpush(.01241) * spin(1.4117) * xpush(0.1241) * cspin(0, 2, 1.1249) * xpush(0.07) * Id;
1343 }
1344 #endif
1345
1346 reg_gmatrix[origin] = make_pair(alt, T);
1347 altmap[alt].emplace_back(origin, T);
1348
1349 celllister cl(origin->c7, 4, 100000, NULL);
1350 for(cell *c: cl.lst) {
1351 hyperpoint h = tC0(relative_matrix(c->master, origin, C0));
1352 cgi.close_distances[bucketer(h)] = cl.getdist(c);
1353 }
1354 }
1355
1356 ld worst_error1, worst_error2;
1357
getOriginhr::reg3::hrmap_h31358 heptagon *getOrigin() override {
1359 return origin;
1360 }
1361
fix_distanceshr::reg3::hrmap_h31362 void fix_distances(heptagon *h, heptagon *h2) {
1363 vector<heptagon*> to_fix;
1364
1365 auto fix_pair = [&] (heptagon *h, heptagon *h2) {
1366 if(!h2) return;
1367 if(h->distance > h2->distance+1) {
1368 h->distance = h2->distance + 1;
1369 to_fix.push_back(h);
1370 }
1371 else if(h2->distance > h->distance+1) {
1372 h2->distance = h->distance + 1;
1373 to_fix.push_back(h2);
1374 }
1375 if(h->alt && h->alt == h2->alt) {
1376 if(altdist(h) > altdist(h2) + 1) {
1377 altdist(h) = altdist(h2) + 1;
1378 to_fix.push_back(h);
1379 }
1380 else if (altdist(h2) > altdist(h) + 1) {
1381 altdist(h2) = altdist(h) + 1;
1382 to_fix.push_back(h2);
1383 }
1384 }
1385 };
1386
1387 if(!h2) to_fix = {h};
1388 else fix_pair(h, h2);
1389
1390 for(int i=0; i<isize(to_fix); i++) {
1391 h = to_fix[i];
1392 for(int j=0; j<S7; j++) fix_pair(h, h->move(j));
1393 }
1394 }
1395
1396 #define DEB 0
1397
counterparthr::reg3::hrmap_h31398 heptagon *counterpart(heptagon *h) {
1399 return quotient_map->allh[h->fieldval];
1400 }
1401
verify_neighborshr::reg3::hrmap_h31402 void verify_neighbors(heptagon *alt, int steps, const hyperpoint& hT) {
1403 ld err;
1404 for(auto& p2: altmap[alt]) if((err = intval(tC0(p2.second), hT)) < 1e-3) {
1405 println(hlog, "FAIL");
1406 exit(3);
1407 }
1408 #if CAP_BT
1409 if(steps) {
1410 dynamicval<eGeometry> g(geometry, gBinary3);
1411 dynamicval<hrmap*> cm(currentmap, binary_map);
1412 for(int i=0; i<alt->type; i++)
1413 verify_neighbors(alt->cmove(i), steps-1, currentmap->iadj(alt, i) * hT);
1414 }
1415 #endif
1416 }
1417
create_stephr::reg3::hrmap_h31418 heptagon *create_step(heptagon *parent, int d) override {
1419 auto& p1 = reg_gmatrix[parent];
1420 if(DEB) println(hlog, "creating step ", parent, ":", d, ", at ", p1.first, tC0(p1.second));
1421 heptagon *alt = p1.first;
1422 #if CAP_FIELD
1423 transmatrix T = p1.second * (quotient_map ? quotient_map->tmatrices[parent->fieldval][d] : cgi.adjmoves[d]);
1424 #else
1425 transmatrix T = p1.second * cgi.adjmoves[d];
1426 #endif
1427 #if CAP_BT
1428 if(hyperbolic) {
1429 dynamicval<eGeometry> g(geometry, gBinary3);
1430 dynamicval<hrmap*> cm(currentmap, binary_map);
1431 binary_map->virtualRebase(alt, T);
1432 }
1433 #endif
1434
1435 fixmatrix(T);
1436 auto hT = tC0(T);
1437
1438 if(DEB) println(hlog, "searching at ", alt, ":", hT);
1439
1440 if(DEB) for(auto& p2: altmap[alt]) println(hlog, "for ", tC0(p2.second), " intval is ", intval(tC0(p2.second), hT));
1441
1442 ld err;
1443
1444 for(auto& p2: altmap[alt]) if((err = intval(tC0(p2.second), hT)) < 1e-3) {
1445 if(err > worst_error1) println(hlog, format("worst_error1 = %lg", double(worst_error1 = err)));
1446 // println(hlog, "YES found in ", isize(altmap[alt]));
1447 if(DEB) println(hlog, "-> found ", p2.first);
1448 int fb = 0;
1449 hyperpoint old = tC0(p1.second);;
1450 #if CAP_FIELD
1451 if(quotient_map) {
1452 p2.first->c.connect(counterpart(parent)->c.spin(d), parent, d, false);
1453 fix_distances(p2.first, parent);
1454 return p2.first;
1455 }
1456 #endif
1457 for(int d2=0; d2<S7; d2++) {
1458 hyperpoint back = p2.second * tC0(cgi.adjmoves[d2]);
1459 if((err = intval(back, old)) < 1e-3) {
1460 if(err > worst_error2) println(hlog, format("worst_error2 = %lg", double(worst_error2 = err)));
1461 if(p2.first->move(d2)) println(hlog, "error: repeated edge");
1462 p2.first->c.connect(d2, parent, d, false);
1463 fix_distances(p2.first, parent);
1464 fb++;
1465 }
1466 }
1467 if(fb != 1) {
1468 println(hlog, "found fb = ", fb);
1469 println(hlog, old);
1470 for(int d2=0; d2<S7; d2++) {
1471 println(hlog, p2.second * tC0(cgi.adjmoves[d2]), " in distance ", intval(p2.second * tC0(cgi.adjmoves[d2]), old));
1472 }
1473 parent->c.connect(d, parent, d, false);
1474 return parent;
1475 }
1476 return p2.first;
1477 }
1478
1479 if(extra_verification) verify_neighbors(alt, extra_verification, hT);
1480
1481 if(DEB) println(hlog, "-> not found");
1482 int d2 = 0, fv = isize(reg_gmatrix);
1483 #if CAP_FIELD
1484 if(quotient_map) {
1485 auto cp = counterpart(parent);
1486 d2 = cp->c.spin(d);
1487 fv = cp->c.move(d)->fieldval;
1488 }
1489 #endif
1490 heptagon *created = init_heptagon(S7);
1491 created->c7 = newCell(S7, created);
1492 #if CAP_FIELD
1493 if(quotient_map) {
1494 created->emeraldval = fv;
1495 created->zebraval = quotient_map->allh[fv]->zebraval;
1496 }
1497 else
1498 #endif
1499 created->zebraval = hrand(10);
1500 created->fieldval = fv;
1501 created->distance = parent->distance + 1;
1502 created->fiftyval = 9999;
1503 fixmatrix(T);
1504 reg_gmatrix[created] = make_pair(alt, T);
1505 altmap[alt].emplace_back(created, T);
1506 created->c.connect(d2, parent, d, false);
1507 return created;
1508 }
1509
~hrmap_h3hr::reg3::hrmap_h31510 ~hrmap_h3() {
1511 #if CAP_BT
1512 if(binary_map) {
1513 dynamicval<eGeometry> g(geometry, gBinary3);
1514 delete binary_map;
1515 }
1516 #endif
1517 if(quotient_map) delete quotient_map;
1518 clearfrom(origin);
1519 }
1520
1521 map<heptagon*, int> reducers;
1522
link_althr::reg3::hrmap_h31523 bool link_alt(heptagon *h, heptagon *alt, hstate firststate, int dir) override {
1524 altdist(h) = 0;
1525 if(firststate != hsOrigin) reducers[h] = dir;
1526 return true;
1527 }
1528
extend_altmaphr::reg3::hrmap_h31529 void extend_altmap(heptagon* h, int levs, bool link_cdata) override {
1530 if(reducers.count(h)) {
1531 heptspin hs(h, reducers[h]);
1532 reducers.erase(h);
1533 hs += wstep;
1534 hs += rev;
1535 altdist(hs.at) = altdist(h) - 1;
1536 hs.at->alt = h->alt;
1537 reducers[hs.at] = hs.spin;
1538 fix_distances(hs.at, NULL);
1539 }
1540 for(int i=0; i<S7; i++) {
1541 auto h2 = h->cmove(i);
1542 if(h2->alt == NULL) {
1543 h2->alt = h->alt;
1544 altdist(h2) = altdist(h) + 1;
1545 fix_distances(h2, NULL);
1546 }
1547 }
1548 }
1549
adjhr::reg3::hrmap_h31550 transmatrix adj(heptagon *h, int d) override {
1551 #if CAP_FIELD
1552 if(quotient_map) return quotient_map->adj(h, d);
1553 else
1554 #endif
1555 return relative_matrix(h->cmove(d), h, C0);
1556 }
1557
relative_matrixhhr::reg3::hrmap_h31558 transmatrix relative_matrixh(heptagon *h2, heptagon *h1, const hyperpoint& hint) override {
1559 auto p1 = reg_gmatrix[h1];
1560 auto p2 = reg_gmatrix[h2];
1561 transmatrix T = Id;
1562 #if CAP_BT
1563 if(hyperbolic) {
1564 dynamicval<eGeometry> g(geometry, gBinary3);
1565 dynamicval<hrmap*> cm(currentmap, binary_map);
1566 T = binary_map->relative_matrix(p2.first, p1.first, hint);
1567 }
1568 #endif
1569 T = inverse(p1.second) * T * p2.second;
1570 if(elliptic && T[LDIM][LDIM] < 0) T = centralsym * T;
1571 return T;
1572 }
1573
get_cellshapehr::reg3::hrmap_h31574 subcellshape& get_cellshape(cell *c) override {
1575 return *cgi.heptshape;
1576 }
1577
strafehr::reg3::hrmap_h31578 cellwalker strafe(cellwalker cw, int j) override {
1579 hyperpoint hfront = tC0(cgi.adjmoves[cw.spin]);
1580 cw.at->cmove(j);
1581 transmatrix T = currentmap->adj(cw.at, j);
1582 for(int i=0; i<S7; i++) if(i != cw.at->c.spin(j))
1583 if(hdist(hfront, T * tC0(cgi.adjmoves[i])) < cgi.strafedist + .01)
1584 return cellwalker(cw.at->cmove(j), i);
1585 throw hr_exception("incorrect strafe");
1586 }
1587
1588 };
1589
1590 struct hrmap_sphere3 : hrmap_closed3 {
1591
1592 vector<transmatrix> locations;
1593
hrmap_sphere3hr::reg3::hrmap_sphere31594 hrmap_sphere3() {
1595 heptagon *h = init_heptagon(S7);
1596 h->s = hsOrigin;
1597
1598 locations.push_back(Id);
1599 allh.push_back(h);
1600
1601 for(int i=0; i<isize(allh); i++) {
1602 tmatrices.emplace_back();
1603 auto& tmi = tmatrices.back();
1604 transmatrix T1 = locations[i];
1605 hyperpoint old = tC0(T1);
1606 for(int d=0; d<S7; d++) {
1607 transmatrix T = T1 * cgi.adjmoves[d];
1608 fixmatrix(T);
1609 auto hT = tC0(T);
1610
1611 bool hopf = stretch::applicable();
1612
1613 if(hopf)
1614 T = stretch::translate(hT);
1615
1616 for(int i1=0; i1<isize(allh); i1++)
1617 if(intval(tC0(locations[i1]), hT) < 1e-3) {
1618 int fb = 0;
1619 for(int d2=0; d2<S7; d2++) {
1620 hyperpoint back = locations[i1] * tC0(cgi.adjmoves[d2]);
1621 if(intval(back, old) < 1e-3) {
1622 allh[i]->c.connect(d, allh[i1], d2, false);
1623 fb++;
1624 tmi.push_back(inverse(T1) * locations[i1]);
1625 }
1626 }
1627 if(fb != 1) throw hr_exception("friend not found");
1628 goto next_d;
1629 }
1630
1631 if(1) {
1632 int d2 = 0;
1633
1634 if(hopf) {
1635 for(d2=0; d2<S7; d2++) {
1636 hyperpoint back = T * tC0(cgi.adjmoves[d2]);
1637 if(intval(back, old) < 1e-3)
1638 break;
1639 }
1640 if(d2 == S7)
1641 throw hr_exception("Hopf connection failed");
1642 }
1643
1644 heptagon *h = init_heptagon(S7);
1645 h->zebraval = hrand(10);
1646 h->fieldval = isize(allh);
1647 h->fiftyval = 9999;
1648 allh.push_back(h);
1649 locations.push_back(T);
1650 if(isnan(T[0][0])) exit(1);
1651
1652 allh[i]->c.connect(d, h, d2, false);
1653 tmi.push_back(inverse(T1) * T);
1654 if(elliptic) fixelliptic(tmi.back());
1655 }
1656 next_d: ;
1657 }
1658 }
1659
1660 make_subconnections();
1661 }
1662
~hrmap_sphere3hr::reg3::hrmap_sphere31663 ~hrmap_sphere3() {
1664 clearfrom(allh[0]);
1665 }
1666
relative_matrixhhr::reg3::hrmap_sphere31667 transmatrix relative_matrixh(heptagon *h2, heptagon *h1, const hyperpoint& hint) override {
1668 return iso_inverse(locations[h1->fieldval]) * locations[h2->fieldval];
1669 }
1670
1671 };
1672
get_sphere_loc(int v)1673 EX const transmatrix& get_sphere_loc(int v) {
1674 return ((hrmap_sphere3*)currentmap)->locations[v];
1675 }
1676
1677 struct hrmap_h3_rule : hrmap {
1678
1679 heptagon *origin;
1680 reg3::hrmap_quotient3 *quotient_map;
1681 reg3::hrmap_quotient3 *emerald_map;
1682
1683 fieldpattern::fpattern fp;
1684
1685 vector<int> root;
1686 string other;
1687 vector<short> children;
1688
1689 vector<int> otherpos;
1690
load_rulesethr::reg3::hrmap_h3_rule1691 void load_ruleset(string fname) {
1692 string buf;
1693 #if ISANDROID || ISIOS
1694 buf = get_asset(fname);
1695 #else
1696 FILE *f = fopen(fname.c_str(), "rb");
1697 if(!f) f = fopen((rsrcdir + fname).c_str(), "rb");
1698 buf.resize(1000000);
1699 int qty = fread(&buf[0], 1, 1000000, f);
1700 buf.resize(qty);
1701 fclose(f);
1702 #endif
1703
1704 shstream ins(decompress_string(buf));
1705 dynamicval<bool> q(fieldpattern::use_quotient_fp, true);
1706 hread_fpattern(ins, fp);
1707
1708 hread(ins, root);
1709 hread(ins, children);
1710 hread(ins, other);
1711 }
1712
1713 /** \brief address = (fieldvalue, state) */
1714 typedef pair<int, int> address;
1715
1716 /** nles[x] lists the addresses from which we can reach address x
1717 * without ever ending in the starting point */
1718
1719 map<address, set<address>> nonlooping_earlier_states;
1720
1721 vector<vector<int>> possible_states;
1722
find_mappingshr::reg3::hrmap_h3_rule1723 void find_mappings() {
1724 auto &nles = nonlooping_earlier_states;
1725 nles.clear();
1726 vector<address> bfs;
1727 int qty = isize(quotient_map->allh);
1728 if(geometry == gSpace535) qty = 1;
1729 for(int i=0; i<qty; i++)
1730 bfs.emplace_back(i, root[i]);
1731 auto mov = [&] (int fv, int d) {
1732 if(geometry == gSpace535) return 0;
1733 return quotient_map->allh[fv]->move(d)->fieldval;
1734 };
1735 int qstate = isize(children) / S7;
1736 DEBB(DF_GEOM, ("qstate = ", qstate));
1737 for(int i=0; i<isize(bfs); i++) {
1738 address last = bfs[i];
1739 int state = last.second;
1740 int fv = last.first;
1741 for(int d=0; d<S7; d++) {
1742 int nstate = children[state*S7+d];
1743 if(nstate < -1) nstate += (1<<16);
1744 if(nstate >= 0) {
1745 address next = {mov(fv, d), nstate};
1746 if(!nles.count(next)) bfs.push_back(next);
1747 nles[next].insert(last);
1748 }
1749 }
1750 }
1751
1752 vector<int> q(qstate, 0);
1753 for(auto p: bfs) q[p.second]++;
1754 vector<int> q2(isize(quotient_map->allh)+1, 0);
1755 for(auto p: q) q2[p]++;
1756 DEBB(DF_GEOM, ("q2 = ", q2));
1757
1758 bfs = {};
1759 for(int i=0; i<qty; i++)
1760 bfs.emplace_back(i, root[i]);
1761 for(int i=0; i<isize(bfs); i++) {
1762 address last = bfs[i];
1763 int state = last.second;
1764 int fv = last.first;
1765 for(int d=0; d<S7; d++) {
1766 int nstate = children[state*S7+d];
1767 if(nstate < -1) nstate += (1<<16);
1768 if(nstate >= 0) {
1769 address next = {mov(fv, d), nstate};
1770 if(!nles.count(next)) continue;
1771 int c = isize(nles[next]);
1772 nles[next].erase(last);
1773 if(nles[next].empty() && c) {
1774 nles.erase(next);
1775 bfs.push_back(next);
1776 }
1777 }
1778 }
1779 }
1780
1781 DEBB(DF_GEOM, ("removed cases = ", isize(bfs)));
1782
1783 // just the number of FV's
1784 int pstable = 0;
1785 for(auto& p: nonlooping_earlier_states)
1786 pstable = max(pstable, p.first.first+1);
1787
1788 println(hlog, "pstable size = ", pstable, " (states: ", qstate, ")");
1789
1790 possible_states.resize(pstable);
1791 for(auto& p: nonlooping_earlier_states)
1792 possible_states[p.first.first].push_back(p.first.second);
1793 }
1794
hrmap_h3_rulehr::reg3::hrmap_h3_rule1795 hrmap_h3_rule() : fp(0) {
1796
1797 load_ruleset(get_rule_filename());
1798
1799 origin = init_heptagon(S7);
1800 heptagon& h = *origin;
1801 h.s = hsOrigin;
1802 h.fiftyval = root[0];
1803 if(PURE) h.c7 = newCell(S7, origin);
1804
1805 int opos = 0;
1806 for(int c: children) {
1807 if(c < -1) c += (1<<16);
1808 if(c >= 0)
1809 otherpos.push_back(-1);
1810 else {
1811 otherpos.push_back(opos);
1812 while(other[opos] != ',') opos++;
1813 opos++;
1814 }
1815 }
1816
1817 quotient_map = nullptr;
1818
1819 if(geometry == gSpace535)
1820 quotient_map = new seifert_weber::hrmap_seifert_cover();
1821 #if CAP_CRYSTAL
1822 else if(geometry == gSpace344)
1823 quotient_map = new hrmap_from_crystal;
1824 #endif
1825 else
1826 quotient_map = new hrmap_field3(&fp);
1827
1828 if(geometry == gSpace535)
1829 emerald_map = new seifert_weber::hrmap_seifert_cover();
1830 #if CAP_CRYSTAL
1831 else if(geometry == gSpace344)
1832 emerald_map = new hrmap_from_crystal;
1833 #endif
1834 else
1835 emerald_map = new hrmap_field3(&currfp);
1836 h.emeraldval = 0;
1837
1838 find_mappings();
1839
1840 if(!PURE) get_cell_at(origin, 0);
1841 }
1842
getOriginhr::reg3::hrmap_h3_rule1843 heptagon *getOrigin() override {
1844 return origin;
1845 }
1846
1847 #define DEB 0
1848
counterparthr::reg3::hrmap_h3_rule1849 heptagon *counterpart(heptagon *h) {
1850 return quotient_map->allh[h->fieldval];
1851 }
1852
1853 vector<short> evmemo;
1854
find_emeraldvalhr::reg3::hrmap_h3_rule1855 void find_emeraldval(heptagon *target, heptagon *parent, int d) {
1856 if(geometry == gSpace535) {
1857 target->emeraldval = target->fieldval;
1858 target->zebraval = 0;
1859 return;
1860 }
1861 generate_cellrotations();
1862 auto& cr = cgi.cellrotations;
1863 if(evmemo.empty()) {
1864 println(hlog, "starting");
1865 map<int, int> matrix_hashtable;
1866 auto matrix_hash = [] (const transmatrix& M) {
1867 return bucketer(M[0][0])
1868 + bucketer(M[0][1]) * 71
1869 + bucketer(M[0][2]) * 113
1870 + bucketer(M[1][0]) * 1301
1871 + bucketer(M[1][1]) * 1703
1872 + bucketer(M[1][2]) * 17031
1873 + bucketer(M[2][2]) * 2307
1874 + bucketer(M[2][0]) * 2311
1875 + bucketer(M[2][1]) * 10311;
1876 };
1877 for(int i=0; i<isize(cr); i++) matrix_hashtable[matrix_hash(cr[i].M)] = cr[i].inverse_id;
1878 println(hlog, "ids size = ", isize(matrix_hashtable));
1879
1880 for(int eid=0; eid<isize(emerald_map->allh); eid++)
1881 for(int k0=0; k0<isize(cr); k0++)
1882 for(int fv=0; fv<isize(quotient_map->allh); fv++) {
1883 for(int d=0; d<S7; d++) {
1884 int ed = cr[k0].mapping[d];
1885 auto cpart = emerald_map->allh[eid];
1886 int eid1 = emerald_map->allh[eid]->move(ed)->fieldval;
1887 const transmatrix& X = cr[cr[k0].inverse_id].M;
1888 transmatrix U = quotient_map->iadj(quotient_map->allh[fv], d) * X * emerald_map->adj(cpart, ed);
1889 int k1 = matrix_hashtable[matrix_hash(U)];
1890 /* for(int ik1=0; ik1<isize(cr); ik1++) {
1891 auto& mX1 = cr[ik1].M;
1892 if(eqmatrix(mX1, U)) k1 = cr[ik1].inverse_id;
1893 } */
1894 evmemo.push_back(eid1 * isize(cr) + k1);
1895 }
1896 }
1897 println(hlog, "generated ", isize(evmemo));
1898 }
1899 int memo_id = parent->emeraldval;
1900 memo_id = memo_id * isize(quotient_map->allh) + parent->fieldval;
1901 memo_id = memo_id * S7 + d;
1902 target->emeraldval = evmemo[memo_id];
1903 target->zebraval = emerald_map->allh[target->emeraldval / isize(cr)]->zebraval;
1904 }
1905
create_stephr::reg3::hrmap_h3_rule1906 heptagon *create_step(heptagon *parent, int d) override {
1907 int id = parent->fiftyval;
1908 if(id < 0) id += (1<<16);
1909
1910 auto cp = counterpart(parent);
1911 int d2 = cp->c.spin(d);
1912 int fv = cp->c.move(d)->fieldval;
1913
1914 // indenter ind(2);
1915
1916 heptagon *res = nullptr;
1917
1918 int id1 = children[S7*id+d];
1919 int pos = otherpos[S7*id+d];
1920 if(id1 < -1) id1 += (1<<16);
1921
1922 if(id1 == -1 && false) {
1923 int kk = pos;
1924 string s;
1925 while(other[kk] != ',') s += other[kk++];
1926 println(hlog, "id=", id, " d=", d, " d2=", d2, " id1=", id1, " pos=", pos, " s = ", s);
1927 }
1928
1929 if(id1 != -1) {
1930 res = init_heptagon(S7);
1931 if(PURE && parent->c7)
1932 res->c7 = newCell(S7, res);
1933 res->fieldval = fv;
1934 res->distance = parent->distance + 1;
1935 res->fiftyval = id1;
1936 find_emeraldval(res, parent, d);
1937 // res->c.connect(d2, parent, d, false);
1938 }
1939
1940 else if(other[pos] == ('A' + d) && other[pos+1] == ',') {
1941 res = init_heptagon(S7);
1942 res->alt = parent->alt;
1943 res->fieldval = fv;
1944 res->distance = parent->distance - 1;
1945 vector<int> possible;
1946 int pfv = parent->fieldval;
1947 if(geometry == gSpace535) pfv = 0;
1948 for(auto s: nonlooping_earlier_states[address{pfv, id}]) possible.push_back(s.second);
1949 id1 = hrand_elt(possible, 0);
1950 res->fiftyval = id1;
1951 find_emeraldval(res, parent, d);
1952 }
1953
1954 else {
1955 heptagon *at = parent;
1956 while(other[pos] != ',') {
1957 int dir = (other[pos++] & 31) - 1;
1958 // println(hlog, "from ", at, " go dir ", dir);
1959 at = at->cmove(dir);
1960 }
1961 res = at;
1962 }
1963
1964 if(!res) throw hr_exception("res missing");
1965
1966 if(res->move(d2)) println(hlog, "res conflict");
1967
1968 res->c.connect(d2, parent, d, false);
1969 return res;
1970 }
1971
~hrmap_h3_rulehr::reg3::hrmap_h3_rule1972 ~hrmap_h3_rule() {
1973 if(quotient_map) delete quotient_map;
1974 clearfrom(origin);
1975 }
1976
adjhr::reg3::hrmap_h3_rule1977 transmatrix adj(heptagon *h, int d) override {
1978 return quotient_map->adj(h, d);
1979 }
1980
relative_matrixhhr::reg3::hrmap_h3_rule1981 transmatrix relative_matrixh(heptagon *h2, heptagon *h1, const hyperpoint& hint) override {
1982 return relative_matrix_recursive(h2, h1);
1983 }
1984
relative_matrixchr::reg3::hrmap_h3_rule1985 transmatrix relative_matrixc(cell *c2, cell *c1, const hyperpoint& hint) override {
1986 if(PURE) return relative_matrix(c2->master, c1->master, hint);
1987 return relative_matrix_via_masters(c2, c1, hint);
1988 }
1989
master_relativehr::reg3::hrmap_h3_rule1990 transmatrix master_relative(cell *c, bool get_inverse) override {
1991 if(PURE) return Id;
1992 int aid = cell_id.at(c);
1993 return quotient_map->master_relative(quotient_map->acells[aid], get_inverse);
1994 }
1995
shvidhr::reg3::hrmap_h3_rule1996 int shvid(cell *c) override {
1997 if(PURE) return 0;
1998 if(!cell_id.count(c)) return quotient_map->shvid(c);
1999 int aid = cell_id.at(c);
2000 return quotient_map->shvid(quotient_map->acells[aid]);
2001 }
2002
wall_offsethr::reg3::hrmap_h3_rule2003 int wall_offset(cell *c) override {
2004 if(PURE) return 0;
2005 if(!cell_id.count(c)) return quotient_map->wall_offset(c); /* necessary because ray samples are from quotient_map */
2006 int aid = cell_id.at(c);
2007 return quotient_map->wall_offset(quotient_map->acells[aid]);
2008 }
2009
adjhr::reg3::hrmap_h3_rule2010 transmatrix adj(cell *c, int d) override {
2011 if(PURE) return adj(c->master, d);
2012 if(!cell_id.count(c)) return quotient_map->adj(c, d); /* necessary because ray samples are from quotient_map */
2013 int aid = cell_id.at(c);
2014 return quotient_map->tmatrices_cell[aid][d];
2015 }
2016
get_cellshapehr::reg3::hrmap_h3_rule2017 subcellshape& get_cellshape(cell *c) override {
2018 if(PURE) return *cgi.heptshape;
2019 int aid = cell_id.at(c);
2020 return quotient_map->get_cellshape(quotient_map->acells[aid]);
2021 }
2022
2023 map<cell*, int> cell_id;
2024 map<pair<heptagon*, int>, cell*> cell_at;
2025
get_cell_athr::reg3::hrmap_h3_rule2026 cell *get_cell_at(heptagon *h, int acell_id) {
2027 pair<heptagon*, int> p(h, acell_id);
2028 auto& ca = cell_at[p];
2029 if(!ca) {
2030 ca = newCell(quotient_map->acells[acell_id]->type, h);
2031 cell_id[ca] = acell_id;
2032 if(!h->c7) h->c7 = ca;
2033 }
2034 return ca;
2035 }
2036
find_cell_connectionhr::reg3::hrmap_h3_rule2037 void find_cell_connection(cell *c, int d) override {
2038 if(PURE) {
2039 auto h = c->master->cmove(d);
2040 c->c.connect(d, h->c7, c->master->c.spin(d), false);
2041 return;
2042 }
2043 int id = cell_id.at(c);
2044 heptagon *h = c->master;
2045 for(int dir: quotient_map->move_sequences[id][d])
2046 h = h->cmove(dir);
2047 auto ac = quotient_map->acells[id];
2048 cell *c1 = get_cell_at(h, quotient_map->local_id[ac->move(d)].first);
2049 c->c.connect(d, c1, ac->c.spin(d), false);
2050 }
2051
ray_iadjhr::reg3::hrmap_h3_rule2052 transmatrix ray_iadj(cell *c, int i) override {
2053 if(PURE) return iadj(c, i);
2054 if(!cell_id.count(c)) return quotient_map->ray_iadj(c, i); /* necessary because ray samples are from quotient_map */
2055 int aid = cell_id.at(c);
2056 return quotient_map->ray_iadj(quotient_map->acells[aid], i);
2057 }
2058
strafehr::reg3::hrmap_h3_rule2059 cellwalker strafe(cellwalker cw, int j) override {
2060
2061 hyperpoint hfront = tC0(cgi.adjmoves[cw.spin]);
2062 cw.at->cmove(j);
2063 transmatrix T = currentmap->adj(cw.at, j);
2064 cellwalker res1;
2065 for(int i=0; i<S7; i++) if(i != cw.at->c.spin(j))
2066 if(hdist(hfront, T * tC0(cgi.adjmoves[i])) < cgi.strafedist + .01)
2067 res1 = cellwalker(cw.at->cmove(j), i);
2068
2069 int aid = PURE ? cw.at->master->fieldval : cell_id.at(cw.at);
2070 auto res = quotient_map->strafe(cellwalker(quotient_map->acells[aid], cw.spin), j);
2071 cellwalker res2 = cellwalker(cw.at->cmove(j), res.spin);
2072
2073 if(PURE && res1 != res2) println(hlog, "h3: ", res1, " vs ", res2);
2074 return res2;
2075 }
2076
get_move_seqhr::reg3::hrmap_h3_rule2077 const vector<int>& get_move_seq(cell *c, int i) override {
2078 int aid = cell_id.at(c);
2079 return quotient_map->get_move_seq(quotient_map->acells[aid], i);
2080 }
2081
2082 virtual bool link_alt(heptagon *h, heptagon *alt, hstate firststate, int dir) override;
2083 };
2084
2085 struct hrmap_h3_rule_alt : hrmap {
2086
2087 heptagon *origin;
2088
hrmap_h3_rule_althr::reg3::hrmap_h3_rule_alt2089 hrmap_h3_rule_alt(heptagon *o) {
2090 origin = o;
2091 }
2092
2093 };
2094
new_alt_map(heptagon * o)2095 EX hrmap *new_alt_map(heptagon *o) {
2096 return new hrmap_h3_rule_alt(o);
2097 }
2098
link_alt(heptagon * h,heptagon * alt,hstate firststate,int dir)2099 bool hrmap_h3_rule::link_alt(heptagon *h, heptagon *alt, hstate firststate, int dir) {
2100 alt->fieldval = h->fieldval;
2101 if(geometry == gSpace535) alt->fieldval = 0;
2102 if(firststate == hsOrigin) {
2103 alt->fiftyval = root[alt->fieldval];
2104 return true;
2105 }
2106 vector<int>& choices = possible_states[alt->fieldval];
2107 vector<int> choices2;
2108 for(auto c: choices) {
2109 bool ok = true;
2110 for(int d=0; d<S7; d++)
2111 if(h->cmove(d)->distance < h->distance)
2112 if(children[S7*c+d] == -1)
2113 ok = false;
2114 if(ok) choices2.push_back(c);
2115 }
2116 alt->fiftyval = hrand_elt(choices2, -1);
2117 return alt->fiftyval != -1;
2118 }
2119
2120 EX bool reg3_rule_available = true;
2121 EX string other_rule = "";
2122
get_rule_filename()2123 EX string get_rule_filename() {
2124 if(other_rule != "") return other_rule;
2125 switch(geometry) {
2126 case gSpace336: return "honeycomb-rules-336.dat";
2127 case gSpace344: return "honeycomb-rules-344.dat";
2128 // case gSpace345: return "honeycomb-rules-345.dat";
2129 case gSpace353: return "honeycomb-rules-353.dat";
2130 case gSpace354: return "honeycomb-rules-354.dat";
2131 // case gSpace355: return "honeycomb-rules-355.dat";
2132 case gSpace435: return "honeycomb-rules-435.dat";
2133 case gSpace436: return "honeycomb-rules-436.dat";
2134 case gSpace534: return "honeycomb-rules-534.dat";
2135 case gSpace535: return "honeycomb-rules-535.dat";
2136 case gSpace536: return "honeycomb-rules-536.dat";
2137
2138 default: return "";
2139 }
2140 }
2141
in_rule()2142 EX bool in_rule() {
2143 return reg3_rule_available && get_rule_filename() != "";
2144 }
2145
rule_get_root(int i)2146 EX int rule_get_root(int i) {
2147 return ((hrmap_h3_rule*)currentmap)->root[i];
2148 }
2149
rule_get_children()2150 EX const vector<short>& rule_get_children() {
2151 return ((hrmap_h3_rule*)currentmap)->children;
2152 }
2153
new_map()2154 EX hrmap* new_map() {
2155 if(geometry == gSeifertCover) return new seifert_weber::hrmap_seifert_cover;
2156 if(geometry == gSeifertWeber) return new seifert_weber::hrmap_singlecell(108*degree);
2157 if(geometry == gHomologySphere) return new seifert_weber::hrmap_singlecell(36*degree);
2158 if(quotient && !sphere) return new hrmap_field3(&currfp);
2159 if(in_rule()) return new hrmap_h3_rule;
2160 if(sphere) return new hrmap_sphere3;
2161 return new hrmap_h3;
2162 }
2163
hypmap()2164 hrmap_h3* hypmap() {
2165 return ((hrmap_h3*) currentmap);
2166 }
2167
quotient_count()2168 EX int quotient_count() {
2169 return isize(hypmap()->quotient_map->allh);
2170 }
2171
2172 /** This is a generalization of hyperbolic_celldistance in expansion.cpp to three dimensions.
2173 It still assumes that there are at most 4 cells around every edge, and that distances from
2174 the origin are known, so it works only in {5,3,4}.
2175 */
2176
celldistance_534(cell * c1,cell * c2)2177 int celldistance_534(cell *c1, cell *c2) {
2178 int d1 = celldist(c1);
2179 int d2 = celldist(c2);
2180
2181 vector<cell*> s1 = {c1};
2182 vector<cell*> s2 = {c2};
2183 int best = 99999999;
2184 int d0 = 0;
2185
2186 auto go_nearer = [&] (vector<cell*>& v, int& d) {
2187 vector<cell*> w;
2188 for(cell *c: v)
2189 forCellEx(c1, c)
2190 if(celldist(c1) < d)
2191 w.push_back(c1);
2192 sort(w.begin(), w.end());
2193 d--; d0++;
2194 auto last = std::unique(w.begin(), w.end());
2195 w.erase(last, w.end());
2196 v = w;
2197 };
2198
2199 while(d0 < best) {
2200 for(cell *a1: s1) for(cell *a2: s2) {
2201 if(a1 == a2) best = min(best, d0);
2202 else if(isNeighbor(a1, a2)) best = min(best, d0+1);
2203 }
2204
2205 if(d1 == 0 && d2 == 0) break;
2206
2207 if(d1 >= d2) go_nearer(s1, d1);
2208 if(d1 < d2) go_nearer(s2, d2);
2209 }
2210
2211 return best;
2212 }
2213
2214
celldistance(cell * c1,cell * c2)2215 EX int celldistance(cell *c1, cell *c2) {
2216 if(c1 == c2) return 0;
2217 if(c1 == currentmap->gamestart()) return c2->master->distance;
2218 if(c2 == currentmap->gamestart()) return c1->master->distance;
2219
2220 if(geometry == gSpace534 && PURE) return celldistance_534(c1, c2);
2221
2222 auto r = hypmap();
2223
2224 hyperpoint h = tC0(r->relative_matrix(c1->master, c2->master, C0));
2225 int b = bucketer(h);
2226 if(cgi.close_distances.count(b)) return cgi.close_distances[b];
2227
2228 if(in_rule())
2229 return clueless_celldistance(c1, c2);
2230
2231 dynamicval<eGeometry> g(geometry, gBinary3);
2232 #if CAP_BT
2233 return 20 + bt::celldistance3(r->reg_gmatrix[c1->master].first, r->reg_gmatrix[c2->master].first);
2234 #else
2235 return 20;
2236 #endif
2237 }
2238
pseudohept(cell * c)2239 EX bool pseudohept(cell *c) {
2240 if(sphere) {
2241 auto m = currentmap;
2242 hyperpoint h = tC0(m->relative_matrix(c->master, m->getOrigin(), C0));
2243 if(S7 == 12) {
2244 hyperpoint h1 = cspin(0, 1, atan2(16, 69) + M_PI/4) * h;
2245 for(int i=0; i<4; i++) if(abs(abs(h1[i]) - .5) > .01) return false;
2246 return true;
2247 }
2248 if(S7 == 8)
2249 return h[3] >= .99 || h[3] <= -.99 || abs(h[3]) < .01;
2250 if(cgi.loop == 3 && cgi.face == 3 && S7 == 4)
2251 return c == m->gamestart();
2252 if(cgi.loop == 4 && cgi.face == 3)
2253 return abs(h[3]) > .9;
2254 if(cgi.loop == 3 && cgi.face == 4)
2255 return abs(h[3]) > .9;
2256 if(cgi.loop == 5 && cgi.face == 3)
2257 return abs(h[3]) > .99 || abs(h[0]) > .99 || abs(h[1]) > .99 || abs(h[2]) > .99;
2258 }
2259 auto m = hypmap();
2260 if(cgflags & qSINGLE) return true;
2261 if(fake::in()) return FPIU(reg3::pseudohept(c));
2262 // chessboard pattern in 534
2263 if(geometry == gField534)
2264 return hr::celldistance(c, currentmap->gamestart()) & 1;
2265 if(geometry == gCrystal344 || geometry == gCrystal534 || geometry == gSeifertCover)
2266 return false;
2267 if(quotient) return false; /* added */
2268 auto mr = dynamic_cast<hrmap_h3_rule*> (currentmap);
2269 if(mr) {
2270 if(geometry == gSpace535)
2271 return c->master->fieldval % 31 == 0;
2272 return c->master->fieldval == 0;
2273 }
2274 if(m && hyperbolic) {
2275 heptagon *h = m->reg_gmatrix[c->master].first;
2276 return (h->zebraval == 1) && (h->distance & 1);
2277 }
2278
2279 return false;
2280 }
2281
generate_cellrotations()2282 EX void generate_cellrotations() {
2283 auto &cr = cgi.cellrotations;
2284 if(isize(cr)) return;
2285
2286 for(int a=0; a<S7; a++)
2287 for(int b=0; b<S7; b++)
2288 for(int c=0; c<S7; c++) {
2289 transmatrix T = build_matrix(cgi.adjmoves[a]*C0, cgi.adjmoves[b]*C0, cgi.adjmoves[c]*C0, C0);
2290 if(abs(det(T)) < 0.001) continue;
2291 transmatrix U = build_matrix(cgi.adjmoves[0]*C0, cgi.adjmoves[1]*C0, cgi.adjmoves[2]*C0, C0);
2292 transmatrix S = U * inverse(T);
2293 if(abs(det(S) - 1) > 0.01) continue;
2294 vector<int> perm(S7);
2295 for(int x=0; x<S7; x++) perm[x] = -1;
2296 for(int x=0; x<S7; x++)
2297 for(int y=0; y<S7; y++)
2298 if(hdist(S * cgi.adjmoves[x] * C0, cgi.adjmoves[y] * C0) < .1) perm[x] = y;
2299 bool bad = false;
2300 for(int x=0; x<S7; x++) if(perm[x] == -1) bad = true;
2301 if(bad) continue;
2302
2303 cr.emplace_back(geometry_information::cellrotation_t{S, perm, 0});
2304 }
2305
2306 int rots = isize(cr);
2307 for(int i=0; i<rots; i++)
2308 for(int j=0; j<rots; j++)
2309 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)
2310 cr[i].inverse_id = j;
2311 }
2312 #endif
2313
2314 #if 0
2315 /* More precise, but very slow distance. Not used/optimized for now */
2316
2317 ld adistance(cell *c) {
2318 hyperpoint h = tC0(regmap()->reg_gmatrix[c->master].second);
2319 h = bt::deparabolic3(h);
2320 return regmap()->reg_gmatrix[c->master].first->distance * log(2) - h[0];
2321 }
2322
2323 map<pair<cell*, cell*>, int> memo;
2324
2325 bool cdd;
2326
2327 int celldistance(cell *c1, cell *c2) {
2328 if(memo.count(make_pair(c1, c2))) return memo[make_pair(c1, c2)];
2329 if(c1 == c2) return 0;
2330 vector<cell*> v[2];
2331 v[0].push_back(c1);
2332 v[1].push_back(c2);
2333
2334 int steps = 0;
2335
2336 map<cell*, int> visited;
2337 visited[c1] = 1;
2338 visited[c2] = 2;
2339
2340 while(true) {
2341 if(cdd) {
2342 println(hlog, "state ", steps, "/",isize(v[0]), "/", isize(v[1]));
2343 println(hlog, " A: ", v[0]);
2344 println(hlog, " B: ", v[1]);
2345 }
2346 for(int i: {0,1}) {
2347 vector<cell*> new_v;
2348 for(cell *c: v[i]) forCellCM(cn, c) if(adistance(cn) < adistance(c)) {
2349 auto &vi = visited[cn];
2350 if((vi&3) == 0) {
2351 vi = 4 * (steps+1);
2352 vi |= (1<<i);
2353 new_v.push_back(cn);
2354 }
2355 else if((vi&3) == 2-i) {
2356 vector<pair<cell*, int>> ca1, ca2;
2357 int b1 = 4*steps-4;
2358 int b2 = ((vi>>2)<<2) - 4;
2359 for(auto p: visited) {
2360 if(cdd) println(hlog, p);
2361 int ps = p.second & 3;
2362 if(ps == 1+i && p.second >= b1)
2363 ca1.emplace_back(p.first, p.second/4);
2364 if(ps == 2-i && p.second >= b2 && p.second <= b2+8)
2365 ca2.emplace_back(p.first, p.second/4);
2366 }
2367 int bound = 1<<16;
2368 for(auto p1: ca1) for(auto p2: ca2) {
2369 hyperpoint h = tC0(relative_matrix(p1.first->master, p2.first->master));
2370 int b = bucketer(h);
2371 if(close_distances.count(b)) {
2372 int d = close_distances[b] + p1.second + p2.second;
2373 if(cdd) println(hlog, "candidate: close=", close_distances[b], p1, p2, "; h = ", h);
2374 if(d < bound) bound = d;
2375 }
2376 else if(cdd) println(hlog, "bucket missing");
2377 }
2378 return memo[make_pair(c1, c2)] = bound;
2379 return bound;
2380 }
2381 }
2382 v[i] = std::move(new_v);
2383 }
2384 steps++;
2385 }
2386 }
2387
2388 cellwalker target;
2389 int tsteps;
2390
2391 int dist_alt(cell *c) {
2392 if(!target.at) {
2393 target = cellwalker(currentmap->gamestart(), 0);
2394 tsteps = 0;
2395 for(int i=0; i<30; i++) target += wstep, target += rev, tsteps++;
2396 }
2397 if(specialland == laCamelot) return reg3::celldistance(c, target.at);
2398 else {
2399 int d = reg3::celldistance(c, target.at) - tsteps;
2400 if(d < 10) target += wstep, target += rev, tsteps++;
2401 return d;
2402 }
2403 }
2404 #endif
2405
2406 // Construct a cellwalker in direction j from cw.at, such that its direction is as close
2407 // as possible to cw.spin. Assume that j and cw.spin are adjacent
2408
2409 #if MAXMDIM >= 4
matrix_order(const transmatrix A)2410 EX int matrix_order(const transmatrix A) {
2411 transmatrix T = A;
2412 int res = 1;
2413 while(!eqmatrix(T, Id)) {
2414 res++; T = T * A;
2415 }
2416 return res;
2417 }
2418
generate_fulls()2419 EX void generate_fulls() {
2420 reg3::generate_cellrotations();
2421
2422 auto cons = [&] (int i0, int i1, int i2) {
2423 transmatrix T = build_matrix(cgi.adjmoves[ 0]*C0, cgi.adjmoves[ 1]*C0, cgi.adjmoves[ 2]*C0, C0);
2424 transmatrix U = build_matrix(cgi.adjmoves[i0]*C0, cgi.adjmoves[i1]*C0, cgi.adjmoves[i2]*C0, C0);
2425 return U * inverse(T);
2426 };
2427
2428 cgi.full_P = cgi.adjmoves[0];
2429 cgi.full_R = S7 == 8 ? cons(1, 7, 0) : S7 == 20 ? cons(1,2,6) : cons(1, 2, 0);
2430 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);
2431
2432 cgi.xp_order = matrix_order(cgi.full_X * cgi.full_P);
2433 cgi.r_order = matrix_order(cgi.full_R);
2434 cgi.rx_order = matrix_order(cgi.full_R * cgi.full_X);
2435 println(hlog, "orders = ", tie(cgi.rx_order, cgi.r_order, cgi.xp_order));
2436 }
2437
construct_relations()2438 EX void construct_relations() {
2439 auto& rels = cgi.rels;
2440 if(!rels.empty()) return;
2441 rels.clear();
2442
2443 reg3::generate_cellrotations();
2444 reg3::generate_fulls();
2445 vector<transmatrix> all;
2446
2447 vector<string> formulas;
2448
2449 formulas.push_back("");
2450
2451 all.push_back(Id);
2452 auto& faces = cgi.heptshape->faces;
2453 hyperpoint v = faces[0][0];
2454 auto add = [&] (transmatrix T) {
2455 for(int i=0; i<isize(all); i++) if(eqmatrix(all[i], T)) return i;
2456 int S = isize(all);
2457 all.push_back(T);
2458 return S;
2459 };
2460
2461 println(hlog, faces);
2462
2463 println(hlog, "cellshape = ", isize(faces));
2464 bool ok = true;
2465 int last_i = -1;
2466 for(auto& v: faces) for(hyperpoint h: v) {
2467 int i = 0, j = 0;
2468 for(auto& uv: faces) for(hyperpoint u: uv) {
2469 if(hdist(h, cgi.full_X*u) < 5e-2) i++;
2470 if(hdist(h, cgi.full_R*u) < 5e-2) j++;
2471 }
2472 if(last_i == -1) last_i = i;
2473 if(i != j || i != last_i) ok = false;
2474 }
2475
2476 if(!ok) { println(hlog, "something wrong"); exit(1); }
2477
2478 add(Id);
2479
2480 auto work = [&] (transmatrix T, int p, char c) {
2481 if(hdist0(tC0(T)) > 5) return;
2482 for(auto& hv: faces) for(hyperpoint h: hv) if(hdist(T * h, v) < 1e-4) goto ok;
2483 return;
2484 ok:
2485 int id = add(T);
2486 // println(hlog, p, " x ", (s0+c), " = ", id);
2487
2488 if(id >= isize(formulas)) formulas.push_back(formulas[p] + c);
2489 else if(id == 0) println(hlog, "reached identity: ", formulas[p]+c);
2490 else if(formulas[p][0] != formulas[id][0])
2491 rels.emplace_back(formulas[p] + c, formulas[id]);
2492 };
2493
2494 for(int i=0; i<isize(all); i++) {
2495 transmatrix T = all[i];
2496 work(T * cgi.full_R, i, 'R');
2497 work(T * cgi.full_X, i, 'X');
2498 work(T * cgi.full_P, i, 'P');
2499 }
2500 }
2501
2502 eVariation target_variation;
2503 flagtype target_coxeter;
2504 int target_subcube_count;
2505
edit_variation()2506 EX void edit_variation() {
2507 cmode = sm::SIDE | sm::MAYDARK;
2508 gamescreen(0);
2509 dialog::init(XLAT("variations"));
2510
2511 dialog::addBoolItem(XLAT("pure"), target_variation == eVariation::pure, 'p');
2512 dialog::add_action([] { target_variation = eVariation::pure; });
2513
2514 dialog::addBoolItem(XLAT("symmetric subdivision"), target_variation == eVariation::coxeter, 't');
2515 dialog::add_action([] { target_variation = eVariation::coxeter; });
2516
2517 if(S7 == 6) {
2518 dialog::addBoolItem(XLAT("sub-cubes"), target_variation == eVariation::subcubes, 'c');
2519 dialog::add_action([] { target_variation = eVariation::subcubes; });
2520
2521 if(!(cgflags & qIDEAL)) {
2522 dialog::addBoolItem(XLAT("dual sub-cubes"), target_variation == eVariation::dual_subcubes, 'd');
2523 dialog::add_action([] { target_variation = eVariation::dual_subcubes; });
2524
2525 dialog::addBoolItem(XLAT("bitruncated sub-cubes"), target_variation == eVariation::bch, 'b');
2526 dialog::add_action([] { target_variation = eVariation::bch; });
2527 }
2528 }
2529
2530 else
2531 dialog::addInfo(XLAT("note: more choices in cubic honeycombs"));
2532
2533 if(is_subcube_based(target_variation)) {
2534 dialog::addBreak(100);
2535 dialog::addSelItem(XLAT("subdivision"), its(target_subcube_count), 'z');
2536 dialog::add_action([] {
2537 dialog::editNumber(target_subcube_count, 1, 8, 1, 2, XLAT("subdivision"), "");
2538 dialog::bound_low(1);
2539 });
2540 }
2541
2542 if(target_variation == eVariation::coxeter) {
2543 dialog::addBreak(100);
2544 dialog::addBoolItem(XLAT("split by original faces"), target_coxeter & cox_othercell, 'f');
2545 dialog::add_action([] { target_coxeter ^= cox_othercell; });
2546 dialog::addBoolItem(XLAT("split by vertex axes"), target_coxeter & cox_vertices, 'v');
2547 dialog::add_action([] { target_coxeter ^= cox_vertices; });
2548 dialog::addBoolItem(XLAT("split by midedges"), target_coxeter & cox_midedges, 'm');
2549 dialog::add_action([] { target_coxeter ^= cox_midedges; });
2550 }
2551
2552 dialog::addBreak(100);
2553 dialog::addItem(XLAT("activate"), 'x');
2554 dialog::add_action([] {
2555 stop_game();
2556 set_variation(target_variation);
2557 subcube_count = target_subcube_count;
2558 coxeter_param = target_coxeter;
2559 start_game();
2560 });
2561 dialog::addBack();
2562 dialog::display();
2563 }
2564
configure_variation()2565 EX void configure_variation() {
2566 target_variation = variation;
2567 target_subcube_count = subcube_count;
2568 target_coxeter = coxeter_param;
2569 pushScreen(edit_variation);
2570 }
2571
2572 EX }
2573 #endif
2574
2575 #if MAXMDIM == 3
2576 EX namespace reg3 {
in()2577 EX bool in() { return false; }
in_rule()2578 EX bool in_rule() { return false; }
2579 EX }
2580 #endif
2581
2582 }
2583
2584