1 // Hyperbolic Rogue -- Euclidean geometry
2 // Copyright (C) 2011-2019 Zeno Rogue, see 'hyper.cpp' for details
3 
4 /** \file euclid.cpp
5  *  \brief Euclidean geometry, including 2D, 3D, and quotient spaces
6  */
7 
8 #include "hyper.h"
9 namespace hr {
10 
11 EX namespace euc {
12 
13   #if HDR
14   struct coord : array<int, 3> {
coordhr::euc::coord15     coord() {}
coordhr::euc::coord16     coord(int x, int y, int z) { self[0] = x; self[1] = y; self[2] = z; }
operator +=hr::euc::coord17     coord& operator += (coord b) { for(int i: {0,1,2}) self[i] += b[i]; return self; }
operator -=hr::euc::coord18     coord& operator -= (coord b) { for(int i: {0,1,2}) self[i] -= b[i]; return self; }
operator +hr::euc::coord19     coord operator + (coord b) const { coord a = self; return a += b; }
operator -hr::euc::coord20     coord operator - (coord b) const { coord a = self; return a -= b; }
operator -hr::euc::coord21     coord operator -() const { return coord(-self[0], -self[1], -self[2]); }
operator +hr::euc::coord22     coord& operator +() { return self; }
operator +hr::euc::coord23     const coord& operator +() const { return self; }
operator *hr::euc::coord24     coord operator *(int x) const { return coord(x*self[0], x*self[1], x*self[2]); }
operator *(int x,const coord & y)25     friend coord operator *(int x, const coord& y) { return coord(x*y[0], x*y[1], x*y[2]); }
26     };
27 
28   typedef array<coord, 3> intmatrix;
29   #endif
30 
31   EX const coord euzero = coord(0,0,0);
32   EX const coord eutester = coord(3,7,0);
33   EX intmatrix euzeroall = make_array<coord>(euzero, euzero, euzero);
34 
35   static const intmatrix main_axes = make_array<coord>(coord(1,0,0), coord(0,1,0), coord(0,0,1));
36 
get_shifttable()37   EX vector<coord> get_shifttable() {
38     static const coord D0 = main_axes[0];
39     static const coord D1 = main_axes[1];
40     static const coord D2 = main_axes[2];
41     vector<coord> shifttable;
42 
43     /* for portal spaces... */
44     auto g = geometry;
45     if(S7 == 6 && WDIM == 3) g = gCubeTiling;
46 
47     switch(g) {
48       case gCubeTiling:
49         shifttable = { +D0, +D1, +D2 };
50         break;
51 
52       case gRhombic3:
53         shifttable = { D0+D1, D0+D2, D1+D2, D1-D2, D0-D2, D0-D1 };
54         break;
55 
56       case gBitrunc3:
57         shifttable = { 2*D0, 2*D1, 2*D2, D0+D1+D2, D0+D1-D2, D0-D1-D2, D0-D1+D2 };
58         break;
59 
60       case gEuclid:
61         shifttable = { D0, D1, D1-D0, -D0, -D1, D0-D1 };
62         break;
63 
64       case gEuclidSquare:
65         shifttable = { D0, D1, -D0, -D1 };
66         break;
67 
68       default:
69         printf("euc::get_shifttable() called in geometry that is not euclid3");
70         exit(1);
71       }
72 
73     // reverse everything
74     int s = isize(shifttable);
75     for(int i=0; i<s; i++) shifttable.push_back(-shifttable[i]);
76     return shifttable;
77     }
78 
79   EX coord basic_canonicalize(coord x);
80 
81   #if HDR
82   struct torus_config {
83     /** periods entered by the user */
84     intmatrix user_axes;
85     /** OR'ed flags: 1 -- flip X in 3D, 2 -- flip Y in 3D, 4 -- flip X/Y in 3D, 8 -- Klein bottle in 2D, 16 -- third turn in 3D, 32 -- Hantzsche-Wendt in 3D */
86     int twisted;
87 
torus_confighr::euc::torus_config88     torus_config() {}
torus_confighr::euc::torus_config89     torus_config(intmatrix user_axes, int twisted) : user_axes(user_axes), twisted(twisted) {}
90     };
91 
92   struct torus_config_full : torus_config {
93     /** optimal representation of the periods */
94     intmatrix optimal_axes;
95     /** regular axes (?) */
96     intmatrix regular_axes;
97     /** in 2D: the period vector which is reflected */
98     gp::loc twisted_vec;
99     /** in 2D: a vector orthogonal to twisted_vec */
100     gp::loc ortho_vec;
101     /** determinant */
102     int det;
103     /** the number of infinite dimensions */
104     int infinite_dims;
105     /** ? */
106     intmatrix inverse_axes;
107     /** for canonicalization on tori */
108     map<coord, int> hash;
109     vector<coord> seq;
110     int index;
111 
resethr::euc::torus_config_full112     void reset() { index = 0; hash.clear(); seq.clear(); }
113 
114     /** add to the tori canonicalization list */
115     void add(coord val);
116 
117     /** get the representative on the tori canonicalization list */
118     coord get(coord x);
119 
120     /** find the equivalence class of coo */
121     coord compute_cat(coord coo);
122 
123     /** canonicalize coord x; in case of twisting, adjust d, M, and mirr accordingly */
124     void canonicalize(coord& x, coord& d, transmatrix& M, bool& mirr);
125     };
126   #endif
127   EX torus_config eu_input, eu_edit;
128   EX torus_config_full eu;
129 
130   struct hrmap_euclidean : hrmap_standard {
131     vector<coord> shifttable;
132     vector<transmatrix> tmatrix;
133     map<coord, heptagon*> spacemap;
134     map<heptagon*, coord> ispacemap;
135     cell *camelot_center;
136 
137     map<gp::loc, struct cdata> eucdata;
138 
compute_tmatrixhr::euc::hrmap_euclidean139     void compute_tmatrix() {
140       shifttable = get_shifttable();
141       tmatrix.resize(S7);
142       for(int i=0; i<S7; i++)
143         tmatrix[i] = eumove(shifttable[i]);
144       }
145 
on_dim_changehr::euc::hrmap_euclidean146     void on_dim_change() override {
147       compute_tmatrix();
148       }
149 
150     vector<cell*> toruscells;
allcellshr::euc::hrmap_euclidean151     vector<cell*>& allcells() override {
152       if(bounded) {
153         if(isize(toruscells) == 0) {
154           celllister cl(getOrigin()->c7, 1000, 1000000, NULL);
155           toruscells = cl.lst;
156           }
157         return toruscells;
158         }
159       return hrmap::allcells();
160       }
161 
hrmap_euclideanhr::euc::hrmap_euclidean162     hrmap_euclidean() {
163       compute_tmatrix();
164       camelot_center = NULL;
165       build_torus3(geometry);
166       #if CAP_IRR
167       if(!valid_irr_torus()) {
168         addMessage(XLAT("Error: period mismatch"));
169         eu_input = irr::base_config;
170         build_torus3(geometry);
171         }
172       #endif
173       }
174 
getOriginhr::euc::hrmap_euclidean175     heptagon *getOrigin() override {
176       return get_at(euzero);
177       }
178 
get_athr::euc::hrmap_euclidean179     heptagon *get_at(coord at) {
180       if(spacemap.count(at))
181         return spacemap[at];
182       else {
183         auto h = init_heptagon(S7);
184         if(!IRREGULAR)
185           h->c7 = newCell(S7, h);
186         #if CAP_IRR
187         else {
188           coord m0 = shifttable[0];
189           transmatrix dummy;
190           bool mirr;
191           auto ati = at;
192           irr::base_config.canonicalize(ati, m0, dummy, mirr);
193           indenter id(2);
194           for(int spin=0; spin<S7; spin++) if(shifttable[spin] == m0) {
195             irr::link_to_base(h, heptspin(((hrmap_euclidean*)irr::base)->get_at(ati), spin, mirr));
196             break;
197             }
198           }
199         #endif
200         if(S7 != 14)
201           h->zebraval = gmod(at[0] + at[1] * 2 + at[2] * 4, 5);
202         else
203           h->zebraval = at[0] & 1;
204         spacemap[at] = h;
205         ispacemap[h] = at;
206 
207         return h;
208         }
209       }
210 
create_stephr::euc::hrmap_euclidean211     heptagon *create_step(heptagon *parent, int d) override {
212       int d1 = (d+S7/2)%S7;
213       bool mirr = false;
214       transmatrix I;
215       auto v = ispacemap[parent] + shifttable[d];
216       auto st = shifttable[d1];
217       eu.canonicalize(v, st, I, mirr);
218       if(eu.twisted)
219         for(int i=0; i<S7; i++) if(shifttable[i] == st) d1 = i;
220 
221       heptagon *h = get_at(v);
222       h->c.connect(d1, parent, d, mirr);
223       return h;
224       }
225 
adjhr::euc::hrmap_euclidean226     transmatrix adj(heptagon *h, int i) override {
227       if(!eu.twisted) return tmatrix[i];
228       transmatrix res = tmatrix[i];
229       coord id = ispacemap[h];
230       id += shifttable[i];
231       auto dummy = euzero;
232       bool dm = false;
233       eu.canonicalize(id, dummy, res, dm);
234 
235       return res;
236       }
237 
adjhr::euc::hrmap_euclidean238     transmatrix adj(cell *c, int i) override {
239       if(WDIM == 3) return adj(c->master, i);
240       else return hrmap_standard::adj(c, i);
241       }
242 
draw_athr::euc::hrmap_euclidean243     void draw_at(cell *at, const shiftmatrix& where) override {
244       dq::clear_all();
245       dq::enqueue_by_matrix(at->master, where * master_relative(centerover, true));
246 
247       while(!dq::drawqueue.empty()) {
248         auto& p = dq::drawqueue.front();
249         heptagon *h = p.first;
250         shiftmatrix V = p.second;
251         dq::drawqueue.pop();
252 
253         cell *c = h->c7;
254 
255         bool draw = drawcell_subs(c, V * spin(master_to_c7_angle()));
256         if(in_wallopt() && isWall3(c) && isize(dq::drawqueue) > 1000 && !hybrid::pmap) continue;
257 
258         if(draw) for(int i=0; i<S7; i++)
259           dq::enqueue_by_matrix(h->move(i), optimized_shift(V * adj(h, i)));
260         }
261       }
262 
relative_matrixhhr::euc::hrmap_euclidean263     transmatrix relative_matrixh(heptagon *h2, heptagon *h1, const hyperpoint& hint) override {
264       if(eu.twisted) {
265         if(h1 == h2) return Id;
266         for(int s=0; s<S7; s++) if(h2 == h1->move(s)) return adj(h1, s);
267         coord c1 = ispacemap[h1];
268         coord c2 = ispacemap[h2];
269         transmatrix T = eumove(c2 - c1);
270 
271         transmatrix I = Id;
272         coord cs = c1;
273         for(int s=0; s<4; s++) {
274           for(int a=-1; a<=1; a++)
275           for(int b=-1; b<=1; b++) {
276             if(b && WDIM == 2) continue;
277             transmatrix T1 = I * eumove((c2 - cs) + a*eu.user_axes[0] + b*eu.user_axes[1]);
278             if(hdist(tC0(T1), hint) < hdist(tC0(T), hint))
279               T = T1;
280             }
281           auto co = eu.user_axes[WDIM-1];
282           cs += co;
283           I = I * eumove(co);
284           auto dummy = euzero;
285           bool dm = false;
286           eu.canonicalize(cs, dummy, I, dm);
287           }
288         return T;
289         }
290       auto d = ispacemap[h2] - ispacemap[h1];
291       d = basic_canonicalize(d);
292       return eumove(d);
293       }
294 
get_cellshapehr::euc::hrmap_euclidean295     subcellshape& get_cellshape(cell* c) override {
296       return *cgi.heptshape;
297       }
298     };
299 
cubemap()300   hrmap_euclidean* cubemap() {
301     if(fake::in()) return FPIU(cubemap());
302     return ((hrmap_euclidean*) currentmap);
303     }
304 
eucmap()305   hrmap_euclidean* eucmap() {
306     return cubemap();
307     }
308 
get_current_shifttable()309   EX vector<coord>& get_current_shifttable() { return cubemap()->shifttable; }
get_spacemap()310   EX map<coord, heptagon*>& get_spacemap() { return cubemap()->spacemap; }
get_ispacemap()311   EX map<heptagon*, coord>& get_ispacemap() { return cubemap()->ispacemap; }
get_camelot_center()312   EX cell *& get_camelot_center() { return cubemap()->camelot_center; }
313 
get_at(coord co)314   EX heptagon* get_at(coord co) { return cubemap()->get_at(co); }
315 
new_map()316   EX hrmap* new_map() {
317     return new hrmap_euclidean;
318     }
319 
move_matrix(heptagon * h,int i)320   EX transmatrix move_matrix(heptagon *h, int i) {
321     return cubemap()->adj(h, i);
322     }
323 
pseudohept(cell * c)324   EX bool pseudohept(cell *c) {
325     if(cgflags & qPORTALSPACE) return false;
326     coord co = cubemap()->ispacemap[c->master];
327     if(S7 == 12) {
328       for(int i=0; i<3; i++) if((co[i] & 1)) return false;
329       }
330     else {
331       for(int i=0; i<3; i++) if(!(co[i] & 1)) return false;
332       }
333     return true;
334     }
335 
dist_alt(cell * c)336   EX int dist_alt(cell *c) {
337     if(WDIM == 2) {
338       auto v = full_coords2(c);
339       return euclidAlt(v.first, v.second);
340       }
341     if(specialland == laCamelot) return dist_relative(c) + roundTableRadius(c);
342     auto v = cubemap()->ispacemap[c->master];
343     if(S7 == 6) return v[2];
344     else if(S7 == 12) return (v[0] + v[1] + v[2]) / 2;
345     else return v[2]/2;
346     }
347 
get_emerald(cell * c)348   EX bool get_emerald(cell *c) {
349     auto v = cubemap()->ispacemap[c->master];
350     int s0 = 0, s1 = 0;
351     for(int i=0; i<3; i++) {
352       v[i] = gmod(v[i], 6);
353       int d = min(v[i], 6-v[i]);;
354       s0 += min(v[i], 6-v[i]);
355       s1 += 3-d;
356       }
357     if(s0 == s1) println(hlog, "equality");
358     return s0 > s1;
359     }
360 
cellvalid(coord v)361   bool cellvalid(coord v) {
362     if(S7 == 6) return true;
363     if(S7 == 12) return (v[0] + v[1] + v[2]) % 2 == 0;
364     if(S7 == 14) return v[0] % 2 == v[1] % 2 && v[0] % 2 == v[2] % 2;
365     return false;
366     }
367 
celldistance(coord v)368   EX int celldistance(coord v) {
369     if(S7 == 6)
370       return abs(v[0]) + abs(v[1]) + abs(v[2]);
371     else {
372       for(int i=0; i<3; i++) v[i] = abs(v[i]);
373       sort(v.begin(), v.end());
374       int dist = 0;
375       if(S7 == 12) {
376         int d = v[1] - v[0]; v[1] -= d; v[2] -= d;
377         dist += d;
378         int m = min((v[2] - v[0]), v[0]);
379         dist += 2 * m;
380         v[0] -= m; v[1] -= m; v[2] -= m * 2;
381         if(v[0])
382           dist += (v[0] + v[1] + v[2]) / 2;
383         else
384           dist += v[2];
385         }
386       else {
387         dist = v[0] + (v[1] - v[0]) / 2 + (v[2] - v[0]) / 2;
388         }
389       return dist;
390       }
391     }
392 
celldistance(cell * c1,cell * c2)393   EX int celldistance(cell *c1, cell *c2) {
394     auto cm = cubemap();
395     if(GDIM == 2)
396       return dist(full_coords2(c1), full_coords2(c2));
397     return celldistance(basic_canonicalize(cm->ispacemap[c1->master] - cm->ispacemap[c2->master]));
398     }
399 
set_land(cell * c)400   EX void set_land(cell *c) {
401     if(cgflags & qPORTALSPACE) return;
402     setland(c, specialland);
403     auto m = cubemap();
404     auto co = m->ispacemap[c->master];
405 
406     int dv = 1;
407     if(geometry != gCubeTiling) dv = 2;
408 
409     int hash = 0;
410     for(int a=0; a<3; a++) hash = 1317 * hash + co[a] / 4;
411 
412     set_euland3(c, co[0]*120, co[1]*120, (co[1]+co[2]) / dv, hash);
413     }
414 
dist_relative(cell * c)415   EX int dist_relative(cell *c) {
416     auto m = cubemap();
417     auto& cc = m->camelot_center;
418     int r = roundTableRadius(NULL);
419     cell *start = m->gamestart();
420     if(!cc) {
421       cc = start;
422       while(euc::celldistance(cc, start) < r + 5)
423         cc = cc->cmove(hrand(cc->type));
424       }
425 
426     return euc::celldistance(cc, c) - r;
427     }
428 
429   /* quotient spaces */
430 
determinant(const intmatrix T)431   int determinant(const intmatrix T) {
432     int det = 0;
433     for(int i=0; i<3; i++)
434       det += T[0][i] * T[1][(i+1)%3] * T[2][(i+2)%3];
435     for(int i=0; i<3; i++)
436       det -= T[0][i] * T[1][(i+2)%3] * T[2][(i+1)%3];
437     return det;
438     }
439 
scaled_inverse(const intmatrix T)440   intmatrix scaled_inverse(const intmatrix T) {
441     intmatrix T2;
442     for(int i=0; i<3; i++)
443     for(int j=0; j<3; j++)
444       T2[j][i] = (T[(i+1)%3][(j+1)%3] * T[(i+2)%3][(j+2)%3] - T[(i+1)%3][(j+2)%3] * T[(i+2)%3][(j+1)%3]);
445     return T2;
446     }
447 
torus3(int x,int y,int z)448   EX torus_config torus3(int x, int y, int z) {
449     intmatrix T0 = euzeroall;
450     tie(T0[0][0], T0[1][1], T0[2][2]) = make_tuple(x, y, z);
451     return {T0, 0};
452     }
453 
clear_torus3()454   EX torus_config clear_torus3() {
455     return {euzeroall, 0};
456     }
457 
compute_cat(coord coo)458   coord torus_config_full::compute_cat(coord coo) {
459     coord cat = euzero;
460     auto& T2 = inverse_axes;
461     for(int i=0; i<3; i++) {
462       int val = T2[0][i] * coo[0] + T2[1][i] * coo[1] + T2[2][i] * coo[2];
463       if(i < WDIM - infinite_dims) val = gmod(val, det);
464       cat += val * main_axes[i];
465       }
466     return cat;
467     }
468 
valid_third_turn(const intmatrix & m)469   EX bool valid_third_turn(const intmatrix& m) {
470     if(m[0][2] != -m[0][0]-m[0][1]) return false;
471     if(m[1][0] != m[0][1]) return false;
472     if(m[1][1] != m[0][2]) return false;
473     if(m[1][2] != m[0][0]) return false;
474     if(m[2][0] != m[2][1]) return false;
475     if(m[2][0] != m[2][2]) return false;
476     return true;
477     }
478 
make_hantzsche_wendt(int v)479   EX torus_config make_hantzsche_wendt(int v) {
480     intmatrix im;
481     for(int i=0; i<3; i++)
482     for(int j=0; j<3; j++) im[i][j] = 0;
483 
484     for(int i=0; i<3; i++) {
485       im[i][i] = v;
486       im[i][(i+1)%3] = v;
487       }
488 
489     return {im, 32};
490     }
491 
valid_hantzsche_wendt(const intmatrix & m)492   EX bool valid_hantzsche_wendt(const intmatrix& m) {
493     return m[0][0] > 0 && m == make_hantzsche_wendt(m[0][0]).user_axes;
494     }
495 
make_third_turn(int a,int b,int c)496   EX torus_config make_third_turn(int a, int b, int c) {
497     intmatrix T0;
498     T0[0][0] = a;
499     T0[0][1] = b;
500     T0[2][0] = c;
501     T0[0][2] = -T0[0][0]-T0[0][1];
502     T0[1][0] = T0[0][1];
503     T0[1][1] = T0[0][2];
504     T0[1][2] = T0[0][0];
505     T0[2][1] = T0[2][2] = c;
506     return {T0, 8};
507     }
508 
make_quarter_turn(int a,int b,int c)509   EX torus_config make_quarter_turn(int a, int b, int c) {
510     intmatrix T0 = euzeroall;
511     T0[0][0] = a;
512     T0[0][1] = b;
513     T0[2][0] = c;
514     return {T0, 5};
515     }
516 
add(coord val)517   void torus_config_full::add(coord val) {
518     auto cat = compute_cat(val); if(hash.count(cat)) return; hash[cat] = isize(seq); seq.push_back(val);
519     }
520 
get(coord x)521   coord torus_config_full::get(coord x) {
522     auto cat = compute_cat(x);
523     auto& st = cubemap()->shifttable;
524     while(!hash.count(cat)) {
525       if(index == isize(seq)) throw hr_exception();
526       auto v = seq[index++];
527       for(auto s: st) add(v + s);
528       }
529     return seq[hash[cat]];
530     }
531 
valid_irr_torus()532   EX bool valid_irr_torus() {
533     #if CAP_IRR
534     if(!IRREGULAR) return true;
535     if(eu.twisted) return false;
536     for(int i=0; i<2; i++) {
537       auto x = eu.user_axes[i];
538       coord dm = eutester;
539       transmatrix dummy = Id;
540       bool mirr = false;
541       irr::base_config.canonicalize(x, dm, dummy, mirr);
542       auto x0 = eu.user_axes[i];
543       auto dm0 = eutester;
544       eu.canonicalize(x0, dm0, dummy, mirr);
545       if(x0 != euzero || dm0 != eutester) return false;
546       }
547     #endif
548     return true;
549     }
550 
build_torus3(eGeometry g)551   EX void build_torus3(eGeometry g) {
552 
553     int dim = ginf[g].g.gameplay_dimension;
554 
555     eu.user_axes = eu_input.user_axes;
556     if(dim == 2) eu.user_axes[2] = euzero;
557 
558     eu.optimal_axes = eu.user_axes;
559 
560     again:
561     for(int i=0; i<dim; i++) if(eu.optimal_axes[i] < euzero) eu.optimal_axes[i] = -eu.optimal_axes[i];
562     if(eu.optimal_axes[0] < eu.optimal_axes[1]) swap(eu.optimal_axes[0], eu.optimal_axes[1]);
563     if(eu.optimal_axes[1] < eu.optimal_axes[dim-1]) swap(eu.optimal_axes[1], eu.optimal_axes[dim-1]);
564     if(eu.optimal_axes[0] < eu.optimal_axes[1]) swap(eu.optimal_axes[0], eu.optimal_axes[1]);
565     for(int i=0; i<3; i++) {
566       int i1 = (i+1) % 3;
567       int i2 = (i+2) % 3;
568       for(int a=-10; a<=10; a++)
569       for(int b=-10; b<=10; b++) {
570         coord cand = eu.optimal_axes[i] + eu.optimal_axes[i1] * a + eu.optimal_axes[i2] * b;
571         if(celldistance(cand) < celldistance(eu.optimal_axes[i])) {
572           eu.optimal_axes[i] = cand;
573           goto again;
574           }
575         }
576       }
577 
578     eu.regular_axes = eu.optimal_axes;
579     eu.infinite_dims = dim;
580     for(int i=0; i<dim; i++) if(eu.optimal_axes[i] != euzero) eu.infinite_dims--;
581 
582     int attempt = 0;
583     next_attempt:
584     for(int i=dim-eu.infinite_dims; i<3; i++)
585       eu.regular_axes[i] = main_axes[(attempt+i)%3];
586 
587     eu.det = determinant(eu.regular_axes);
588     if(eu.det == 0) {
589       attempt++;
590       if(attempt == 3) {
591         println(hlog, "weird singular!\n");
592         exit(1);
593         }
594       goto next_attempt;
595       }
596 
597     if(eu.det < 0) eu.det = -eu.det;
598 
599     eu.inverse_axes = scaled_inverse(eu.regular_axes);
600     eu.reset();
601     eu.add(euzero);
602 
603     eu.twisted = eu_input.twisted;
604     if(dim == 3) {
605       auto &T0 = eu.user_axes;
606       if(valid_third_turn(eu.user_axes)) {
607         eu.twisted &= 16;
608         if(g == gRhombic3 && (T0[2][2]&1)) eu.twisted = 0;
609         if(g == gBitrunc3 && (T0[0][0]&1)) eu.twisted = 0;
610         if(g == gBitrunc3 && (T0[1][1]&1)) eu.twisted = 0;
611         }
612       else if(valid_hantzsche_wendt(eu.user_axes)) {
613         eu.twisted &= 32;
614         if(g == gBitrunc3 && (T0[0][0]&1)) eu.twisted = 0;
615         }
616       else {
617         eu.twisted &= 7;
618         if(g != gCubeTiling && ((T0[0][0]+T0[2][2]) & 1)) eu.twisted &=~ 1;
619         if(g != gCubeTiling && ((T0[1][1]+T0[2][2]) & 1)) eu.twisted &=~ 2;
620         for(int i=0; i<3; i++) for(int j=0; j<3; j++)
621           if(i != j && T0[i][j]) eu.twisted = 0;
622         if(T0[2][2] == 0) eu.twisted = 0;
623         if(T0[0][0] != T0[1][1]) eu.twisted &= 3;
624         }
625       }
626     else {
627       eu.twisted &= 8;
628       eu.twisted_vec = to_loc(eu.user_axes[1]);
629       eu.ortho_vec = to_loc(eu.user_axes[0]);
630       if(eu.twisted_vec == gp::loc{0,0}) eu.twisted = 0;
631       if(chiral(eu.twisted_vec)) eu.twisted = 0;
632       if(dscalar(eu.twisted_vec, eu.ortho_vec))
633         eu.twisted = 0;
634       }
635 
636     set_flag(ginf[g].flags, qANYQ, eu.infinite_dims < dim);
637     set_flag(ginf[g].flags, qBOUNDED, eu.infinite_dims == 0);
638     set_flag(ginf[g].flags, qSMALL, eu.infinite_dims == 0 && eu.det <= 4096);
639     bool nonori = false;
640     if(eu.twisted&1) nonori = !nonori;
641     if(eu.twisted&2) nonori = !nonori;
642     if(eu.twisted&4) nonori = !nonori;
643     if(eu.twisted&8) nonori = !nonori;
644     set_flag(ginf[g].flags, qNONORIENTABLE, nonori);
645     }
646 
build_torus3()647   EX void build_torus3() {
648     for(eGeometry g: { gEuclid, gEuclidSquare, gCubeTiling, gRhombic3, gBitrunc3})
649       build_torus3(g);
650     }
651 
swap01(transmatrix & M)652   void swap01(transmatrix& M) {
653     for(int i=0; i<4; i++) swap(M[i][0], M[i][1]);
654     }
655 
ort1()656   gp::loc ort1() { return (S3 == 3 ? gp::loc(1, -2) : gp::loc(0, 1)); }
657 
diagonal_cross(const coord & a,const coord & b)658   int diagonal_cross(const coord& a, const coord& b) {
659     return a[0]*b[1] + a[1]*b[2] + a[2]*b[0]
660          - b[0]*a[1] - b[1]*a[2] - b[2]*a[0];
661     }
662 
canonicalize(coord & x,coord & d,transmatrix & M,bool & mirr)663   void torus_config_full::canonicalize(coord& x, coord& d, transmatrix& M, bool& mirr) {
664     if(!twisted) {
665       if(infinite_dims == WDIM) return;
666       if(infinite_dims == WDIM-1) {
667         auto& o = optimal_axes;
668         while(celldistance(x + o[0]) <= celldistance(x)) x += o[0];
669         while(celldistance(x - o[0]) <  celldistance(x)) x -= o[0];
670         return;
671         }
672       x = get(x);
673       return;
674       }
675     auto& T0 = user_axes;
676     if(twisted & 32) {
677       int period = T0[0][0];
678       auto& coo = x;
679 
680       while(true) {
681         restart:
682         /* These coordinates cause the algorithm below to go in circles. We simply break if they are detected */
683         if(coo[0] >= 0 && coo[1] == period - coo[0] && coo[2] == -coo[1] && coo[0]*2 > period && coo[0] < period) return;
684         if(coo[0]*2 <= -period && coo[0] >= -period && coo[2] == period+coo[0] && coo[2] == -coo[1]) return;
685 
686         /* apply periods */
687         for(int i=0; i<3; i++) {
688           int j = (i+1) % 3;
689           int k = (i+2) % 3;
690           int v1 = coo[i] + coo[j];
691           int v2 = coo[i] - coo[j];
692           if(v1 >= period) {
693             coo[i] -= period; coo[j] -= period;
694             }
695           else if(v1 < -period) {
696             coo[i] += period; coo[j] += period;
697             }
698           else if(v2 >= period) {
699             coo[i] -= period; coo[j] += period;
700             }
701           else if(v2 < -period) {
702             coo[i] += period; coo[j] -= period;
703             }
704           else continue;
705           d[j] = -d[j]; d[k] = -d[k];
706           coo[j] = -coo[j]; coo[k] = -coo[k];
707           transmatrix S = Id;
708           S[j][j] = -1; S[k][k] = -1;
709           M = M * S;
710           goto restart;
711           }
712         return;
713         }
714       }
715     #if MAXMDIM >= 4
716     if(twisted & 16) {
717       int period = T0[2][2];
718       transmatrix RotYZX = Zero;
719       RotYZX[1][0] = 1;
720       RotYZX[2][1] = 1;
721       RotYZX[0][2] = 1;
722       RotYZX[3][3] = 1;
723       auto& coo = x;
724       while(true) {
725         auto coosum = coo[0] + coo[1] + coo[2];
726         if(coosum >= 3 * period) {
727           coo[0] -= period, coo[1] -= period, coo[2] -= period;
728           tie(d[0], d[1], d[2]) = make_tuple(d[1], d[2], d[0]);
729           tie(coo[0], coo[1], coo[2]) = make_tuple(coo[1], coo[2], coo[0]);
730           M = M * RotYZX;
731           }
732         else if(coosum < 0) {
733           coo[0] += period, coo[1] += period, coo[2] += period;
734           tie(d[0], d[1], d[2]) = make_tuple(d[2], d[0], d[1]);
735           tie(coo[0], coo[1], coo[2]) = make_tuple(coo[2], coo[0], coo[1]);
736           M = M * RotYZX * RotYZX;
737           }
738         else break;
739         }
740       if(T0[0] != euzero) {
741         while(diagonal_cross(coo, T0[1]) < 0) coo -= T0[0];
742         while(diagonal_cross(coo, T0[1]) > 0) coo += T0[0];
743         while(diagonal_cross(coo, T0[0]) > 0) coo -= T0[1];
744         while(diagonal_cross(coo, T0[0]) < 0) coo += T0[1];
745         }
746       return;
747       }
748     #endif
749     if(WDIM == 3) {
750       #if MAXMDIM >= 4
751       auto& coo = x;
752       while(coo[2] >= T0[2][2]) {
753         coo[2] -= T0[2][2];
754         if(twisted & 1) coo[0] *= -1, d[0] *= -1, M = M * MirrorX;
755         if(twisted & 2) coo[1] *= -1, d[1] *= -1, M = M * MirrorY;
756         if(twisted & 4) swap(coo[0], coo[1]), swap01(M), swap(d[0], d[1]);
757         }
758       while(coo[2] < 0) {
759         coo[2] += T0[2][2];
760         if(twisted & 4) swap(coo[0], coo[1]), swap(d[0], d[1]), swap01(M);
761         if(twisted & 1) coo[0] *= -1, d[0] *= -1, M = M * MirrorX;
762         if(twisted & 2) coo[1] *= -1, d[1] *= -1, M = M * MirrorY;
763         }
764       for(int i: {0,1})
765         if(T0[i][i]) coo[i] = gmod(coo[i], T0[i][i]);
766       return;
767       #endif
768       }
769     else {
770       gp::loc coo = to_loc(x);
771       gp::loc ort = ort1() * twisted_vec;
772       int dsc = dscalar(twisted_vec, twisted_vec);
773       gp::loc d0 (d[0], d[1]);
774       hyperpoint h = eumove(to_coord(twisted_vec)) * C0;
775       while(true) {
776         int dsx = dscalar(coo, twisted_vec);
777         if(dsx >= dsc) coo = coo - twisted_vec;
778         else if (dsx < 0) coo = coo + twisted_vec;
779         else break;
780         M = M * spintox(h) * MirrorY * rspintox(h);
781         auto s = ort * dscalar(d0, ort) * 2;
782         auto v = dscalar(ort, ort);
783         s.first /= v;
784         s.second /= v;
785         d0 = d0 - s;
786         s = ort * dscalar(coo, ort) * 2;
787         s.first /= v;
788         s.second /= v;
789         coo = coo - s;
790         mirr = !mirr;
791         }
792       if(ortho_vec != gp::loc{0,0}) {
793         int osc = dscalar(ortho_vec, ortho_vec);
794         while(true) {
795           int dsx = dscalar(coo, ortho_vec);
796           if(dsx >= osc) coo = coo - ortho_vec;
797           else if(dsx < 0) coo = coo + ortho_vec;
798           else break;
799           }
800         }
801       d[0] = d0.first; d[1] = d0.second;
802       x = to_coord(coo);
803       return;
804       }
805     }
806 
basic_canonicalize(coord x)807   coord basic_canonicalize(coord x) {
808     transmatrix M = Id;
809     auto dummy = euzero;
810     bool dm = false;
811     eu.canonicalize(x, dummy, M, dm);
812     return x;
813     }
814 
prepare_torus3()815   EX void prepare_torus3() {
816     eu_edit = eu_input;
817     }
818 
show_fundamental()819   EX void show_fundamental() {
820     initquickqueue();
821     shiftmatrix M = ggmatrix(cwt.at);
822     shiftpoint h0 = M*C0;
823     auto& T_edit = eu_edit.user_axes;
824     hyperpoint ha = M.T*(eumove(T_edit[0]) * C0 - C0) / 2;
825     hyperpoint hb = M.T*(eumove(T_edit[1]) * C0 - C0) / 2;
826     if(WDIM == 3) {
827       hyperpoint hc = M.T*(eumove(T_edit[2]) * C0 - C0) / 2;
828       for(int d:{-1,1}) for(int e:{-1,1}) {
829         queueline(h0+d*ha+e*hb-hc, h0+d*ha+e*hb+hc, 0xFFFFFFFF);
830         queueline(h0+d*hb+e*hc-ha, h0+d*hb+e*hc+ha, 0xFFFFFFFF);
831         queueline(h0+d*hc+e*ha-hb, h0+d*hc+e*ha+hb, 0xFFFFFFFF);
832         }
833       }
834     else {
835       queueline(h0+ha+hb, h0+ha-hb, 0xFFFFFFFF);
836       queueline(h0-ha+hb, h0-ha-hb, 0xFFFFFFFF);
837       queueline(h0+ha+hb, h0-ha+hb, 0xFFFFFFFF);
838       queueline(h0+ha-hb, h0-ha-hb, 0xFFFFFFFF);
839       }
840 
841     quickqueue();
842     }
843 
on_periods(gp::loc a,gp::loc b)844   intmatrix on_periods(gp::loc a, gp::loc b) {
845     intmatrix res;
846     for(int i=0; i<3; i++) for(int j=0; j<3; j++) res[i][j] = 0;
847     res[0][0] = a.first;
848     res[0][1] = a.second;
849     res[1][0] = b.first;
850     res[1][1] = b.second;
851     res[2][2] = 1;
852     return res;
853     }
854 
single_row_torus(int qty,int dy)855   torus_config single_row_torus(int qty, int dy) {
856     return { on_periods(gp::loc{dy, -1}, gp::loc{qty, 0}), 0 };
857     }
858 
regular_torus(gp::loc p)859   torus_config regular_torus(gp::loc p) {
860     return { on_periods(p, gp::loc(0,1) * p), 0 };
861     }
862 
rectangular_torus(int x,int y,bool klein)863   EX torus_config rectangular_torus(int x, int y, bool klein) {
864     if(S3 == 3) y /= 2;
865     return { on_periods(ort1() * gp::loc(y,0), gp::loc(x,0)), klein?8:0 };
866     }
867 
torus_config_option(string name,char key,torus_config tc)868   void torus_config_option(string name, char key, torus_config tc) {
869     dialog::addBoolItem(name, eu_edit.user_axes == tc.user_axes && eu_edit.twisted == tc.twisted && PURE, key);
870     dialog::add_action([tc] {
871       stop_game();
872       eu_input = eu_edit = tc;
873       set_variation(eVariation::pure);
874       start_game();
875       });
876     }
877 
878   EX int quotient_size = 2;
879 
show_torus3()880   EX void show_torus3() {
881     int dim = WDIM;
882     auto& T_edit = eu_edit.user_axes;
883     auto& twisted_edit = eu_edit.twisted;
884     cmode = sm::SIDE | sm::MAYDARK | sm::TORUSCONFIG;
885     gamescreen(1);
886     dialog::init(XLAT("Euclidean quotient spaces"));
887 
888     for(int y=0; y<dim+1; y++)
889       dialog::addBreak(100);
890 
891     dialog::addInfo(XLAT("columns specify periods"));
892     dialog::addInfo(XLAT("(vectors you need to take to get back to start)"));
893 
894     dialog::addBreak(50);
895 
896     show_fundamental();
897     if(dim == 3) {
898       bool nondiag = false;
899       for(int i=0; i<dim; i++)
900         for(int j=0; j<dim; j++)
901           if(T_edit[i][j] && i != j) nondiag = true;
902 
903       if(valid_third_turn(T_edit)) {
904         auto g = geometry;
905         if(g == gCubeTiling ||
906           (g == gRhombic3 && T_edit[2][2] % 2 == 0) ||
907           (g == gBitrunc3 && T_edit[0][0] % 2 == 0 && T_edit[1][1] % 2 == 0))
908           dialog::addBoolItem(XLAT("third-turn space"), twisted_edit & 16, 'x');
909         else
910           dialog::addBoolItem(XLAT("make it even"), twisted_edit & 16, 'x');
911         dialog::add_action([] { eu_edit.twisted ^= 16; });
912         }
913 
914       if(valid_hantzsche_wendt(T_edit)) {
915         auto g = geometry;
916         if(g == gCubeTiling || g == gRhombic3 || (g == gBitrunc3 && T_edit[0][0] % 2 == 0))
917           dialog::addBoolItem(XLAT("Hantzsche-Wendt space"), twisted_edit & 32, 'x');
918         else
919           dialog::addBoolItem(XLAT("make it even"), twisted_edit & 32, 'x');
920         dialog::add_action([] { eu_edit.twisted ^= 32; });
921         }
922 
923       if(nondiag) {
924         dialog::addInfo(XLAT("twisting implemented only for diagonal matrices"));
925         dialog::addInfo(XLAT("or for columns : (A,B,C), (B,C,A), (D,D,D) where A+B+C=0"));
926         dialog::addBreak(200);
927         }
928       else if(T_edit[dim-1][dim-1] == 0) {
929         dialog::addInfo(XLAT("nothing to twist"));
930         dialog::addInfo(XLAT("change the bottom left corner"));
931         dialog::addBreak(100);
932         }
933       else {
934         auto g = geometry;
935         if(g == gCubeTiling || (T_edit[0][0]+T_edit[2][2]) % 2 == 0)
936           dialog::addBoolItem(XLAT("flip X coordinate"), twisted_edit & 1, 'x');
937         else
938           dialog::addBoolItem(XLAT("flipping X impossible"), twisted_edit & 1, 'x');
939         dialog::add_action([] { eu_edit.twisted ^= 1; });
940 
941         if(g == gCubeTiling || (T_edit[1][1]+T_edit[2][2]) % 2 == 0)
942           dialog::addBoolItem(XLAT("flip Y coordinate"), twisted_edit & 2, 'y');
943         else
944           dialog::addBoolItem(XLAT("flipping Y impossible"), twisted_edit & 2, 'y');
945         dialog::add_action([] { eu_edit.twisted ^= 2; });
946 
947         if(T_edit[0][0] == T_edit[1][1])
948           dialog::addBoolItem(XLAT("swap X and Y"), twisted_edit & 4, 'z');
949         else
950           dialog::addBoolItem(XLAT("swapping impossible"), twisted_edit & 4, 'z');
951         dialog::add_action([] { eu_edit.twisted ^= 4; });
952         }
953       dialog::addBreak(50);
954       dialog::addItem("special manifolds", 'S');
955       dialog::add_action([] {
956         dialog::editNumber(quotient_size, 1, 12, 1, 2, "special manifold size", "");
957         dialog::extra_options = [] {
958           auto q = quotient_size;
959           torus_config_option(XLAT("third-turn space"), 'A', make_third_turn(q,0,q));
960           torus_config_option(XLAT("quarter-turn space"), 'B', make_quarter_turn(q,0,q));
961           torus_config_option(XLAT("Hantzsche-Wendt space"), 'C', make_hantzsche_wendt(q));
962           };
963         });
964       }
965     else {
966       if(T_edit[1][0] == 0 && T_edit[1][1] == 0)
967         dialog::addInfo(XLAT("change the second column for Möbius bands and Klein bottles"));
968       else if(chiral(to_loc(T_edit[1])))
969         dialog::addInfo(XLAT("second period is chiral -- cannot be mirrored"));
970       else if(dscalar(to_loc(T_edit[1]), to_loc(T_edit[0])))
971         dialog::addInfo(XLAT("periods must be orthogonal for mirroring"));
972       else {
973         dialog::addBoolItem(XLAT("mirror flip in the second period"), twisted_edit & 8, 'x');
974         dialog::add_action([] { eu_edit.twisted ^= 8; });
975         }
976 
977       dialog::addBreak(50);
978       torus_config_option(XLAT("single-cell torus"), 'A', regular_torus(gp::loc{1,0}));
979       torus_config_option(XLAT("large regular torus"), 'B', regular_torus(gp::loc{12, 0}));
980       torus_config_option(XLAT("Klein bottle"), 'C', rectangular_torus(12, 6, true));
981       torus_config_option(XLAT("cylinder"), 'D', rectangular_torus(6, 0, false));
982       torus_config_option(XLAT("Möbius band"), 'E', rectangular_torus(6, 0, true));
983       if(S3 == 3) torus_config_option(XLAT("seven-colorable torus"), 'F', regular_torus(gp::loc{1,2}));
984       if(S3 == 3) torus_config_option(XLAT("HyperRogue classic torus"), 'G', single_row_torus(381, -22));
985       torus_config_option(XLAT("no quotient"), 'H', rectangular_torus(0, 0, false));
986       }
987 
988     dialog::addBreak(50);
989     dialog::addBoolItem(XLAT("standard rotation"), eqmatrix(models::euclidean_spin, Id), 's');
990     dialog::add_action([] { rotate_view(models::euclidean_spin); });
991 
992 #if CAP_RUG
993     if(GDIM == 2) {
994       dialog::addBoolItem(XLAT("hypersian rug mode"), (rug::rugged), 'u');
995       dialog::add_action(rug::select);
996       }
997 #endif
998 
999     dialog::addBreak(50);
1000 
1001     char xch = 'p';
1002     for(eGeometry g: {gCubeTiling, gRhombic3, gBitrunc3}) {
1003       if(dim == 2) g = geometry;
1004       dialog::addItem(XLAT(ginf[g].menu_displayed_name), xch++);
1005       dialog::add_action([g] {
1006         stop_game();
1007         set_geometry(g);
1008         eu_input = eu_edit;
1009         start_game();
1010         });
1011       if(dim == 2) break;
1012       }
1013 
1014     dialog::addBreak(50);
1015     dialog::addBack();
1016     dialog::display();
1017 
1018     int i = -1;
1019     for(auto& v: dialog::items) if(v.type == dialog::diBreak) {
1020       if(i >= 0 && i < dim) {
1021         for(int j=0; j < dim; j++) {
1022           char ch = 'a' + i * 3 + j;
1023           if(displayfr(dialog::dcenter + dialog::dfspace * 4 * (j-(dim-1.)/2), v.position, 2, dialog::dfsize, its(T_edit[j][i]), 0xFFFFFF, 8))
1024             getcstat = ch;
1025           dialog::add_key_action(ch, [i, j] {
1026             auto& T_edit = eu_edit.user_axes;
1027             dialog::editNumber(T_edit[j][i], -10, +10, 1, 0, "", XLAT(
1028               "This matrix lets you play on the quotient spaces of three-dimensional. "
1029               "Euclidean space. Every column specifies a translation vector which "
1030               "takes you back to the starting point. For example, if you put "
1031               "set 2, 6, 0 on the diagonal, you get back to the starting point "
1032               "if you move 2 steps in the X direction, 6 steps in the Y direction "
1033               "(the quotient space is infinite in the Z direction).\n\n"
1034               "You can also introduce twists for diagonal matrices: after going "
1035               "the given number of steps in the Z direction, the space is also "
1036               "mirrored or rotated. (More general 'twisted' spaces are currently "
1037               "not implemented.)"
1038               )
1039               );
1040             dialog::extra_options = show_fundamental;
1041             });
1042           }
1043         }
1044       i++;
1045       }
1046     }
1047 
1048   #if CAP_COMMANDLINE
euArgs()1049   int euArgs() {
1050     using namespace arg;
1051 
1052     if(0) ;
1053     else if(argis("-t3")) {
1054       PHASEFROM(2);
1055       stop_game();
1056       auto& T0 = eu_input.user_axes;
1057       for(int i=0; i<3; i++)
1058       for(int j=0; j<3; j++) {
1059         shift(); T0[i][j] = argi();
1060         }
1061       build_torus3();
1062       }
1063     else if(argis("-t2")) {
1064       PHASEFROM(2);
1065       stop_game();
1066       auto& T0 = eu_input.user_axes;
1067       for(int i=0; i<2; i++)
1068       for(int j=0; j<2; j++) {
1069         shift(); T0[i][j] = argi();
1070         }
1071       shift(); eu_input.twisted = argi();
1072       build_torus3();
1073       }
1074     else if(argis("-twistthird")) {
1075       PHASEFROM(2);
1076       stop_game();
1077       shift(); int a = argi();
1078       shift(); int b = argi();
1079       shift(); int c = argi();
1080       eu_input = make_third_turn(a, b, c);
1081       build_torus3();
1082       }
1083     else if(argis("-twist3")) {
1084       PHASEFROM(2);
1085       stop_game();
1086       auto& T0 = eu_input.user_axes;
1087       for(int i=0; i<3; i++)
1088       for(int j=0; j<3; j++) T0[i][j] = 0;
1089 
1090       for(int i=0; i<3; i++) {
1091         shift(); T0[i][i] = argi();
1092         }
1093       shift(); eu_input.twisted = argi();
1094       build_torus3();
1095       }
1096     else if(argis("-hw")) {
1097       PHASEFROM(2);
1098       stop_game();
1099       shift();
1100       eu_input = make_hantzsche_wendt(argi());
1101       build_torus3();
1102       }
1103     else if(argis("-twisttest")) {
1104       start_game();
1105       celllister cl(cwt.at, 10000, 10000, NULL);
1106       for(cell *c: cl.lst) {
1107         heptagon *h = c->master;
1108         for(int i=0; i<S7; i++)
1109         for(int j=0; j<S7; j++)
1110         for(int k=0; k<S7; k++)
1111         for(int l=0; l<S7; l++)
1112           if(h->move(i) && c->move(k) && h->move(i)->move(j) == h->move(k)->move(l) && h->move(i)->move(j)) {
1113             transmatrix T1 = move_matrix(h, i) * move_matrix(h->move(i), j);
1114             transmatrix T2 = move_matrix(h, k) * move_matrix(h->move(k), l);
1115             if(!eqmatrix(T1, T2)) {
1116               println(hlog, c, " @ ", cubemap()->ispacemap[c->master], " : ", i, "/", j, "/", k, "/", l, " :: ", T1, " vs ", T2);
1117               exit(1);
1118               }
1119             }
1120         }
1121       }
1122 
1123     else return 1;
1124     return 0;
1125     }
1126 
1127   auto euhook = addHook(hooks_args, 100, euArgs);
1128   #endif
1129 
dscalar(gp::loc e1,gp::loc e2)1130 EX int dscalar(gp::loc e1, gp::loc e2) {
1131   return 2 * (e1.first * e2.first + e1.second*e2.second) + (S3 == 3 ? e1.first*e2.second + e2.first * e1.second : 0);
1132   }
1133 
dsquare(gp::loc e)1134 EX int dsquare(gp::loc e) { return dscalar(e, e)/2; }
1135 
dcross(gp::loc e1,gp::loc e2)1136 EX int dcross(gp::loc e1, gp::loc e2) {
1137   return e1.first * e2.second - e1.second*e2.first;
1138   }
1139 
full_coords2(cell * c)1140 EX gp::loc full_coords2(cell *c) {
1141   if(INVERSE) {
1142     cell *c1 = gp::get_mapped(c);
1143     return UIU(full_coords2(c1));
1144     }
1145   auto ans = eucmap()->ispacemap[c->master];
1146   if(S7 == 4 && BITRUNCATED) {
1147     if(c == c->master->c7) return to_loc(ans) * gp::loc(1,1);
1148     else {
1149       auto res = full_coords2(c->cmove(0)) + full_coords2(c->cmove(4));
1150       res.first /= 2;
1151       res.second /= 2;
1152       return res;
1153       }
1154     }
1155   if(BITRUNCATED)
1156     return to_loc(ans) * gp::loc(1,1) + (c == c->master->c7 ? gp::loc(0,0) : gp::eudir((c->c.spin(0)+4)%6));
1157   if(GOLDBERG) {
1158     auto li = gp::get_local_info(c);
1159     gp::loc shift(0,0);
1160     if(li.first_dir >= 0) shift = gp::eudir(li.last_dir) * li.relative;
1161     return to_loc(ans) * gp::param + shift;
1162     }
1163   return to_loc(ans);
1164   }
1165 
1166 /** this is slow, but we use it only for small p's */
at(gp::loc p)1167 EX cell* at(gp::loc p) {
1168   cellwalker cw(currentmap->gamestart());
1169   while(p.first--) cw += revstep;
1170   cw ++;
1171   while(p.second--) cw += revstep;
1172   return cw.at;
1173   }
1174 
to_coord(gp::loc p)1175 EX coord to_coord(gp::loc p) { return coord(p.first, p.second, 0); }
1176 
sdxy()1177 EX gp::loc sdxy() { return to_loc(eu.user_axes[1]) * gp::univ_param(); }
1178 
coord_display(const shiftmatrix & V,cell * c)1179 EX pair<bool, string> coord_display(const shiftmatrix& V, cell *c) {
1180   if(c != c->master->c7) return {false, ""};
1181   hyperpoint hx = eumove(main_axes[0]) * C0;
1182   hyperpoint hy = eumove(main_axes[1]) * C0;
1183   hyperpoint hz = WDIM == 2 ? C0 : eumove(main_axes[2]) * C0;
1184   hyperpoint h = kz(inverse(build_matrix(hx, hy, hz, C03)) * inverse_shift(ggmatrix(cwt.at->master->c7), V) * C0);
1185 
1186   if(WDIM == 3)
1187     return {true, fts(h[0]) + "," + fts(h[1]) + "," + fts(h[2]) };
1188   else
1189     return {true, fts(h[0]) + "," + fts(h[1]) };
1190   }
1191 
to_loc(const coord & v)1192 EX gp::loc to_loc(const coord& v) { return gp::loc(v[0], v[1]); }
1193 
get_cdata()1194 EX map<gp::loc, cdata>& get_cdata() { return eucmap()->eucdata; }
1195 
eumove(coord co)1196 EX transmatrix eumove(coord co) {
1197   const double q3 = sqrt(double(3));
1198   if(WDIM == 3) {
1199     return eupush3(co[0], co[1], co[2]);
1200     }
1201   transmatrix Mat = Id;
1202   if(a4) {
1203     Mat[0][LDIM] += co[0] * cgi.tessf;
1204     Mat[1][LDIM] += co[1] * cgi.tessf;
1205     }
1206   else {
1207     Mat[0][LDIM] += (co[0] + co[1] * .5) * cgi.tessf;
1208     Mat[1][LDIM] += co[1] * q3 /2 * cgi.tessf;
1209     }
1210   return Mat;
1211   }
1212 
eumove(gp::loc co)1213 EX transmatrix eumove(gp::loc co) { return eumove(to_coord(co)); }
1214 
chiral(gp::loc g)1215 EX bool chiral(gp::loc g) {
1216   int x = g.first;
1217   int y = g.second;
1218   if(x == 0) return false;
1219   if(y == 0) return false;
1220   if(x+y == 0) return false;
1221   if(x==y) return false;
1222   if(S3 == 3 && y == -2*x) return false;
1223   if(S3 == 3 && x == -2*y) return false;
1224   return true;
1225   }
1226 
twist_once(gp::loc coo)1227 EX void twist_once(gp::loc coo) {
1228   coo = coo - eu.twisted_vec * gp::univ_param();
1229   if(eu.twisted&8) {
1230     gp::loc ort = ort1() * eu.twisted_vec * gp::univ_param();
1231     auto s = ort * dscalar(coo, ort) * 2;
1232     auto v = dscalar(ort, ort);
1233     s.first /= v;
1234     s.second /= v;
1235     coo = coo - s;
1236     }
1237   }
1238 
dist(int sx,int sy,bool reduce IS (true))1239 EX int dist(int sx, int sy, bool reduce IS(true)) {
1240   int z0 = abs(sx);
1241   int z1 = abs(sy);
1242   if(a4 && BITRUNCATED)
1243     return (z0 == z1 && z0 > 0 && !reduce) ? z0+1: max(z0, z1);
1244   if(a4) return z0 + z1;
1245   int z2 = abs(sx+sy);
1246   return max(max(z0,z1), z2);
1247   }
1248 
dist(gp::loc a,gp::loc b)1249 EX int dist(gp::loc a, gp::loc b) {
1250   return dist(a.first-b.first, a.second-b.second, (a.first ^ a.second)&1);
1251   }
1252 
cyldist(gp::loc a,gp::loc b)1253 EX int cyldist(gp::loc a, gp::loc b) {
1254   a = to_loc(basic_canonicalize(to_coord(a)));
1255   b = to_loc(basic_canonicalize(to_coord(b)));
1256 
1257   if(!quotient) return dist(a, b);
1258 
1259   int best = 0;
1260   for(int sa=0; sa<16; sa++) {
1261     auto _a = a, _b = b;
1262     if(sa&1) twist_once(_a);
1263     if(sa&2) twist_once(_b);
1264     if(sa&4) _a = _a + eu.ortho_vec * gp::univ_param();
1265     if(sa&8) _b = _b + eu.ortho_vec * gp::univ_param();
1266     int val = dist(_a, _b);
1267     if(sa == 0 || val < best) best = val;
1268     }
1269 
1270   return best;
1271   }
1272 
generate()1273 EX void generate() {
1274 
1275   #if MAXMDIM >= 4
1276   if(fake::in()) {
1277     fake::generate();
1278     return;
1279     }
1280 
1281   auto v = euc::get_shifttable();
1282 
1283   auto& hsh = get_hsh();
1284 
1285   auto& cs = hsh.faces;
1286 
1287   cgi.loop = 4;
1288   cgi.schmid = 3;
1289 
1290   cs.clear();
1291   cs.resize(S7);
1292 
1293   if(S7 == 6) {
1294     cgi.adjcheck = 1;
1295     cgi.face = 4;
1296     for(int w=0; w<6; w++) {
1297       for(int a=0; a<4; a++) {
1298         int t[3];
1299         t[0] = (w>=3) ? -1 : 1;
1300         t[1] = among(a, 0, 3) ? -1 : 1;
1301         t[2] = among(a, 2, 3) ? -1 : 1;
1302         int x = w%3;
1303         int y = (x+2)%3;
1304         int z = (y+2)%3;
1305         cs[w].push_back(hpxy3(t[x]/2., t[y]/2., t[z]/2.));
1306         }
1307       }
1308     }
1309 
1310   if(S7 == 12) {
1311     cgi.adjcheck = sqrt(2);
1312     cgi.face = 4;
1313     for(int w=0; w<12; w++) {
1314       auto co = v[w];
1315       vector<int> valid;
1316       for(int c=0; c<3; c++) if(co[c]) valid.push_back(c);
1317       int third = 3 - valid[1] - valid[0];
1318       hyperpoint v0 = cpush0(valid[0], co[valid[0]] > 0 ? 1 : -1);
1319       hyperpoint v1 = cpush0(valid[1], co[valid[1]] > 0 ? 1 : -1);
1320       cs[w].push_back(v0);
1321       cs[w].push_back(v0/2 + v1/2 + cpush0(third, .5) - C0);
1322       cs[w].push_back(v1);
1323       cs[w].push_back(v0/2 + v1/2 + cpush0(third, -.5) - C0);
1324       }
1325     }
1326 
1327   if(S7 == 14) {
1328     cgi.adjcheck = 2;
1329     cgi.face = 4; /* the first face */
1330     auto v = euc::get_shifttable();
1331     for(int w=0; w<14; w++) {
1332       if(w%7 < 3) {
1333         int z = w>=7?-1:1;
1334         cs[w].push_back(cpush0(w%7, z) + cpush0((w%7+1)%3, 1/2.) - C0);
1335         cs[w].push_back(cpush0(w%7, z) + cpush0((w%7+2)%3, 1/2.) - C0);
1336         cs[w].push_back(cpush0(w%7, z) + cpush0((w%7+1)%3,-1/2.) - C0);
1337         cs[w].push_back(cpush0(w%7, z) + cpush0((w%7+2)%3,-1/2.) - C0);
1338         }
1339       else {
1340         auto t = v[w];
1341         ld x = t[0], y = t[1], z = t[2];
1342         for(hyperpoint h: {
1343           hpxy3(x, y/2, 0), hpxy3(x/2, y, 0), hpxy3(0, y, z/2),
1344           hpxy3(0, y/2, z), hpxy3(x/2, 0, z), hpxy3(x, 0, z/2)
1345           }) cs[w].push_back(h);
1346         }
1347       }
1348     }
1349 
1350   hsh.compute_hept();
1351   #endif
1352   }
1353 
1354 /** @brief returns true if the current geometry is based on this module
1355  *  (For example, Archimedean, kite, or fake with underlying non-Euclidean geometry returns false)
1356  */
in()1357 EX bool in() {
1358   if(fake::in()) return FPIU(in());
1359   return euclid && standard_tiling();
1360   }
1361 
in(int dim)1362 EX bool in(int dim) { return in() && WDIM == dim; }
in(int dim,int s7)1363 EX bool in(int dim, int s7) { return in(dim) && S7 == s7; }
1364 
1365 EX }
1366 
euc2_coordinates(cell * c)1367 EX gp::loc euc2_coordinates(cell *c) {
1368   if(euc::in()) return euc::full_coords2(c);
1369   hyperpoint h = calc_relative_matrix(c, currentmap->gamestart(), C0) * C0;
1370   return gp::loc(floor(h[0]), floor(h[1]));
1371   }
1372 
1373 }
1374