1 // Hyperbolic Rogue -- advanced geometry
2 // Copyright (C) 2011-2019 Zeno Rogue, see 'hyper.cpp' for details
3
4 /** \file geometry2.cpp
5 * \brief Matrices to transform between coordinates of various cells, coordinates of cell corners, etc.
6 */
7
8 #include "hyper.h"
9 namespace hr {
10
11 shiftmatrix &ggmatrix(cell *c);
12
fixelliptic(transmatrix & at)13 EX void fixelliptic(transmatrix& at) {
14 if(elliptic && at[LDIM][LDIM] < 0) {
15 for(int i=0; i<MXDIM; i++) for(int j=0; j<MXDIM; j++)
16 at[i][j] = -at[i][j];
17 }
18 }
19
fixelliptic(hyperpoint & h)20 EX void fixelliptic(hyperpoint& h) {
21 if(elliptic && h[LDIM] < 0)
22 for(int i=0; i<MXDIM; i++) h[i] = -h[i];
23 }
24
25 /** find relative_matrix via recursing the tree structure */
relative_matrix_recursive(heptagon * h2,heptagon * h1)26 EX transmatrix relative_matrix_recursive(heptagon *h2, heptagon *h1) {
27 if(gmatrix0.count(h2->c7) && gmatrix0.count(h1->c7))
28 return inverse_shift(gmatrix0[h1->c7], gmatrix0[h2->c7]);
29 transmatrix gm = Id, where = Id;
30 while(h1 != h2) {
31 for(int i=0; i<h1->type; i++) {
32 if(h1->move(i) == h2) {
33 return gm * currentmap->adj(h1, i) * where;
34 }
35 }
36 if(h1->distance > h2->distance) {
37 for(int i=0; i<h1->type; i++) if(h1->move(i) && h1->move(i)->distance < h1->distance) {
38 gm = gm * currentmap->adj(h1, i);
39 h1 = h1->move(i);
40 goto again;
41 }
42 }
43 else {
44 for(int i=0; i<h2->type; i++) if(h2->move(i) && h2->move(i)->distance < h2->distance) {
45 where = currentmap->iadj(h2, 0) * where;
46 h2 = h2->move(i);
47 goto again;
48 }
49 }
50 again: ;
51 }
52 return gm * where;
53 }
54
master_relative(cell * c,bool get_inverse)55 transmatrix hrmap_standard::master_relative(cell *c, bool get_inverse) {
56 if(0) ;
57 #if CAP_IRR
58 else if(IRREGULAR) {
59 int id = irr::cellindex[c];
60 ld alpha = 2 * M_PI / S7 * irr::periodmap[c->master].base.spin;
61 return get_inverse ? irr::cells[id].rpusher * spin(-alpha-master_to_c7_angle()): spin(alpha + master_to_c7_angle()) * irr::cells[id].pusher;
62 }
63 #endif
64 #if CAP_GP
65 else if(GOLDBERG) {
66 if(c == c->master->c7) {
67 return spin((get_inverse?-1:1) * master_to_c7_angle());
68 }
69 else {
70 auto li = gp::get_local_info(c);
71 transmatrix T = spin(master_to_c7_angle()) * cgi.gpdata->Tf[li.last_dir][li.relative.first&GOLDBERG_MASK][li.relative.second&GOLDBERG_MASK][gp::fixg6(li.total_dir)];
72 if(get_inverse) T = iso_inverse(T);
73 return T;
74 }
75 }
76 #endif
77 else if(BITRUNCATED) {
78 if(c == c->master->c7)
79 return Id;
80 return (get_inverse?cgi.invhexmove:cgi.hexmove)[c->c.spin(0)];
81 }
82 else if(WDIM == 3)
83 return Id;
84 else
85 return pispin * Id;
86 }
87
calc_relative_matrix(cell * c2,cell * c1,const hyperpoint & hint)88 EX transmatrix calc_relative_matrix(cell *c2, cell *c1, const hyperpoint& hint) {
89 return currentmap->relative_matrix(c2, c1, hint);
90 }
91
92 // target, source, direction from source to target
93
94 #if CAP_GP
95 namespace gp { extern gp::local_info draw_li; }
96 #endif
97
adj(heptagon * h,int d)98 transmatrix hrmap_standard::adj(heptagon *h, int d) {
99 if(inforder::mixed()) {
100 int t0 = h->type;
101 int t1 = h->cmove(d)->type;
102 int sp = h->c.spin(d);
103 return spin(-d * 2 * M_PI / t0) * xpush(spacedist(h->c7, d)) * spin(M_PI + 2*M_PI*sp/t1);
104 }
105 transmatrix T = cgi.heptmove[d];
106 if(h->c.mirror(d)) T = T * Mirror;
107 int sp = h->c.spin(d);
108 if(sp) T = T * spin(2*M_PI*sp/S7);
109 return T;
110 }
111
relative_matrix_via_masters(cell * c2,cell * c1,const hyperpoint & hint)112 EX transmatrix relative_matrix_via_masters(cell *c2, cell *c1, const hyperpoint& hint) {
113 heptagon *h1 = c1->master;
114 transmatrix gm = currentmap->master_relative(c1, true);
115 heptagon *h2 = c2->master;
116 transmatrix where = currentmap->master_relative(c2);
117
118 transmatrix U = currentmap->relative_matrix(h2, h1, hint);
119
120 return gm * U * where;
121 }
122
relative_matrixc(cell * c2,cell * c1,const hyperpoint & hint)123 transmatrix hrmap_standard::relative_matrixc(cell *c2, cell *c1, const hyperpoint& hint) {
124 return relative_matrix_via_masters(c2, c1, hint);
125 }
126
relative_matrixh(heptagon * h2,heptagon * h1,const hyperpoint & hint)127 transmatrix hrmap_standard::relative_matrixh(heptagon *h2, heptagon *h1, const hyperpoint& hint) {
128
129 transmatrix gm = Id, where = Id;
130 // always add to last!
131 //bool hsol = false;
132 //transmatrix sol;
133
134 set<heptagon*> visited;
135 map<ld, vector<pair<heptagon*, transmatrix>>> hbdist;
136
137 int steps = 0;
138 while(h1 != h2) {
139 steps++; if(steps > 10000) {
140 println(hlog, "not found"); return Id;
141 }
142 if(bounded) {
143 transmatrix T;
144 ld bestdist = 1e9;
145 for(int d=0; d<S7; d++) {
146 auto hm = h1->move(d);
147 if(!hm) continue;
148 transmatrix S = adj(h1, d);
149 if(hm == h2) {
150 transmatrix T1 = gm * S * where;
151 auto curdist = hdist(tC0(T1), hint);
152 if(curdist < bestdist) T = T1, bestdist = curdist;
153 }
154 if(geometry != gMinimal) for(int e=0; e<S7; e++) if(hm->move(e) == h2) {
155 transmatrix T1 = gm * S * adj(hm, e) * where;
156 auto curdist = hdist(tC0(T1), hint);
157 if(curdist < bestdist) T = T1, bestdist = curdist;
158 }
159 }
160 if(bestdist < 1e8) return T;
161 }
162 for(int d=0; d<h1->type; d++) if(h1->move(d) == h2) {
163 return gm * adj(h1, d) * where;
164 }
165 if(among(geometry, gFieldQuotient, gBring, gMacbeath)) {
166 int bestdist = 1000000, bestd = 0;
167 for(int d=0; d<S7; d++) {
168 int dist = celldistance(h1->cmove(d)->c7, h2->c7);
169 if(dist < bestdist) bestdist = dist, bestd = d;
170 }
171 gm = gm * adj(h1, bestd);
172 h1 = h1->move(bestd);
173 }
174 #if CAP_CRYSTAL
175 else if(cryst) {
176 for(int d3=0; d3<S7; d3++) {
177 auto hm = h1->cmove(d3);
178 if(visited.count(hm)) continue;
179 visited.insert(hm);
180 ld dist = crystal::space_distance(hm->c7, h2->c7);
181 hbdist[dist].emplace_back(hm, gm * adj(h1, d3));
182 }
183 auto &bestv = hbdist.begin()->second;
184 tie(h1, gm) = bestv.back();
185 bestv.pop_back();
186 if(bestv.empty()) hbdist.erase(hbdist.begin());
187 }
188 #endif
189 else if(h1->distance < h2->distance) {
190 where = iadj(h2, 0) * where;
191 h2 = h2->move(0);
192 }
193 else {
194 gm = gm * adj(h1, 0);
195 h1 = h1->move(0);
196 }
197 }
198 return gm * where;
199 }
200
ggmatrix(cell * c)201 EX shiftmatrix &ggmatrix(cell *c) {
202 shiftmatrix& t = gmatrix[c];
203 if(t[LDIM][LDIM] == 0) {
204 t.T = actual_view_transform * View * calc_relative_matrix(c, centerover, C0);
205 t.shift = 0;
206 }
207 return t;
208 }
209
210 #if HDR
211 struct horo_distance {
212 ld a, b;
213
214 void become(hyperpoint h1);
horo_distancehr::horo_distance215 horo_distance(hyperpoint h) { become(h); }
216 horo_distance(shiftpoint h1, const shiftmatrix& T);
217 bool operator < (const horo_distance z) const;
print(hstream & hs,horo_distance x)218 friend void print(hstream& hs, horo_distance x) { print(hs, "[", x.a, ":", x.b, "]"); }
219 };
220 #endif
221
become(hyperpoint h1)222 void horo_distance::become(hyperpoint h1) {
223 #if CAP_SOLV
224 if(sn::in()) {
225 a = abs(h1[2]);
226 if(asonov::in()) h1 = asonov::straighten * h1;
227 b = hypot_d(2, h1);
228 }
229 #else
230 if(0) {}
231 #endif
232 #if CAP_BT
233 else if(bt::in()) {
234 b = intval(h1, C0);
235 a = abs(bt::horo_level(h1));
236 }
237 #endif
238 else if(hybri)
239 a = 0, b = hdist(h1, C0);
240 else
241 a = 0, b = intval(h1, C0);
242 }
243
horo_distance(shiftpoint h1,const shiftmatrix & T)244 horo_distance::horo_distance(shiftpoint h1, const shiftmatrix& T) {
245 #if CAP_BT
246 if(bt::in()) become(inverse_shift(T, h1));
247 else
248 #endif
249 if(sn::in() || hybri || nil) become(inverse_shift(T, h1));
250 else
251 a = 0, b = intval(h1.h, unshift(tC0(T), h1.shift));
252 }
253
operator <(const horo_distance z) const254 bool horo_distance::operator < (const horo_distance z) const {
255 #if CAP_BT
256 if(bt::in() || sn::in()) {
257 if(a < z.a-1e-6) return true;
258 if(a > z.a+1e-6) return false;
259 }
260 #endif
261 return b < z.b - 1e-4;
262 }
263
264 template<class T, class U>
virtualRebase_cell(cell * & base,T & at,const U & check)265 void virtualRebase_cell(cell*& base, T& at, const U& check) {
266 horo_distance currz(check(at));
267 T best_at = at;
268 while(true) {
269 cell *newbase = NULL;
270 forCellIdCM(c2, i, base) {
271 transmatrix V2 = currentmap->iadj(base, i);
272 T cand_at = V2 * at;
273 horo_distance newz(check(cand_at));
274 if(newz < currz) {
275 currz = newz;
276 best_at = cand_at;
277 newbase = c2;
278 }
279 if(arb::in()) forCellIdCM(c3, j, c2) {
280 transmatrix V3 = currentmap->iadj(c2, j);
281 T cand_at3 = V3 * cand_at;
282 horo_distance newz3(check(cand_at3));
283 if(newz3 < currz) {
284 currz = newz3;
285 best_at = cand_at3;
286 newbase = c3;
287 }
288 }
289 }
290 if(!newbase) break;
291 base = newbase;
292 at = best_at;
293 }
294 #if MAXMDIM >= 4
295 if(reg3::ultra_mirror_in()) {
296 again:
297 for(auto& v: cgi.ultra_mirrors) {
298 T cand_at = v * at;
299 horo_distance newz(check(cand_at));
300 if(newz < currz) {
301 currz = newz;
302 at = cand_at;
303 goto again;
304 }
305 }
306 }
307 #endif
308 }
309
310 template<class T, class U>
virtualRebase(cell * & base,T & at,const U & check)311 void virtualRebase(cell*& base, T& at, const U& check) {
312
313 if(nil) {
314 hyperpoint h = check(at);
315 auto step = [&] (int i) {
316 at = currentmap->iadj(base, i) * at;
317 base = base->cmove(i);
318 h = check(at);
319 };
320
321 auto& nw = nilv::nilwidth;
322
323 bool ss = S7 == 6;
324
325 while(h[1] < -0.5 * nw) step(ss ? 1 : 2);
326 while(h[1] >= 0.5 * nw) step(ss ? 4 : 6);
327 while(h[0] < -0.5 * nw) step(0);
328 while(h[0] >= 0.5 * nw) step(ss ? 3 : 4);
329 while(h[2] < -0.5 * nw * nw) step(ss ? 2 : 3);
330 while(h[2] >= 0.5 * nw * nw) step(ss ? 5 : 7);
331 return;
332 }
333
334 if(prod) {
335 auto d = product_decompose(check(at)).first;
336 while(d > cgi.plevel / 2) {
337 at = currentmap->iadj(base, base->type-1) * at;
338 base = base->cmove(base->type-1); d -= cgi.plevel;
339 }
340 while(d < -cgi.plevel / 2) {
341 at = currentmap->iadj(base, base->type-2) * at;
342 base = base->cmove(base->type-2); d += cgi.plevel;
343 }
344 auto w = hybrid::get_where(base);
345 at = mscale(at, -d);
346 PIU( virtualRebase(w.first, at, check) );
347 at = mscale(at, +d);
348 base = hybrid::get_at(w.first, w.second);
349 return;
350 }
351
352 virtualRebase_cell(base, at, check);
353 }
354
virtualRebase(cell * & base,transmatrix & at)355 EX void virtualRebase(cell*& base, transmatrix& at) {
356 virtualRebase(base, at, tC0_t);
357 }
358
virtualRebase(cell * & base,hyperpoint & h)359 EX void virtualRebase(cell*& base, hyperpoint& h) {
360 // we perform fixing in check, so that it works with larger range
361 virtualRebase(base, h, [] (const hyperpoint& h) {
362 if(hyperbolic && GDIM == 2) return hpxy(h[0], h[1]);
363 if(hyperbolic && GDIM == 3) return hpxy3(h[0], h[1], h[2]);
364 return h;
365 });
366 }
367
virtualRebase(heptagon * & base,transmatrix & at)368 void hrmap_hyperbolic::virtualRebase(heptagon*& base, transmatrix& at) {
369
370 while(true) {
371
372 double currz = at[LDIM][LDIM];
373
374 heptagon *h = base;
375
376 heptagon *newbase = NULL;
377
378 transmatrix bestV {};
379
380 for(int d=0; d<S7; d++) {
381 heptspin hs(h, d, false);
382 heptspin hs2 = hs + wstep;
383 transmatrix V2 = iadj(h, d) * at;
384 double newz = V2[LDIM][LDIM];
385 if(newz < currz) {
386 currz = newz;
387 bestV = V2;
388 newbase = hs2.at;
389 }
390 }
391
392 if(newbase) {
393 base = newbase;
394 at = bestV;
395 continue;
396 }
397
398 return;
399 }
400 }
401
no_easy_spin()402 EX bool no_easy_spin() {
403 return NONSTDVAR || arcm::in() || WDIM == 3 || bt::in() || kite::in();
404 }
405
spin_angle(cell * c,int d)406 ld hrmap_standard::spin_angle(cell *c, int d) {
407 if(WDIM == 3) return SPIN_NOT_AVAILABLE;
408 ld hexshift = 0;
409 if(c == c->master->c7 && (S7 % 2 == 0) && BITRUNCATED) hexshift = cgi.hexshift + 2*M_PI/c->type;
410 else if(cgi.hexshift && c == c->master->c7) hexshift = cgi.hexshift;
411 #if CAP_IRR
412 if(IRREGULAR) {
413 auto id = irr::cellindex[c];
414 auto& vs = irr::cells[id];
415 if(d < 0 || d >= c->type) return 0;
416 auto& p = vs.jpoints[vs.neid[d]];
417 return -atan2(p[1], p[0]) - hexshift;
418 }
419 #endif
420 return M_PI - d * 2 * M_PI / c->type - hexshift;
421 }
422
423 EX transmatrix ddspin(cell *c, int d, ld bonus IS(0)) { return currentmap->spin_to(c, d, bonus); }
424 EX transmatrix iddspin(cell *c, int d, ld bonus IS(0)) { return currentmap->spin_from(c, d, bonus); }
cellgfxdist(cell * c,int d)425 EX ld cellgfxdist(cell *c, int d) { return currentmap->spacedist(c, d); }
426
427 EX transmatrix ddspin_side(cell *c, int d, ld bonus IS(0)) {
428 if(kite::in()) {
429 hyperpoint h1 = get_corner_position(c, gmod(d, c->type), 3);
430 hyperpoint h2 = get_corner_position(c, gmod(d+1, c->type) , 3);
431 hyperpoint hm = mid(h1, h2);
432 return rspintox(hm) * spin(bonus);
433 }
434 return currentmap->spin_to(c, d, bonus);
435 }
436
437 EX transmatrix iddspin_side(cell *c, int d, ld bonus IS(0)) {
438 if(kite::in()) {
439 hyperpoint h1 = get_corner_position(c, gmod(d, c->type), 3);
440 hyperpoint h2 = get_corner_position(c, gmod(d+1, c->type) , 3);
441 hyperpoint hm = mid(h1, h2);
442 return spintox(hm) * spin(bonus);
443 }
444 return currentmap->spin_from(c, d, bonus);
445 }
446
spacedist(cell * c,int i)447 double hrmap_standard::spacedist(cell *c, int i) {
448 if(NONSTDVAR || WDIM == 3) return hrmap::spacedist(c, i);
449 if(inforder::mixed()) {
450 int t0 = c->type;
451 int t1 = c->cmove(i)->type;
452 auto halfmove = [] (int i) {
453 if(i == 1) return 0.0;
454 if(i == 2) return 0.1;
455 return edge_of_triangle_with_angles(0, M_PI/i, M_PI/i);
456 };
457 ld tessf0 = halfmove(t0);
458 ld tessf1 = halfmove(t1);
459 return (tessf0 + tessf1) / 2;
460 }
461 if(!BITRUNCATED) return cgi.tessf;
462 if(c->type == S6 && (i&1)) return cgi.hexhexdist;
463 return cgi.crossf;
464 }
465
neighborId(heptagon * h1,heptagon * h2)466 int neighborId(heptagon *h1, heptagon *h2) {
467 for(int i=0; i<h1->type; i++) if(h1->move(i) == h2) return i;
468 return -1;
469 }
470
adj(cell * c,int i)471 transmatrix hrmap_standard::adj(cell *c, int i) {
472 if(GOLDBERG) {
473 transmatrix T = master_relative(c, true);
474 transmatrix U = master_relative(c->cmove(i), false);
475 heptagon *h = c->master, *h1 = c->cmove(i)->master;
476 static bool first = true;
477 if(h == h1)
478 return T * U;
479 else if(gp::do_adjm) {
480 if(gp::gp_adj.count(make_pair(c,i))) {
481 return T * gp::get_adj(c,i) * U;
482 }
483 if(first) { first = false; println(hlog, "no gp_adj"); }
484 }
485 else for(int i=0; i<h->type; i++) if(h->move(i) == h1)
486 return T * adj(h, i) * U;
487 if(first) {
488 first = false;
489 println(hlog, "not adjacent");
490 }
491 }
492 if(NONSTDVAR || WDIM == 3) {
493 return calc_relative_matrix(c->cmove(i), c, C0);
494 }
495 double d = cellgfxdist(c, i);
496 transmatrix T = ddspin(c, i) * xpush(d);
497 if(c->c.mirror(i)) T = T * Mirror;
498 cell *c1 = c->cmove(i);
499 T = T * iddspin(c1, c->c.spin(i), M_PI);
500 return T;
501 }
502
randd()503 EX double randd() { return (rand() + .5) / (RAND_MAX + 1.); }
504
randomPointIn(int t)505 EX hyperpoint randomPointIn(int t) {
506 if(NONSTDVAR || arcm::in() || kite::in()) {
507 // Let these geometries be less confusing.
508 // Also easier to implement ;)
509 return xspinpush0(2 * M_PI * randd(), asinh(randd() / 20));
510 }
511 while(true) {
512 hyperpoint h = xspinpush0(2*M_PI*(randd()-.5)/t, asinh(randd()));
513 double d =
514 PURE ? cgi.tessf : t == 6 ? cgi.hexhexdist : cgi.crossf;
515 if(hdist0(h) < hdist0(xpush(-d) * h))
516 return spin(2*M_PI/t * (rand() % t)) * h;
517 }
518 }
519
520 /** /brief get the coordinates of the vertex of cell c indexed with cid
521 * the two vertices c and c->move(cid) share are indexed cid and gmod(cid+1, c->type)
522 * cf=3 is the vertex itself; larger values are closer to the center
523 */
524
525 EX hyperpoint get_corner_position(cell *c, int cid, ld cf IS(3)) {
526 return currentmap->get_corner(c, cid, cf);
527 }
528
get_corner(cell * c,int cid,ld cf)529 hyperpoint hrmap_standard::get_corner(cell *c, int cid, ld cf) {
530 #if CAP_GP
531 if(GOLDBERG) return gp::get_corner_position(c, cid, cf);
532 #endif
533 #if CAP_IRR
534 if(IRREGULAR) {
535 auto& vs = irr::cells[irr::cellindex[c]];
536 return mid_at_actual(vs.vertices[cid], 3/cf);
537 }
538 #endif
539 if(PURE) {
540 return ddspin(c,cid,M_PI/S7) * xpush0(cgi.hcrossf * 3 / cf);
541 }
542 if(BITRUNCATED) {
543 if(!ishept(c))
544 return ddspin(c,cid,M_PI/S6) * xpush0(cgi.hexvdist * 3 / cf);
545 else
546 return ddspin(c,cid,M_PI/S7) * xpush0(cgi.rhexf * 3 / cf);
547 }
548 return C0;
549 }
550
551 EX bool approx_nearcorner = false;
552
553 /** /brief get the coordinates of the center of c->move(i) */
554
nearcorner(cell * c,int i)555 EX hyperpoint nearcorner(cell *c, int i) {
556 if(GOLDBERG_INV) {
557 i = gmod(i, c->type);
558 cellwalker cw(c, i);
559 cw += wstep;
560 transmatrix cwm = currentmap->adj(c, i);
561 if(elliptic && cwm[2][2] < 0) cwm = centralsym * cwm;
562 return cwm * C0;
563 }
564 #if CAP_IRR
565 if(IRREGULAR) {
566 auto& vs = irr::cells[irr::cellindex[c]];
567 hyperpoint nc = vs.jpoints[vs.neid[i]];
568 return mid_at(C0, nc, .94);
569 }
570 #endif
571 #if CAP_ARCM
572 if(arcm::in()) {
573 if(PURE) {
574 auto &ac = arcm::current;
575 auto& t = ac.get_triangle(c->master, i-1);
576 int id = arcm::id_of(c->master);
577 int id1 = ac.get_adj(ac.get_adj(c->master, i-1), -2).first;
578 return xspinpush0(-t.first - M_PI / c->type, ac.inradius[id/2] + ac.inradius[id1/2] + (ac.real_faces == 0 ? 2 * M_PI / (ac.N == 2 ? 2.1 : ac.N) : 0));
579 }
580 if(BITRUNCATED) {
581 auto &ac = arcm::current;
582 auto& t = ac.get_triangle(c->master, i);
583 return xspinpush0(-t.first, t.second);
584 }
585 if(DUAL) {
586 auto &ac = arcm::current;
587 auto& t = ac.get_triangle(c->master, i * 2);
588 return xspinpush0(-t.first, t.second);
589 }
590 }
591 #endif
592 #if CAP_BT
593 if(geometry == gBinary4) {
594 ld yx = log(2) / 2;
595 ld yy = yx;
596 hyperpoint neis[5];
597 neis[0] = bt::get_horopoint(2*yy, -0.5);
598 neis[1] = bt::get_horopoint(2*yy, +0.5);
599 neis[2] = bt::get_horopoint(0, 1);
600 neis[3] = bt::get_horopoint(-2*yy, c->master->zebraval ? -0.25 : +0.25);
601 neis[4] = bt::get_horopoint(0, -1);
602 return neis[i];
603 }
604 if(geometry == gTernary) {
605 ld yx = log(3) / 2;
606 ld yy = yx;
607 hyperpoint neis[6];
608 neis[0] = bt::get_horopoint(2*yy, -1);
609 neis[1] = bt::get_horopoint(2*yy, +0);
610 neis[2] = bt::get_horopoint(2*yy, +1);
611 neis[3] = bt::get_horopoint(0, 1);
612 neis[4] = bt::get_horopoint(-2*yy, c->master->zebraval / 3.);
613 neis[5] = bt::get_horopoint(0, -1);
614 return neis[i];
615 }
616 if(kite::in()) {
617 if(approx_nearcorner)
618 return currentmap->get_corner(c, i, 3) + currentmap->get_corner(c, i+1, 3) - C0;
619 else
620 return calc_relative_matrix(c->cmove(i), c, C0) * C0;
621 }
622 if(bt::in()) {
623 if(WDIM == 3) {
624 println(hlog, "nearcorner called");
625 return Hypc;
626 }
627 ld yx = log(2) / 2;
628 ld yy = yx;
629 // ld xx = 1 / sqrt(2)/2;
630 hyperpoint neis[7];
631 neis[0] = bt::get_horopoint(0, 1);
632 neis[1] = bt::get_horopoint(yy*2, 1);
633 neis[2] = bt::get_horopoint(yy*2, 0);
634 neis[3] = bt::get_horopoint(yy*2, -1);
635 neis[4] = bt::get_horopoint(0, -1);
636 if(c->type == 7)
637 neis[5] = bt::get_horopoint(-yy*2, -.5),
638 neis[6] = bt::get_horopoint(-yy*2, +.5);
639 else
640 neis[5] = bt::get_horopoint(-yy*2, 0);
641 return neis[i];
642 }
643 #endif
644 double d = cellgfxdist(c, i);
645 return ddspin(c, i) * xpush0(d);
646 }
647
648 /** /brief get the coordinates of the another vertex of c->move(i)
649 * this is useful for tessellation remapping TODO COMMENT
650 */
651
farcorner(cell * c,int i,int which)652 EX hyperpoint farcorner(cell *c, int i, int which) {
653 #if CAP_GP
654 if(GOLDBERG_INV) {
655 cellwalker cw(c, i);
656 cw += wstep;
657 if(!cw.mirrored) cw += (which?-1:2);
658 else cw += (which?2:-1);
659 transmatrix cwm = currentmap->adj(c, i);
660 if(gp::variation_for(gp::param) == eVariation::goldberg) {
661 auto li1 = gp::get_local_info(cw.at);
662 return cwm * get_corner_position(li1, cw.spin);
663 }
664 else {
665 return cwm * get_corner_position(cw.at, cw.spin, 3);
666 }
667 }
668 #endif
669 #if CAP_IRR
670 if(IRREGULAR) {
671 auto& vs = irr::cells[irr::cellindex[c]];
672 int neid = vs.neid[i];
673 int spin = vs.spin[i];
674 auto &vs2 = irr::cells[neid];
675 int cor2 = isize(vs2.vertices);
676 transmatrix rel = vs.rpusher * vs.relmatrices[vs2.owner] * vs2.pusher;
677
678 if(which == 0) return rel * vs2.vertices[(spin+2)%cor2];
679 if(which == 1) return rel * vs2.vertices[(spin+cor2-1)%cor2];
680 }
681 #endif
682 #if CAP_BT
683 if(bt::in() || kite::in())
684 return nearcorner(c, (i+which) % c->type); // lazy
685 #endif
686 #if CAP_ARCM
687 if(arcm::in()) {
688 if(PURE) {
689 auto &ac = arcm::current;
690 auto& t = ac.get_triangle(c->master, i-1);
691 int id = arcm::id_of(c->master);
692 auto id1 = ac.get_adj(ac.get_adj(c->master, i-1), -2).first;
693 int n1 = isize(ac.adjacent[id1]);
694 return spin(-t.first - M_PI / c->type) * xpush(ac.inradius[id/2] + ac.inradius[id1/2]) * xspinpush0(M_PI + M_PI/n1*(which?3:-3), ac.circumradius[id1/2]);
695 }
696 if(BITRUNCATED || DUAL) {
697 int mul = DUALMUL;
698 auto &ac = arcm::current;
699 auto adj = ac.get_adj(c->master, i * mul);
700 heptagon h; cell cx; cx.master = &h;
701 arcm::id_of(&h) = adj.first;
702 arcm::parent_index_of(&h) = adj.second;
703
704 auto& t1 = arcm::current.get_triangle(c->master, i);
705
706 auto& t2 = arcm::current.get_triangle(adj);
707
708 return spin(-t1.first) * xpush(t1.second) * spin(M_PI + t2.first) * get_corner_position(&cx, which ? -mul : 2*mul);
709 }
710 }
711 #endif
712
713 cellwalker cw(c, i);
714 cw += wstep;
715 if(!cw.mirrored) cw.spin += (which?-1:2);
716 else cw.spin += (which?2:-1);
717 return currentmap->adj(c, i) * get_corner_position(c->move(i), cw.spin);
718 }
719
midcorner(cell * c,int i,ld v)720 EX hyperpoint midcorner(cell *c, int i, ld v) {
721 auto hcor = farcorner(c, i, 0);
722 auto tcor = get_corner_position(c, i, 3);
723 return mid_at(tcor, hcor, v);
724 }
725
get_warp_corner(cell * c,int cid)726 EX hyperpoint get_warp_corner(cell *c, int cid) {
727 // midcorner(c, cid, .5) but sometimes easier versions exist
728 #if CAP_GP
729 if(GOLDBERG) return gp::get_corner_position(c, cid, 2);
730 #endif
731 #if CAP_IRR || CAP_ARCM
732 if(IRREGULAR || arcm::in()) return midcorner(c, cid, .5);
733 #endif
734 return ddspin(c,cid,M_PI/S7) * xpush0(cgi.tessf/2);
735 }
736
737 EX map<cell*, map<cell*, vector<transmatrix>>> brm_structure;
738
generate_brm(cell * c1)739 EX void generate_brm(cell *c1) {
740 set<unsigned> visited_by_matrix;
741 queue<pair<cell*, transmatrix>> q;
742 map<cell*, ld> cutoff;
743 auto& res = brm_structure[c1];
744
745 auto enqueue = [&] (cell *c, const transmatrix& T) {
746 auto b = bucketer(tC0(T));
747 if(visited_by_matrix.count(b)) return;
748 visited_by_matrix.insert(b);
749 q.emplace(c, T);
750 };
751
752 enqueue(c1, Id);
753 while(!q.empty()) {
754 cell *c2;
755 transmatrix T;
756 tie(c2,T) = q.front();
757 q.pop();
758
759 ld mindist = HUGE_VAL, maxdist = 0;
760
761 if(WDIM == 2) {
762 for(int i=0; i<c1->type; i++)
763 for(int j=0; j<c2->type; j++) {
764 ld d = hdist(get_corner_position(c1, i), T * get_corner_position(c2, j));
765 if(d < mindist) mindist = d;
766 if(d > maxdist) maxdist = d;
767 }
768 }
769 else {
770 auto& ss1 = currentmap->get_cellshape(c1);
771 auto& ss2 = currentmap->get_cellshape(c2);
772 for(auto v: ss1.vertices_only)
773 for(auto w: ss2.vertices_only) {
774 ld d = hdist(v, T*w);
775 if(d < mindist) mindist = d;
776 if(d > maxdist) maxdist = d;
777 }
778 }
779
780 auto& cu = cutoff[c2];
781 if(cu == 0 || cu > maxdist)
782 cu = maxdist;
783
784 if(mindist >= cu) continue;
785 res[c2].push_back(T);
786
787 forCellIdCM(c3, i, c2) enqueue(c3, T * currentmap->adj(c2, i));
788 }
789
790 vector<int> cts;
791 for(auto& p: res) cts.push_back(isize(p.second));
792 }
793
794 /** this function exhaustively finds the best transmatrix from (c1,h1) to (c2,h2) */
brm_get(cell * c1,hyperpoint h1,cell * c2,hyperpoint h2)795 EX const transmatrix& brm_get(cell *c1, hyperpoint h1, cell *c2, hyperpoint h2) {
796 if(!brm_structure.count(c1))
797 generate_brm(c1);
798 transmatrix *result = nullptr;
799 ld best = HUGE_VAL;
800 for(auto& t: brm_structure[c1][c2]) {
801 ld d = hdist(h1, t * h2);
802 if(d < best) best = d, result = &t;
803 }
804 return *result;
805 }
806
__anonecadbb9b0502() 807 int brm_hook = addHook(hooks_clearmemory, 0, []() {
808 brm_structure.clear();
809 });
810
exhaustive_distance_appropriate()811 EX bool exhaustive_distance_appropriate() {
812 if(euclid && (kite::in() || arcm::in() || arb::in() || quotient)) return true;
813 #if MAXMDIM >= 4
814 if(nil && quotient) return true;
815 #endif
816 #if CAP_SOLV
817 if(asonov::in() && asonov::period_xy && asonov::period_xy <= 256) return true;
818 #endif
819
820 if(bounded) return true;
821
822 return false;
823 }
824
825 #if HDR
826 struct pathgen {
827 cellwalker start;
828 cellwalker last;
829 vector<cell*> path;
830 bignum full_id_0;
831 int last_id;
832 };
833 #endif
834
generate_random_path_randomdir(cellwalker start,int length,bool for_yendor)835 EX pathgen generate_random_path_randomdir(cellwalker start, int length, bool for_yendor) {
836 start.spin = hrand(start.at->type);
837 return generate_random_path(start, length, for_yendor, false);
838 }
839
generate_random_path(cellwalker start,int length,bool for_yendor,bool randomdir)840 EX pathgen generate_random_path(cellwalker start, int length, bool for_yendor, bool randomdir) {
841 pathgen p;
842 p.start = start;
843 p.path.resize(length+1);
844 p.path[0] = start.at;
845 p.last_id = 0;
846
847 int turns = 0;
848
849 if(exhaustive_distance_appropriate()) {
850 permanent_long_distances(start.at);
851 int dist = max_saved_distance(start.at);
852 dist = min(dist, length);
853 auto at = random_in_distance(start.at, dist);
854 permanent_long_distances(at);
855 for(int a=length-1; a>=0; a--) {
856 p.path[a+1] = at;
857 vector<cell*> prev;
858 forCellCM(c2, at) if(celldistance(start.at, c2) == a) prev.push_back(c2);
859 if(isize(prev)) at = prev[hrand(isize(prev))];
860 }
861 p.path[0] = start.at;
862 p.last = p.path.back();
863 }
864
865 else if(hybri) {
866 /* I am lazy */
867 for(int i=1; i<=length; i++) p.path[i] = p.path[i-1]->cmove(p.path[i-1]->type-1);
868 p.last = p.path.back();
869 }
870
871 else {
872 int t = -1;
873 bignum full_id;
874 bool onlychild = true;
875 bool launched = false;
876
877 cellwalker ycw = start;
878 if(for_yendor) setdist(p.path[0], 7, NULL);
879
880 for(int i=0; i<length; i++) {
881
882 if(for_yendor && yendor::control(p, i, ycw)) { }
883
884 else if(bt::in()) {
885 // make it challenging
886 vector<int> ds;
887 for(int d=0; d<ycw.at->type; d++) {
888 bool increase;
889 if(sol)
890 increase = i < YDIST / 4 || i > 3 * YDIST / 4;
891 else
892 increase = i < YDIST/2;
893 if(increase) {
894 if(celldistAlt((ycw+d).cpeek()) < celldistAlt(ycw.at))
895 ds.push_back(d);
896 }
897 else {
898 if(celldistAlt((ycw+d).cpeek()) > celldistAlt(ycw.at) && (ycw+d).cpeek() != p.path[i-1])
899 ds.push_back(d);
900 }
901 }
902 if(isize(ds)) ycw += ds[hrand(isize(ds))];
903 }
904
905 else if(currentmap->strict_tree_rules()) {
906 if(for_yendor && i < arb::current.yendor_backsteps) {
907 println(hlog, i, " < ", arb::current.yendor_backsteps);
908 ycw.spin = 0;
909 }
910
911 else {
912 if(!launched) {
913 t = ycw.at->master->fieldval;
914 bignum b = expansion.get_descendants(length-i, t);
915 if(!full_id.approx_int()) goto stupid;
916 p.full_id_0 = full_id = hrand(b);
917 /* it may happen that the subtree dies out */
918 launched = true;
919 }
920
921 ycw.spin = 0;
922
923 auto& r = rulegen::treestates[t];
924 for(int ri=0; ri<isize(r.rules); ri++) {
925 int tch = r.rules[ri];
926 if(tch < 0) continue;
927 auto& sub_id = expansion.get_descendants(length-1-i, tch);
928 if(full_id < sub_id) {
929 t = tch; ycw += ri; break;
930 }
931 full_id.addmul(sub_id, -1);
932 }
933 }
934 }
935
936 else if(trees_known() && WDIM == 2) {
937 auto sdist = [start] (cell *c) { return celldistance(start.at, c); };
938 if(i == 0) {
939 t = type_in(expansion, randomdir ? start.at : start.cpeek(), sdist);
940 ycw--;
941 if(valence() == 3) ycw--;
942 bignum b = expansion.get_descendants(randomdir ? length : length-1, t);
943 p.full_id_0 = full_id = hrand(b);
944 }
945
946 #if DEBUG_YENDORGEN
947 printf("#%3d t%d %s / %s\n", i, t, full_id.get_str(100).c_str(), expansion.get_descendants(length-i, t).get_str(100).c_str());
948 for(int tch: expansion.children[t]) {
949 printf(" t%d %s\n", tch, expansion.get_descendants(length-i-1, t).get_str(100).c_str());
950 }
951 #endif
952
953 if(i == 1)
954 onlychild = true;
955 if(!onlychild) ycw++;
956 if(valence() == 3) ycw++;
957
958 onlychild = false;
959
960 for(int tch: expansion.children[t]) {
961 ycw++;
962 if(i < 2) tch = type_in(expansion, ycw.cpeek(), sdist);
963 auto& sub_id = expansion.get_descendants(length-1-i, tch);
964 if(full_id < sub_id) { t = tch; break; }
965
966 full_id.addmul(sub_id, -1);
967 onlychild = true;
968 }
969 }
970
971 else if(WDIM == 3) {
972 cell *prev = p.path[max(i-3, 0)];
973 int d = celldistance(prev, ycw.at);
974 vector<int> next;
975 forCellIdCM(c, i, ycw.at) if(celldistance(prev, c) > d) next.push_back(i);
976 if(!isize(next)) {
977 println(hlog, "error: no more cells for i=", i);
978 ycw.spin = hrand(ycw.at->type);
979 }
980 else {
981 ycw.spin = hrand_elt(next);
982 }
983 }
984
985 else {
986 stupid:
987 // stupid
988 ycw += rev;
989 // well, make it a bit more clever on bitruncated a4 grids
990 if(a4 && BITRUNCATED && S7 <= 5) {
991 if(ycw.at->type == 8 && ycw.cpeek()->type != 8)
992 ycw++;
993 if(hrand(100) < 10) {
994 if(euclid ? (turns&1) : (hrand(100) < 50))
995 ycw+=2;
996 else
997 ycw-=2;
998 turns++;
999 }
1000 }
1001 }
1002
1003 if(for_yendor) while(p.last_id < i && (p.path[p.last_id]->land == laMirror || inmirror(p.path[p.last_id]))) {
1004 p.last_id++;
1005 setdist(p.path[p.last_id], 7, nullptr);
1006 }
1007
1008 if(for_yendor && inmirror(ycw.at)) ycw = mirror::reflect(ycw);
1009 ycw += wstep;
1010 p.path[i+1] = ycw.at;
1011 }
1012 p.last = ycw + rev;
1013 }
1014 return p;
1015 }
1016
1017 }
1018