1 // Hyperbolic Rogue -- Goldberg-Coxeter construction
2 // Copyright (C) 2011-2019 Zeno Rogue, see 'hyper.cpp' for details
3 
4 /** \file goldberg.cpp
5  *  \brief Goldberg-Coxeter construction
6  *
7  *  This is generally not used for standard pure and bitruncated tilings, even though they are technically Goldberg too.
8  */
9 
10 #include "hyper.h"
11 namespace hr {
12 
13 #if HDR
14 struct hrmap;
15 extern hrmap *currentmap;
16 #endif
17 
18 EX namespace gp {
19 
20   #if HDR
21   struct loc : pair<int, int> {
lochr::gp::loc22     loc() {}
23 
lochr::gp::loc24     loc(int x, int y) : pair<int,int> (x,y) {}
25 
operator +hr::gp::loc26     loc operator+(loc e2) {
27       return loc(first+e2.first, second+e2.second);
28       }
29 
operator -hr::gp::loc30     loc operator-(loc e2) {
31       return loc(first-e2.first, second-e2.second);
32       }
33 
operator *hr::gp::loc34     loc operator*(loc e2) {
35       return loc(first*e2.first-second*e2.second,
36         first*e2.second + e2.first*second + (S3 == 3 ? second*e2.second : 0));
37       }
38 
operator *hr::gp::loc39     loc operator*(int i) {
40       return loc(first*i, second*i);
41       }
42 
operator %hr::gp::loc43     int operator %(int i) {
44       return gmod(first, i) + gmod(second, i);
45       }
46 
operator /hr::gp::loc47     loc operator /(int i) {
48       return loc(first/i, second/i);
49       }
50 
conjhr::gp::loc51     loc conj() {
52       if(S3 == 4) return loc(first, -second);
53       return loc(first+second, -second);
54       }
55 
56     };
57 
58   struct local_info {
59     int last_dir;
60     loc relative;
61     int first_dir;
62     int total_dir;
63     };
64   #endif
65 
66   EX local_info current_li;
67   EX cell *li_for;
68 
eudir(int d)69   EX loc eudir(int d) {
70     if(S3 == 3) {
71       d %= 6; if (d < 0) d += 6;
72       switch(d) {
73         case 0: return loc(1, 0);
74         case 1: return loc(0, 1);
75         case 2: return loc(-1, 1);
76         case 3: return loc(-1, 0);
77         case 4: return loc(0, -1);
78         case 5: return loc(1, -1);
79         default: return loc(0, 0);
80         }
81       }
82     else switch(d&3) {
83       case 0: return loc(1, 0);
84       case 1: return loc(0, 1);
85       case 2: return loc(-1, 0);
86       case 3: return loc(0, -1);
87       default: return loc(0, 0);
88       }
89     }
90 
length(loc p)91   EX int length(loc p) {
92     return euc::dist(p.first, p.second);
93     }
94 
95 #if CAP_GP
96   EX loc param = loc(1, 0);
97 
98   EX hyperpoint next;
99 
100   struct goldberg_mapping_t {
101     cellwalker cw;
102     signed char rdir;
103     signed char mindir;
104     loc start;
105     transmatrix adjm;
106     };
107 
fixg6(int x)108   EX int fixg6(int x) { return gmod(x, SG6); }
109 
110   const int GOLDBERG_LIMIT_HALF = GOLDBERG_LIMIT/2;
111   const int GOLDBERG_MASK_HALF = GOLDBERG_MASK/2;
112 
get_code(const local_info & li)113   EX int get_code(const local_info& li) {
114     return
115       ((li.relative.first & GOLDBERG_MASK_HALF) << 0) +
116       ((li.relative.second & GOLDBERG_MASK_HALF) << (GOLDBERG_BITS-1)) +
117       ((fixg6(li.total_dir)) << (2*GOLDBERG_BITS-2)) +
118       ((li.last_dir & 15) << (2*GOLDBERG_BITS+2));
119     }
120 
get_local_info(cell * c)121   EX local_info get_local_info(cell *c) {
122     if(INVERSE) {
123       c = get_mapped(c);
124       return UIU(get_local_info(c));
125       }
126     local_info li;
127     if(c == c->master->c7) {
128       li.relative = loc(0,0);
129       li.first_dir = -1;
130       li.last_dir = -1;
131       li.total_dir = -1;
132       }
133     else {
134       vector<int> dirs;
135       while(c != c->master->c7) {
136         dirs.push_back(c->c.spin(0));
137         c = c->move(0);
138         }
139       li.first_dir = dirs[0];
140       li.last_dir = dirs.back();
141 
142       loc at(0,0);
143       int dir = 0;
144       at = at + eudir(dir);
145       dirs.pop_back();
146       while(dirs.size()) {
147         dir += dirs.back() + SG3;
148         dirs.pop_back();
149         at = at + eudir(dir);
150         }
151       li.relative = at;
152       li.total_dir = dir + SG3;
153       }
154     return li;
155     }
156 
last_dir(cell * c)157   EX int last_dir(cell *c) {
158     return get_local_info(c).last_dir;
159     }
160 
get_coord(cell * c)161   EX loc get_coord(cell *c) {
162     return get_local_info(c).relative;
163     }
164 
pseudohept_val(cell * c)165   EX int pseudohept_val(cell *c) {
166     loc v = get_coord(c);
167     return gmod(v.first - v.second, 3);
168     }
169 
170   // mapping of the local equilateral triangle
171   // goldberg_map[y][x].cw is the cellwalker in this triangle at position (x,y)
172   // facing local direction 0
173 
174   goldberg_mapping_t goldberg_map[GOLDBERG_LIMIT][GOLDBERG_LIMIT];
clear_mapping()175   void clear_mapping() {
176     for(int y=0; y<GOLDBERG_LIMIT; y++) for(int x=0; x<GOLDBERG_LIMIT; x++) {
177       goldberg_map[y][x].cw.at = NULL;
178       goldberg_map[y][x].rdir = -1;
179       goldberg_map[y][x].mindir = 0;
180       }
181     }
182 
get_mapping(loc c)183   goldberg_mapping_t& get_mapping(loc c) {
184     return goldberg_map[c.second&GOLDBERG_MASK][c.first&GOLDBERG_MASK];
185     }
186 
187   int spawn;
188 
peek(cellwalker cw)189   cell*& peek(cellwalker cw) {
190     return cw.at->move(cw.spin);
191     }
192 
get_localwalk(const goldberg_mapping_t & wc,int dir)193   cellwalker get_localwalk(const goldberg_mapping_t& wc, int dir) {
194     if(dir < wc.mindir) dir += SG6;
195     if(dir >= wc.mindir + SG6) dir -= SG6;
196     return wc.cw + dir;
197     }
198 
set_localwalk(goldberg_mapping_t & wc,int dir,const cellwalker & cw)199   void set_localwalk(goldberg_mapping_t& wc, int dir, const cellwalker& cw) {
200     if(dir < wc.mindir) dir += SG6;
201     if(dir >= wc.mindir + SG6) dir -= SG6;
202     wc.cw = cw - dir;
203     }
204 
pull(loc at,int dir)205   bool pull(loc at, int dir) {
206     auto& wc = get_mapping(at);
207     auto at1 = at + eudir(dir);
208     int dir1 = fixg6(dir+SG3);
209     cellwalker wcw = get_localwalk(wc, dir);
210     auto& wc1= get_mapping(at1);
211     if(wc1.cw.at) {
212       if(peek(wcw)) {
213         auto wcw1 = get_localwalk(wc1, dir1);
214         if(wcw + wstep != wcw1) {
215           DEBB(DF_GP, (at1, " : ", (wcw+wstep), " / ", wcw1, " (pull error from ", at, " :: ", wcw, ")") );
216           exit(1);
217           }
218         if(do_adjm) wc1.adjm = wc.adjm * get_adj(wcw.at, wcw.spin);
219         }
220       return false;
221       }
222     if(peek(wcw)) {
223       set_localwalk(wc1, dir1, wcw + wstep);
224       DEBB(DF_GP, (at1, " :", wcw+wstep, " (pulled from ", at, " :: ", wcw, ")"));
225       if(do_adjm) wc1.adjm = wc.adjm * get_adj(wcw.at, wcw.spin);
226       return true;
227       }
228     return false;
229     }
230 
231   EX bool do_adjm;
232 
conn1(loc at,int dir,int dir1)233   void conn1(loc at, int dir, int dir1) {
234     auto& wc = get_mapping(at);
235     auto wcw = get_localwalk(wc, dir);
236     auto& wc1 = get_mapping(at + eudir(dir));
237     DEBB0(DF_GP, (format("  md:%02d s:%d", wc.mindir, wc.cw.spin)); )
238     DEBB0(DF_GP, ("  connection ", at, "/", dir, " ", wc.cw+dir, "=", wcw, " ~ ", at+eudir(dir), "/", dir1); )
239     if(!wc1.cw.at) {
240       wc1.start = wc.start;
241       if(peek(wcw)) {
242         DEBB0(DF_GP, ("(pulled) "); )
243         set_localwalk(wc1, dir1, wcw + wstep);
244         if(do_adjm) wc1.adjm = wc.adjm * get_adj(wcw.at, wcw.spin);
245         }
246       else {
247         peek(wcw) = newCell(SG6, wc.cw.at->master);
248         wcw.at->c.setspin(wcw.spin, 0, false);
249         set_localwalk(wc1, dir1, wcw + wstep);
250         if(do_adjm) wc1.adjm = wc.adjm;
251         spawn++;
252         DEBB0(DF_GP, ("(created) "); )
253         }
254       }
255     DEBB0(DF_GP, (wc1.cw+dir1, " "));
256     auto wcw1 = get_localwalk(wc1, dir1);
257     if(peek(wcw)) {
258       if(wcw+wstep != wcw1) {
259         DEBB(DF_GP, ("FAIL: ", wcw, " / ", wcw1); exit(1); )
260         }
261       else {
262         DEBB(DF_GP, ("(was there)"));
263         }
264       }
265     else {
266       DEBB(DF_GP, ("ok"));
267       peek(wcw) = wcw1.at;
268       wcw.at->c.setspin(wcw.spin, wcw1.spin, wcw.mirrored != wcw1.mirrored);
269       if(wcw+wstep != wcw1) {
270         DEBB(DF_GP | DF_ERROR, ("assertion failed"));
271         exit(1);
272         }
273       }
274     if(do_adjm) {
275       get_adj(wcw.at, wcw.spin) = inverse(wc.adjm) * wc1.adjm;
276       get_adj(wcw1.at, wcw1.spin) = inverse(wc1.adjm) * wc.adjm;
277       }
278     }
279 
conn(loc at,int dir)280   void conn(loc at, int dir) {
281     conn1(at, fixg6(dir), fixg6(dir+SG3));
282     conn1(at + eudir(dir), fixg6(dir+SG3), fixg6(dir));
283     }
284 
285   EX map<pair<cell*, int>, transmatrix> gp_adj;
286 
get_adj(cell * c,int i)287   EX transmatrix& get_adj(cell *c, int i) { return gp_adj[make_pair(c,i)]; }
288 
set_heptspin(loc at,heptspin hs)289   goldberg_mapping_t& set_heptspin(loc at, heptspin hs) {
290     auto& ac0 = get_mapping(at);
291     ac0.cw = cellwalker(hs.at->c7, hs.spin, hs.mirrored);
292     ac0.start = at;
293     DEBB(DF_GP, (at, " : ", ac0.cw));
294     return ac0;
295     }
296 
extend_map(cell * c,int d)297   EX void extend_map(cell *c, int d) {
298     DEBB(DF_GP, ("EXTEND ",c, " ", d));
299     if(c->master->c7 != c) {
300       while(c->master->c7 != c) {
301         DEBB(DF_GP, (c, " direction 0 corresponds to ", c->move(0), " direction ", c->c.spin(0)); )
302         d = c->c.spin(0);
303         c = c->move(0);
304         }
305       // c move 0 equals c' move spin(0)
306       extend_map(c, d);
307       extend_map(c, c->c.fix(d-1));
308       extend_map(c, c->c.fix(d+1));
309       if(S3 == 4 && !c->move(d))
310         for(int i=0; i<S7; i++)
311         for(int j=0; j<S7; j++)
312           extend_map(createStep(c->master, i)->c7, j);
313       return;
314       }
315 
316     if(S3 == 4 && param.first <= param.second) { d--; if(d<0) d += S7; }
317     clear_mapping();
318 
319     // we generate a local map from an Euclidean grid to the
320     // hyperbolic grid we build.
321 
322     // we fill the equilateral triangle with the following vertices:
323 
324     loc vc[4];
325     vc[0] = loc(0,0);
326     vc[1] = param;
327     if(S3 == 3)
328       vc[2] = param * loc(0,1);
329     else
330       vc[2] = param * loc(1,1),
331       vc[3] = param * loc(0,1);
332 
333     heptspin hs(c->master, d, false);
334 
335     auto& ac0 = set_heptspin(vc[0], hs);
336     ac0.mindir = -1;
337     auto& ac1 = set_heptspin(vc[1], hs + wstep - SG3);
338     ac1.mindir = 0;
339     auto& ac2 = set_heptspin(vc[S3-1], S3 == 3 ? hs + 1 + wstep - 4 : hs + 1 + wstep + 1);
340     ac2.mindir = S3 == 3 ? 1 : -2;
341     if(S3 == 4) {
342       set_heptspin(vc[2], hs + wstep - 1 + wstep + 1).mindir = -3;
343       }
344 
345     do_adjm = quotient || sphere;
346     if(do_adjm) {
347       auto m = (hrmap_standard*)currentmap;
348       get_mapping(vc[0]).adjm = Id;
349       get_mapping(vc[1]).adjm = m->adj(c->master, d);
350       get_mapping(vc[S3-1]).adjm = m->adj(c->master, (d+1)%c->master->type);
351       if(S3 == 4) {
352         heptspin hs1 = hs + wstep - 1;
353         get_mapping(vc[2]).adjm = m->adj(c->master, d) * m->adj(hs1.at, hs1.spin);
354         }
355       }
356 
357     if(S3 == 4 && param == loc(1,1)) {
358       conn(loc(0,0), 1);
359       conn(loc(0,1), 0);
360       conn(loc(0,1), 1);
361       conn(loc(0,1), 2);
362       conn(loc(0,1), 3);
363       return;
364       }
365 
366     if(nonorientable && param.first == param.second) {
367       int x = param.first;
368       if(ac1.cw.mirrored != hs.mirrored) ac1.cw--;
369       if(ac2.cw.mirrored != hs.mirrored) ac2.cw--;
370 
371       for(int d=0; d<3; d++) for(int k=0; k<3; k++)
372       for(int i=0; i<x; i++) {
373         int dd = (2*d+k);
374         loc cx = vc[d] + eudir(dd) * i;
375         if(!pull(cx, dd)) break;
376         }
377 
378       for(int i=0; i<=2*x; i++)
379       for(int d=0; d<3; d++) {
380         loc cx = vc[d] + eudir(1+2*d) * i;
381         if(i < 2*x) conn(cx, 1+2*d);
382         int jmax = x-i, drev = 2*d;
383         if(jmax < 0) drev += 3, jmax = -jmax;
384         for(int j=0; j<jmax; j++) {
385           loc cy = cx + eudir(drev) * j;
386           conn(cy, drev);
387           conn(cy, drev+1);
388           conn(cy, drev+2);
389           }
390         }
391 
392       return;
393       }
394 
395     // then we set the edges of our big equilateral triangle (in a symmetric way)
396     for(int i=0; i<S3; i++) {
397       loc start = vc[i];
398       loc end = vc[(i+1)%S3];
399       DEBB(DF_GP, ("from ", start, " to ", end); )
400       loc rel = param;
401       auto build = [&] (loc& at, int dx, bool forward) {
402         int dx1 = dx + SG2*i;
403         DEBB(DF_GP, (at, " .. ", make_pair(at + eudir(dx1), fixg6(dx1+SG3))));
404         conn(at, dx1);
405         if(forward) get_mapping(at).rdir = fixg6(dx1);
406         else get_mapping(at+eudir(dx1)).rdir = fixg6(dx1+SG3);
407         at = at + eudir(dx1);
408         };
409       while(rel.first >= 2 && (S3 == 3 ? rel.first >= 2 - rel.second : true)) {
410         build(start, 0, true);
411         build(end, SG3, false);
412         rel.first -= 2;
413         }
414       while(rel.second >= 2) {
415         build(start, 1, true);
416         build(end, 1+SG3, false);
417         rel.second -= 2;
418         }
419       while(rel.second <= -2 && S3 == 3) {
420         build(start, 5, true);
421         build(end, 2, false);
422         rel.second += 2;
423         rel.first -= 2;
424         }
425       if(S3 == 3) while((rel.first>0 && rel.second > 0) | (rel.first > 1 && rel.second < 0)) {
426         build(start, 0, true);
427         build(end, 3, false);
428         rel.first -= 2;
429         }
430       if(S3 == 4 && rel == loc(1,1)) {
431         if(param == loc(3,1) || param == loc(5,1)) {
432           build(start, 1, true);
433           build(end, 2, false);
434           rel.first--;
435           rel.second--;
436           }
437         else {
438           build(start, 0, true);
439           build(end, 3, false);
440           rel.first--;
441           rel.second--;
442           }
443         }
444       for(int k=0; k<SG6; k++)
445         if(start + eudir(k+SG2*i) == end)
446           build(start, k, true);
447       if(start != end) { DEBB(DF_GP | DF_ERROR, ("assertion failed: start ", start, " == end ", end)); exit(1); }
448       }
449 
450     // now we can fill the interior of our big equilateral triangle
451     loc at = vc[0];
452     int maxstep = 3000;
453     while(true) {
454       maxstep--; if(maxstep < 0) { DEBB(DF_GP | DF_ERROR, ("maxstep exceeded")); exit(1); }
455       auto& wc = get_mapping(at);
456       int dx = wc.rdir;
457       auto at1 = at + eudir(dx);
458       auto& wc1 = get_mapping(at1);
459       DEBB(DF_GP, (make_pair(at, dx), " ", make_pair(at1, wc1.rdir)));
460       int df = wc1.rdir - dx;
461       if(df < 0) df += SG6;
462       if(df == SG3) break;
463       if(S3 == 3) switch(df) {
464         case 0:
465         case 4:
466         case 5:
467           at = at1;
468           continue;
469         case 2: {
470           conn(at, dx+1);
471           wc.rdir = (dx+1) % 6;
472           break;
473           }
474         case 1: {
475           auto at2 = at + eudir(dx+1);
476           auto& wc2 = get_mapping(at2);
477           if(wc2.cw.at) { at = at1; continue; }
478           wc.rdir = (dx+1) % 6;
479           conn(at, (dx+1) % 6);
480           conn(at1, (dx+2) % 6);
481           conn(at2, (dx+0) % 6);
482           wc1.rdir = -1;
483           wc2.rdir = dx;
484           break;
485           }
486         default:
487           println(hlog, "case unhandled ", df);
488           exit(1);
489         }
490       else switch(df) {
491         case 0:
492         case 3:
493           at = at1;
494           continue;
495         case 1:
496           auto at2 = at + eudir(dx+1);
497           auto& wc2 = get_mapping(at2);
498           if(wc2.cw.at) {
499             auto at3 = at1 + eudir(wc1.rdir);
500             auto& wc3 = get_mapping(at3);
501             auto at4 = at3 + eudir(wc3.rdir);
502             if(at4 == at2) {
503               wc.rdir = (dx+1)%4;
504               wc1.rdir = -1;
505               wc3.rdir = -1;
506               conn(at, (dx+1)%4);
507               }
508             else {
509               at = at1;
510               }
511             }
512           else {
513             wc.rdir = (dx+1)%4;
514             wc1.rdir = -1;
515             wc2.rdir = dx%4;
516             int bdir = -1;
517             int bdist = 100;
518             for(int d=0; d<4; d++) {
519               auto &wcm = get_mapping(at2 + eudir(d));
520               if(wcm.cw.at && length(wcm.start - at2) < bdist)
521                 bdist = length(wcm.start - at2), bdir = d;
522               }
523             if(bdir != -1) conn(at2 + eudir(bdir), bdir ^ 2);
524             conn(at, (dx+1)%4);
525             conn(at2, dx%4);
526 
527             at = param * loc(1,0) + at * loc(0, 1);
528             }
529           break;
530         }
531       }
532 
533     DEBB(DF_GP, ("DONE"))
534     }
535 
loctoh_ort(loc at)536   EX hyperpoint loctoh_ort(loc at) {
537     return point3(at.first, at.second, 1);
538     }
539 
540   hyperpoint corner_coords6[7] = {
541     point3(2, -1, 0),
542     point3(1, 1, 0),
543     point3(-1, 2, 0),
544     point3(-2, 1, 0),
545     point3(-1, -1, 0),
546     point3(1, -2, 0),
547     point3(0, 0, 0) // center, not a corner
548     };
549 
550   hyperpoint corner_coords4[7] = {
551     point3(1.5, -1.5, 0),
552 //    point3(1, 0, 0),
553     point3(1.5, 1.5, 0),
554 //    point3(0, 1, 0),
555     point3(-1.5, 1.5, 0),
556 //    point3(-1, 0, 0),
557     point3(-1.5, -1.5, 0),
558 //    point3(0, -1, 0),
559     point3(0, 0, 0),
560     point3(0, 0, 0),
561     point3(0, 0, 0)
562     };
563 
564   #define corner_coords (S3==3 ? corner_coords6 : corner_coords4)
565 
cornmul(const transmatrix & corners,const hyperpoint & c)566   hyperpoint cornmul(const transmatrix& corners, const hyperpoint& c) {
567     if(sphere && S3 == 3) {
568       ld cmin = c[0] * c[1] * c[2] * (6 - S7);
569       return corners * point3(c[0] + cmin, c[1] + cmin, c[2] + cmin);
570       }
571     else return corners * c;
572     }
573 
atz(const transmatrix & T,const transmatrix & corners,loc at,int cornerid=6,ld cf=3)574   hyperpoint atz(const transmatrix& T, const transmatrix& corners, loc at, int cornerid = 6, ld cf = 3) {
575     int sp = 0;
576     again:
577     auto corner = corners * (loctoh_ort(at) + (corner_coords[cornerid] / cf));
578     if(corner[1] < -1e-6 || corner[2] < -1e-6) {
579       at = at * eudir(1);
580       if(cornerid < SG6) cornerid = (1 + cornerid) % SG6;
581       sp++;
582       goto again;
583       }
584     if(sp>SG3) sp -= SG6;
585 
586     return normalize(spin(2*M_PI*sp/S7) * cornmul(T, corner));
587     }
588 
dir_matrix(int i)589   transmatrix dir_matrix(int i) {
590     auto ddspin = [] (int d) -> transmatrix {
591       return spin(M_PI - d * 2 * M_PI / S7 - cgi.hexshift);
592       };
593     return spin(-cgi.gpdata->alpha) * build_matrix(
594       C0,
595       ddspin(i) * xpush0(cgi.tessf),
596       ddspin(i+1) * xpush0(cgi.tessf),
597       C03
598       );
599     }
600 
prepare_matrices()601   void prepare_matrices() {
602     cgi.gpdata->corners = inverse(build_matrix(
603       loctoh_ort(loc(0,0)),
604       loctoh_ort(param),
605       loctoh_ort(param * loc(0,1)),
606       C03
607       ));
608     cgi.gpdata->Tf.resize(S7);
609     for(int i=0; i<S7; i++) {
610       transmatrix T = dir_matrix(i);
611       for(int x=-GOLDBERG_LIMIT_HALF; x<GOLDBERG_LIMIT_HALF; x++)
612       for(int y=-GOLDBERG_LIMIT_HALF; y<GOLDBERG_LIMIT_HALF; y++)
613       for(int d=0; d<(S3==3?6:4); d++) {
614         loc at = loc(x, y);
615 
616         hyperpoint h = atz(T, cgi.gpdata->corners, at, 6);
617         hyperpoint hl = atz(T, cgi.gpdata->corners, at + eudir(d), 6);
618         cgi.gpdata->Tf[i][x&GOLDBERG_MASK][y&GOLDBERG_MASK][d] = rgpushxto0(h) * rspintox(gpushxto0(h) * hl) * spin(M_PI);
619         }
620       }
621     }
622 
623   EX hyperpoint get_corner_position(const local_info& li, int cid, ld cf IS(3)) {
624     int i = li.last_dir;
625     if(i == -1)
626       return atz(dir_matrix(cid), cgi.gpdata->corners, li.relative, 0, cf);
627     else {
628       auto& cellmatrix = cgi.gpdata->Tf[i][li.relative.first&GOLDBERG_MASK][li.relative.second&GOLDBERG_MASK][fixg6(li.total_dir)];
629       return inverse(cellmatrix) * atz(dir_matrix(i), cgi.gpdata->corners, li.relative, fixg6(cid + li.total_dir), cf);
630       }
631     }
632 
633   EX hyperpoint get_corner_position(cell *c, int cid, ld cf IS(3)) {
634     return get_corner_position(get_local_info(c), cid, cf);
635     }
636 
637   map<pair<int, int>, loc> center_locs;
638 
compute_geometry(bool inv)639   EX void compute_geometry(bool inv) {
640     center_locs.clear();
641     if(GOLDBERG_INV || inv) {
642       if(!cgi.gpdata) cgi.gpdata = make_shared<geometry_information::gpdata_t>();
643       gp::clear_plainshapes();
644       int x = param.first;
645       int y = param.second;
646       if(S3 == 3)
647         cgi.gpdata->area = ((2*x+y) * (2*x+y) + y*y*3) / 4;
648       else
649         cgi.gpdata->area = x * x + y * y;
650       next = point3(x+y/2., -y * sqrt(3) / 2, 0);
651       ld scale = 1 / hypot_d(2, next);
652       if(!GOLDBERG) scale = 1;
653       cgi.crossf *= scale;
654       cgi.hepvdist *= scale;
655       cgi.hexhexdist *= scale;
656       cgi.hexvdist *= scale;
657       cgi.rhexf *= scale;
658 //    spin = spintox(next);
659 //    ispin = rspintox(next);
660       cgi.gpdata->alpha = -atan2(next[1], next[0]) * 6 / S7;
661       if(S3 == 3)
662         cgi.base_distlimit = (cgi.base_distlimit + log(scale) / log(2.618)) / scale;
663       else
664         cgi.base_distlimit = 3 * max(param.first, param.second) + 2 * min(param.first, param.second);
665       if(S7 == 12)
666         cgi.base_distlimit = 2 * param.first + 2 * param.second + 1;
667       if(cgi.base_distlimit > SEE_ALL)
668         cgi.base_distlimit = SEE_ALL;
669       prepare_matrices();
670       DEBB(DF_GEOM | DF_POLY, ("scale = ", scale));
671       }
672     }
673 
674   loc config;
675 
internal_representation(loc v)676   loc internal_representation(loc v) {
677     int& x = v.first, &y = v.second;
678     while(x < 0 || y < 0 || (x == 0 && y > 0))
679       v = v * loc(0, 1);
680     if(x > 8) x = 8;
681     if(y > 8) y = 8;
682     if(S3 == 3 && y > x) v = v * loc(1, -1);
683     return v;
684     }
685 
human_representation(loc v)686   EX loc human_representation(loc v) {
687     int& x = v.first, &y = v.second;
688     if(S3 == 3) while(x < 0 || y < 0 || (x == 0 && y > 0))
689       v = v * loc(0, 1);
690     return v;
691     }
692 
variation_for(loc xy)693   EX eVariation variation_for(loc xy) {
694     if(xy.first == 1 && xy.second == 0)
695       return eVariation::pure;
696     if(xy.first == 1 && xy.second == 1 && S3 == 3)
697       return eVariation::bitruncated;
698     return eVariation::goldberg;
699     }
700 
whirl_set(loc xy)701   void whirl_set(loc xy) {
702     xy = internal_representation(xy);
703     if(xy.second && xy.second != xy.first && nonorientable) {
704       addMessage(XLAT("This does not work in non-orientable geometries"));
705       xy.second = 0;
706       }
707     config = human_representation(xy);
708     auto g = screens;
709     if(xy.first == 0 && xy.second == 0) xy.first = 1;
710     stop_game();
711     param = xy;
712     if(xy.first == 1 && xy.second == 0) {
713       set_variation(eVariation::pure);
714       }
715     else if(xy.first == 1 && xy.second == 1 && S3 == 3) {
716       set_variation(eVariation::bitruncated);
717       }
718     else
719       set_variation(eVariation::goldberg);
720     start_game();
721     screens = g;
722     }
723 
helptext()724   string helptext() {
725     return XLAT(
726       "Goldberg polyhedra are obtained by adding extra hexagons to a dodecahedron. "
727       "GP(x,y) means that, to get to a nearest non-hex from any non-hex, you should move x "
728       "cells in any direction, turn right 60 degrees, and move y cells. "
729       "HyperRogue generalizes this to any tesselation with 3 faces per vertex. "
730       "By default HyperRogue uses bitruncation, which corresponds to GP(1,1)."
731       );
732     }
733 
show()734   void show() {
735     cmode = sm::SIDE | sm::MAYDARK;
736     gamescreen(0);
737     dialog::init(XLAT("variations"));
738 
739     int min_quality_chess = 0;
740 
741     int min_quality = 0;
742 #if CAP_TEXTURE
743     if((texture::config.tstate == texture::tsActive) && (S7 % 2 == 1)) {
744       if(texture::cgroup == cpFootball || texture::cgroup == cpThree) min_quality = 1;
745       }
746 
747     if((texture::config.tstate == texture::tsActive) && (S7 % 2 == 1) && (S3 == 4)) {
748       if(texture::cgroup == cpChess) min_quality = 1;
749       }
750 #endif
751     if(min_quality == 0 && min_quality_chess == 0) {
752       dialog::addBoolItem(XLAT("pure"), PURE || (GOLDBERG && univ_param() == loc(1,0)), 'a');
753       dialog::lastItem().value = "GP(1,0)";
754       dialog::add_action_confirmed([] { whirl_set(loc(1, 0)); });
755       }
756 
757     if(min_quality_chess == 0) {
758       dialog::addBoolItem(XLAT("bitruncated"), BITRUNCATED, 'b');
759       dialog::add_action_confirmed([] {
760         if(S3 == 4) {
761           if(!BITRUNCATED) {
762             stop_game();
763             set_variation(eVariation::bitruncated);
764             start_game();
765             }
766           }
767         else
768           whirl_set(loc(1, 1));
769         });
770       }
771 
772     dialog::lastItem().value = S3 == 3 ? "GP(1,1)" : ONOFF(BITRUNCATED);
773 
774     if(min_quality == 0 || min_quality_chess) {
775       dialog::addBoolItem(S3 == 3 ? XLAT("chamfered") : XLAT("expanded"), univ_param() == loc(2,0) && GOLDBERG, 'c');
776       dialog::lastItem().value = "GP(2,0)";
777       dialog::add_action_confirmed([] {
778         whirl_set(loc(2, 0));
779         });
780       }
781 
782     if(S3 == 3) {
783       dialog::addBoolItem(XLAT("2x bitruncated"), GOLDBERG && univ_param() == loc(3,0), 'd');
784       dialog::lastItem().value = "GP(3,0)";
785       dialog::add_action_confirmed([] {
786         whirl_set(loc(3, 0));
787         });
788       }
789     else {
790       dialog::addBoolItem(XLAT("rectified"), param == loc(1,1) && GOLDBERG, 'd');
791       dialog::lastItem().value = "GP(1,1)";
792       dialog::add_action_confirmed([] {
793         whirl_set(loc(1, 1));
794         });
795       }
796 
797     dialog::addBreak(100);
798     dialog::addSelItem("x", its(config.first), 'x');
799     dialog::add_action([] { dialog::editNumber(config.first, 0, 8, 1, 1, "x", helptext()); });
800     dialog::addSelItem("y", its(config.second), 'y');
801     dialog::add_action([] { dialog::editNumber(config.second, 0, 8, 1, 1, "y", helptext()); });
802 
803     if(config.second && config.second != config.first && nonorientable) {
804       dialog::addInfo(XLAT("This does not work in non-orientable geometries"));
805       }
806     else if((config.first-config.second)%3 && min_quality)
807       dialog::addInfo(XLAT("This pattern needs x-y divisible by 3"));
808     else if((config.first-config.second)%2 && min_quality_chess)
809       dialog::addInfo(XLAT("This pattern needs x-y divisible by 2"));
810     else {
811       dialog::addBoolItem(XLAT("select"), param == internal_representation(config) && !IRREGULAR && !INVERSE, 'f');
812       dialog::lastItem().value = "GP(x,y)";
813       }
814     dialog::add_action_confirmed([] { whirl_set(config); });
815 
816     dialog::addBreak(100);
817 
818     #if CAP_IRR
819     if(irr::supports(geometry)) {
820       dialog::addBoolItem(XLAT("irregular"), IRREGULAR, 'i');
821       dialog::add_action(dialog::add_confirmation([=] () {
822         if(min_quality && !irr::bitruncations_requested) irr::bitruncations_requested++;
823         if(euclid && (!bounded || nonorientable)) {
824           println(hlog, XLAT("To create Euclidean irregular tesselations, first enable a torus"));
825           return;
826           }
827         if(!IRREGULAR) irr::visual_creator();
828         }));
829       }
830     #endif
831 
832     dialog::addBreak(100);
833     int style = 0;
834     auto v0 = variation_for(param);
835     bool bad_bi = BITRUNCATED && a4;
836     if(!bad_bi) {
837       dynamicval<eVariation> v(variation, v0);
838       if(geosupport_football() == 2) style = 3;
839       if(geosupport_chessboard()) style = 2;
840       }
841     if(style == 2) {
842       dialog::addBoolItem(XLAT("inverse rectify"), UNRECTIFIED, 'r');
843       dialog::add_action_confirmed([v0] {
844         param = univ_param();
845         if(UNRECTIFIED) set_variation(v0);
846         else set_variation(eVariation::unrectified);
847         start_game();
848         config = human_representation(univ_param());
849         });
850       }
851     else if(style == 3) {
852       dialog::addBoolItem(XLAT("inverse truncate"), UNTRUNCATED, 't');
853       dialog::add_action_confirmed([v0] {
854         param = univ_param();
855         if(UNTRUNCATED) set_variation(v0);
856         else set_variation(eVariation::untruncated);
857         start_game();
858         });
859       dialog::addBoolItem(XLAT("warped version"), WARPED, 'w');
860       dialog::add_action_confirmed([v0] {
861         param = univ_param();
862         if(WARPED) set_variation(v0);
863         else set_variation(eVariation::warped);
864         start_game();
865         });
866       }
867 
868     dialog::addBreak(100);
869     dialog::addItem(XLAT("swap x and y"), 'z');
870     dialog::add_action([] { swap(config.first, config.second); });
871 
872     bool have_dual = !bad_bi && !IRREGULAR && !WARPED;
873     if(S3 == 3 && UNTRUNCATED && (univ_param()*loc(1,1)) % 3) have_dual = false;
874     if(S3 == 4 && UNRECTIFIED && (univ_param()*loc(1,1)) % 2) have_dual = false;
875 
876     if(have_dual) {
877       dialog::addItem(XLAT("dual of current"), 'D');
878       dialog::add_action([] {
879         auto p = univ_param();
880         if(S3 == 3 && !UNTRUNCATED) {
881           println(hlog, "set param to ", p * loc(1,1));
882           whirl_set(p * loc(1, 1));
883           set_variation(eVariation::untruncated);
884           start_game();
885           config = human_representation(univ_param());
886           }
887         else if(S3 == 4 && !UNRECTIFIED) {
888           whirl_set(p * loc(1, 1));
889           set_variation(eVariation::unrectified);
890           start_game();
891           config = human_representation(univ_param());
892           }
893         else if(S3 == 3 && UNTRUNCATED) {
894           println(hlog, "whirl_set to ", (p * loc(1,1)) / 3);
895           whirl_set((p * loc(1,1)) / 3);
896           config = human_representation(univ_param());
897           }
898         else if(S3 == 4 && UNRECTIFIED) {
899           whirl_set((p * loc(1,1)) / 2);
900           config = human_representation(univ_param());
901           }
902         });
903       }
904 
905     dialog::addBreak(100);
906     dialog::addHelp();
907     dialog::add_action([] { gotoHelp(helptext()); });
908     dialog::addBack();
909     dialog::display();
910     }
911 
univ_param()912   EX loc univ_param() {
913     if(GOLDBERG_INV) return param;
914     else if(PURE) return loc(1,0);
915     else return loc(1,1);
916     }
917 
configure()918   EX void configure() {
919     auto l = univ_param();
920     param = l;
921     config = human_representation(l);
922     pushScreen(gp::show);
923     }
924 
be_in_triangle(local_info & li)925   EX void be_in_triangle(local_info& li) {
926     int sp = 0;
927     auto& at = li.relative;
928     again:
929     auto corner = cgi.gpdata->corners * loctoh_ort(at);
930     if(corner[1] < -1e-6 || corner[2] < -1e-6) {
931       at = at * eudir(1);
932       sp++;
933       goto again;
934       }
935     if(sp>SG3) sp -= SG6;
936     li.last_dir = gmod(li.last_dir - sp, S7);
937     }
938 
939   // from some point X, (0,0) is in distance dmain, param is in distance d0, and param*z is in distance d1
940   // what is the distance of at from X?
941 
solve_triangle(int dmain,int d0,int d1,loc at)942   EX int solve_triangle(int dmain, int d0, int d1, loc at) {
943     loc centerloc(0, 0);
944     auto rel = make_pair(d0-dmain, d1-dmain);
945     if(center_locs.count(rel))
946       centerloc = center_locs[rel];
947     else {
948       bool found = false;
949       for(int y=-20; y<=20; y++)
950       for(int x=-20; x<=20; x++) {
951         loc c(x, y);
952         int cc = length(c);
953         int c0 = length(c - param);
954         int c1 = length(c - param*loc(0,1));
955         if(c0-cc == d0-dmain && c1-cc == d1-dmain)
956           found = true, centerloc = c;
957         }
958       if(!found && !quotient) {
959         println(hlog, "Warning: centerloc not found: ", make_tuple(dmain, d0, d1));
960         }
961       center_locs[rel] = centerloc;
962       }
963 
964     return dmain + length(centerloc-at) - length(centerloc);
965     }
966 
solve_quad(int dmain,int d0,int d1,int dx,loc at)967   int solve_quad(int dmain, int d0, int d1, int dx, loc at) {
968     loc centerloc(0, 0);
969     auto rel = make_pair(d0-dmain, (d1-dmain) + 1000 * (dx-dmain) + 1000000);
970     if(center_locs.count(rel))
971       centerloc = center_locs[rel];
972     else {
973       bool found = false;
974       for(int y=-20; y<=20; y++)
975       for(int x=-20; x<=20; x++) {
976         loc c(x, y);
977         int cc = length(c);
978         int c0 = length(c - param);
979         int c1 = length(c - param*loc(0,1));
980         int c2 = length(c - param*loc(1,1));
981         if(c0-cc == d0-dmain && c1-cc == d1-dmain && c2-cc == dx-dmain)
982           found = true, centerloc = c;
983         }
984       if(!found && !quotient) {
985         println(hlog, "Warning: centerloc not found: ", make_tuple(dmain, d0, d1, dx));
986         }
987       center_locs[rel] = centerloc;
988       }
989 
990     return dmain + length(centerloc-at) - length(centerloc);
991     }
992 
get_master_coordinates(cell * c)993   EX hyperpoint get_master_coordinates(cell *c) {
994     auto li = get_local_info(c);
995     be_in_triangle(li);
996     return cgi.gpdata->corners * loctoh_ort(li.relative);
997     }
998 
compute_dist(cell * c,int master_function (cell *))999   EX int compute_dist(cell *c, int master_function(cell*)) {
1000     if(!GOLDBERG) return master_function(c);
1001     auto li = get_local_info(c);
1002     be_in_triangle(li);
1003 
1004     cell *cm = c->master->c7;
1005 
1006     int i = li.last_dir;
1007     auto at = li.relative;
1008 
1009     auto dmain = master_function(cm);
1010     auto d0 = master_function(createStep(cm->master, i)->c7);
1011     auto d1 = master_function(createStep(cm->master, cm->c.fix(i+1))->c7);
1012 
1013     if(S3 == 4) {
1014       heptspin hs(cm->master, i);
1015       hs += wstep; hs+=-1; hs += wstep;
1016       auto d2 = master_function(hs.at->c7);
1017       return solve_quad(dmain, d0, d1, d2, at);
1018       }
1019 
1020     return solve_triangle(dmain, d0, d1, at);
1021     }
1022 
dist_2()1023   EX int dist_2() {
1024     return length(univ_param());
1025     }
1026 
dist_3()1027   EX int dist_3() {
1028     return length(univ_param() * loc(1,1));
1029     }
1030 
dist_1()1031   EX int dist_1() {
1032     return dist_3() - dist_2();
1033     }
1034 #else
dist_1()1035   EX int dist_1() { return 1; }
dist_2()1036   EX int dist_2() { return BITRUNCATED ? 2 : 1; }
dist_3()1037   EX int dist_3() { return BITRUNCATED ? 3 : 2; }
1038 #endif
1039 
get_masters(cell * c)1040   EX array<heptagon*, 3> get_masters(cell *c) {
1041     if(0);
1042     #if CAP_GP
1043     else if(INVERSE) {
1044       c = get_mapped(c);
1045       return UIU(get_masters(c));
1046       }
1047     else if(GOLDBERG) {
1048       auto li = get_local_info(c);
1049       be_in_triangle(li);
1050       auto cm = c->master;
1051       int i = li.last_dir;
1052       return make_array(cm, cm->cmove(i), cm->cmodmove(i+1));
1053       }
1054     #endif
1055     #if CAP_IRR
1056     else if(IRREGULAR)
1057       return irr::get_masters(c);
1058     #endif
1059     else
1060       return make_array(c->cmove(0)->master, c->cmove(2)->master, c->cmove(4)->master);
1061     }
1062 
operation_name()1063   EX string operation_name() {
1064     if(0);
1065     #if CAP_IRR
1066     else if(IRREGULAR)
1067       return XLAT("irregular");
1068     #endif
1069     else if(DUAL)
1070       return XLAT("dual");
1071     else if(PURE)
1072       return XLAT("pure");
1073     else if(BITRUNCATED)
1074       return XLAT("bitruncated");
1075     #if CAP_GP
1076     else if(GOLDBERG && param == loc(1, 0))
1077       return XLAT("pure");
1078     else if(GOLDBERG && param == loc(1, 1) && S3 == 3)
1079       return XLAT("bitruncated");
1080     else if(GOLDBERG && param == loc(1, 1) && S3 == 4)
1081       return XLAT("rectified");
1082     else if(UNRECTIFIED && param == loc(1, 1) && S3 == 4)
1083       return XLAT("dual");
1084     else if(UNTRUNCATED && param == loc(1, 1) && S3 == 3)
1085       return XLAT("dual");
1086     else if(GOLDBERG && param == loc(2, 0))
1087       return S3 == 3 ? XLAT("chamfered") : XLAT("expanded");
1088     else if(GOLDBERG && param == loc(3, 0) && S3 == 3)
1089       return XLAT("2x bitruncated");
1090     else if(variation == eVariation::subcubes)
1091       return XLAT("subcubed") + "(" + its(reg3::subcube_count) + ")";
1092     else if(variation == eVariation::dual_subcubes)
1093       return XLAT("dual-subcubed") + "(" + its(reg3::subcube_count) + ")";
1094     else if(variation == eVariation::bch)
1095       return XLAT("bitruncated-subcubed") + "(" + its(reg3::subcube_count) + ")";
1096     else if(variation == eVariation::coxeter)
1097       return XLAT("subdivided") + "(" + its(reg3::coxeter_param) + ")";
1098     else {
1099       auto p = human_representation(param);
1100       string s = "GP(" + its(p.first) + "," + its(p.second) + ")";
1101       if(UNRECTIFIED) return XLAT("unrectified") + " " + s;
1102       if(WARPED) return XLAT("warped") + " " + s;
1103       if(UNTRUNCATED) return XLAT("untruncated") + " " + s;
1104       return s;
1105       }
1106     #else
1107     else return "UNSUPPORTED";
1108     #endif
1109     }
1110 
1111   /* inverse map */
1112 
1113   EX hrmap *pmap;
1114   // EX geometry_information *underlying_cgip;
1115 
1116   struct hrmap_inverse : hrmap {
1117     hrmap *underlying_map;
1118 
1119     map<cell*, cell*> mapping;
1120     map<cell*, int> shift;
1121 
in_underlyinghr::gp::hrmap_inverse1122     template<class T> auto in_underlying(const T& t) -> decltype(t()) {
1123       dynamicval<hrmap*> gpm(pmap, this);
1124       dynamicval<eVariation> gva(variation, variation_for(param));
1125       dynamicval<hrmap*> gu(currentmap, underlying_map);
1126       // dynamicval<geometry_information*> gc(cgip, underlying_cgip);
1127       return t();
1128       }
1129 
get_mappedhr::gp::hrmap_inverse1130     cell* get_mapped(cell *underlying_cell, int set_shift) {
1131       if(mapping.count(underlying_cell))
1132         return mapping[underlying_cell];
1133       int d = underlying_cell->type;
1134       if(UNTRUNCATED) d /= 2;
1135       if(WARPED && set_shift < 2) d /= 2;
1136       cell *c = newCell(d, underlying_cell->master);
1137       mapping[underlying_cell] = c;
1138       if(!UNRECTIFIED) shift[c] = set_shift;
1139       mapping[c] = underlying_cell;
1140       return c;
1141       }
1142 
relative_matrixhhr::gp::hrmap_inverse1143     transmatrix relative_matrixh(heptagon *h2, heptagon *h1, const hyperpoint& hint) override {
1144       return in_underlying([&] { return currentmap->relative_matrix(h2, h1, hint); });
1145       }
1146 
relative_matrixchr::gp::hrmap_inverse1147     transmatrix relative_matrixc(cell *c2, cell *c1, const hyperpoint& hint) override {
1148       c1 = mapping[c1];
1149       c2 = mapping[c2];
1150       return in_underlying([&] { return currentmap->relative_matrix(c2, c1, hint); });
1151       }
1152 
~hrmap_inversehr::gp::hrmap_inverse1153     ~hrmap_inverse() {
1154       in_underlying([this] { delete underlying_map; });
1155       }
1156 
getOriginhr::gp::hrmap_inverse1157     heptagon *getOrigin() override { return in_underlying([this] { return underlying_map->getOrigin(); }); }
1158 
1159     cell *gs;
1160 
gamestarthr::gp::hrmap_inverse1161     cell* gamestart() override {
1162       return gs;
1163       }
1164 
hrmap_inversehr::gp::hrmap_inverse1165     hrmap_inverse() {
1166       if(0) {
1167         println(hlog, "making ucgi");
1168         dynamicval<eVariation> gva(variation, variation_for(param));
1169         check_cgi();
1170         cgi.require_basics();
1171         // underlying_cgip = cgip;
1172         println(hlog, "done ucgi");
1173         }
1174       bool warped = WARPED;
1175       in_underlying([&,this] {
1176         initcells();
1177         underlying_map = currentmap;
1178         gs = currentmap->gamestart();
1179         if(!warped) gs = gs->cmove(0);
1180         });
1181       if(UNTRUNCATED) gs = get_mapped(gs, 1);
1182       else gs = get_mapped(gs, 2);
1183       for(hrmap*& m: allmaps) if(m == underlying_map) m = NULL;
1184       }
1185 
create_movehr::gp::hrmap_inverse1186     cell *create_move(cell *parent, int d) {
1187       if(UNRECTIFIED) {
1188         cellwalker cw(mapping[parent], d);
1189         in_underlying([&] {
1190           cw += wstep;
1191           cw --;
1192           cw += wstep;
1193           cw --;
1194           });
1195         cw.at = get_mapped(cw.at, 0);
1196         parent->c.connect(d, cw.at, cw.spin, cw.mirrored);
1197         return cw.at;
1198         }
1199       if(UNTRUNCATED) {
1200         cellwalker cw(mapping[parent], 2*d+shift[parent]);
1201         in_underlying([&] {
1202           cw += wstep;
1203           });
1204         cw.at = get_mapped(cw.at, cw.spin & 1);
1205         parent->c.connect(d, cw.at, cw.spin / 2, cw.mirrored);
1206         return cw.at;
1207         }
1208       if(WARPED) {
1209         int sh = shift[parent];
1210         if(sh == 2) {
1211           cellwalker cw(mapping[parent], d);
1212           in_underlying([&] { cw += wstep; });
1213           cw.at = get_mapped(cw.at, cw.spin & 1);
1214           parent->c.connect(d, cw.at, cw.spin / 2, cw.mirrored);
1215           return cw.at;
1216           }
1217         else {
1218           cellwalker cw(mapping[parent], 2*d+sh);
1219           in_underlying([&] {
1220             cw += wstep;
1221             });
1222           cw.at = get_mapped(cw.at, 2);
1223           parent->c.connect(d, cw.at, cw.spin, cw.mirrored);
1224           return cw.at;
1225           }
1226         }
1227       throw hr_exception("unimplemented");
1228       }
1229 
adjhr::gp::hrmap_inverse1230     transmatrix adj(cell *c, int d) override {
1231       transmatrix T;
1232       if(UNRECTIFIED) {
1233         cellwalker cw(mapping[c], d);
1234         in_underlying([&] {
1235           T = currentmap->adj(cw.at, cw.spin);
1236           cw += wstep;
1237           cw --;
1238           T = T * currentmap->adj(cw.at, cw.spin);
1239           });
1240         }
1241       if(UNTRUNCATED) {
1242         cellwalker cw(mapping[c], 2*d+shift[c]);
1243         in_underlying([&] { T = currentmap->adj(cw.at, cw.spin); });
1244         }
1245       if(WARPED) {
1246         int sh = shift[c];
1247         if(sh == 2) {
1248           auto c1 = mapping[c];
1249           in_underlying([&] { T = currentmap->adj(c1, d); });
1250           }
1251         else {
1252           cellwalker cw(mapping[c], 2*d+shift[c]);
1253           in_underlying([&] { T = currentmap->adj(cw.at, cw.spin); });
1254           }
1255         }
1256       return T;
1257       }
1258 
draw_athr::gp::hrmap_inverse1259     void draw_at(cell *at, const shiftmatrix& where) override {
1260 
1261       dq::clear_all();
1262 
1263       auto enqueue = (quotient ? dq::enqueue_by_matrix_c : dq::enqueue_c);
1264       enqueue(at, where);
1265 
1266       while(!dq::drawqueue_c.empty()) {
1267         auto& p = dq::drawqueue_c.front();
1268         cell *c = p.first;
1269         shiftmatrix V = p.second;
1270         auto c1 = get_mapped(c, 0);
1271 
1272         in_underlying([&] {
1273           if(GOLDBERG) {
1274             gp::current_li = gp::get_local_info(c1);
1275             }
1276           else {
1277             gp::current_li.relative.first = shvid(c1);
1278             gp::current_li.relative.second = shift[c];
1279             }
1280           });
1281 
1282 
1283         dq::drawqueue_c.pop();
1284 
1285         if(!do_draw(c, V)) continue;
1286         drawcell(c, V);
1287 
1288         for(int i=0; i<c->type; i++) if(c->cmove(i))
1289           enqueue(c->move(i), optimized_shift(V * adj(c, i)));
1290         }
1291       }
1292 
find_cell_connectionhr::gp::hrmap_inverse1293     void find_cell_connection(cell *c, int d) override {
1294       inverse_move(c, d);
1295       }
1296 
shvidhr::gp::hrmap_inverse1297     int shvid(cell *c) override {
1298       return gp::get_plainshape_id(c);
1299       }
1300 
get_cornerhr::gp::hrmap_inverse1301     hyperpoint get_corner(cell *c, int cid, ld cf) override {
1302       if(UNTRUNCATED) {
1303         cell *c1 = gp::get_mapped(c);
1304         cellwalker cw(c1, cid*2);
1305         if(!gp::untruncated_shift(c)) cw--;
1306         hyperpoint h = UIU(nearcorner(c1, cw.spin));
1307         return mid_at_actual(h, 3/cf);
1308         }
1309       if(UNRECTIFIED) {
1310         cell *c1 = gp::get_mapped(c);
1311         hyperpoint h = UIU(nearcorner(c1, cid));
1312         return mid_at_actual(h, 3/cf);
1313         }
1314       if(WARPED) {
1315         int sh = gp::untruncated_shift(c);
1316         cell *c1 = gp::get_mapped(c);
1317         if(sh == 2) {
1318           cellwalker cw(c, cid);
1319           hyperpoint h1 = UIU(tC0(currentmap->adj(c1, cid)));
1320           cw--;
1321           hyperpoint h2 = UIU(tC0(currentmap->adj(c1, cw.spin)));
1322           hyperpoint h = mid(h1, h2);
1323           return mid_at_actual(h, 3/cf);
1324           }
1325         else {
1326           cellwalker cw(c1, cid*2);
1327           if(!gp::untruncated_shift(c)) cw--;
1328           hyperpoint h = UIU(nearcorner(c1, cw.spin));
1329           h = mid(h, C0);
1330           return mid_at_actual(h, 3/cf);
1331           }
1332         }
1333       return C0;
1334       }
1335     };
1336 
new_inverse()1337   EX hrmap* new_inverse() { return new hrmap_inverse; }
1338 
inv_map()1339   hrmap_inverse* inv_map() { return (hrmap_inverse*)currentmap; }
1340 
get_underlying_map()1341   EX hrmap* get_underlying_map() { return inv_map()->underlying_map; }
get_mapped(cell * c)1342   EX cell* get_mapped(cell *c) { return inv_map()->get_mapped(c, 0); }
untruncated_shift(cell * c)1343   EX int untruncated_shift(cell *c) { return inv_map()->shift[c]; }
1344 
delete_mapped(cell * c)1345   EX void delete_mapped(cell *c) {
1346     if(!pmap) return;
1347     auto i = (hrmap_inverse*) pmap;
1348     if(i->mapping.count(c))
1349       destroy_cell(i->mapping[c]);
1350     }
1351 
inverse_move(cell * c,int d)1352   EX cell *inverse_move(cell *c, int d) { return inv_map()->create_move(c, d); }
1353 
1354   #if HDR
in_underlying_geometry(const T & f)1355   template<class T> auto in_underlying_geometry(const T& f) -> decltype(f()) {
1356     if(!INVERSE) return f();
1357     dynamicval<hrmap*> gpm(pmap, currentmap);
1358     dynamicval<eVariation> gva(variation, variation_for(param));
1359     dynamicval<hrmap*> gu(currentmap, get_underlying_map());
1360     // dynamicval<geometry_information*> gc(cgip, underlying_cgip);
1361     return f();
1362     }
1363 
1364   #define UIU(x) hr::gp::in_underlying_geometry([&] { return (x); })
1365   #endif
1366 
1367 
1368 
1369   }}
1370