1 // Hyperbolic Rogue -- nonisotropic spaces (Solv and Nil)
2 // Copyright (C) 2011-2019 Zeno Rogue, see 'hyper.cpp' for details
3
4 /** \file nonisotropic.cpp
5 * \brief nonisotropic spaces (Solv and Nil)
6 */
7
8 #include "hyper.h"
9 namespace hr {
10
11 EX namespace nisot {
12
13 #if HDR
local_perspective_used()14 inline bool local_perspective_used() { return nonisotropic || prod; }
15 #endif
16
17 EX bool geodesic_movement = true;
18
19 EX transmatrix translate(hyperpoint h, ld co IS(1)) {
20 if(sl2 || sphere)
21 return co > 0 ? stretch::translate(h) : stretch::itranslate(h);
22 transmatrix T = Id;
23 for(int i=0; i<GDIM; i++) T[i][LDIM] = h[i];
24 if(sol && nih) {
25 T[0][0] = pow(2, -h[2]);
26 T[1][1] = pow(3, h[2]);
27 }
28 else if(sol) {
29 T[0][0] = exp(-h[2]);
30 T[1][1] = exp(+h[2]);
31 }
32 else if(nih) {
33 T[0][0] = pow(2, h[2]);
34 T[1][1] = pow(3, h[2]);
35 }
36 if(nil)
37 T[2][1] = h[0];
38 if(co < 0) return iso_inverse(T);
39 return T;
40 }
41
42 EX }
43
44 #if !CAP_SOLV
45 EX namespace sn {
in()46 EX always_false in() { return always_false(); }
47 EX }
48 #endif
49
50 #if CAP_SOLV
51 EX namespace sn {
52
in()53 EX bool in() { return cgclass == gcSolNIH; }
54
geom()55 EX eGeometry geom() {
56 if(asonov::in()) return gSol;
57 else return geometry;
58 }
59
60 #if HDR
61 typedef array<float, 3> compressed_point;
62
decompress(compressed_point p)63 inline hyperpoint decompress(compressed_point p) { return point3(p[0], p[1], p[2]); }
compress(hyperpoint h)64 inline compressed_point compress(hyperpoint h) { return make_array<float>(h[0], h[1], h[2]); }
65
66 struct tabled_inverses {
67 int PRECX, PRECY, PRECZ;
68 vector<compressed_point> tab;
69 string fname;
70 bool loaded;
71
72 void load();
73 hyperpoint get(ld ix, ld iy, ld iz, bool lazy);
74
get_inthr::sn::tabled_inverses75 compressed_point& get_int(int ix, int iy, int iz) { return tab[(iz*PRECY+iy)*PRECX+ix]; }
76
77 GLuint texture_id;
78 bool toload;
79
80 GLuint get_texture_id();
81
tabled_inverseshr::sn::tabled_inverses82 tabled_inverses(string s) : fname(s), texture_id(0), toload(true) {}
83 };
84 #endif
85
load()86 void tabled_inverses::load() {
87 if(loaded) return;
88 FILE *f = fopen(fname.c_str(), "rb");
89 if(!f) f = fopen((rsrcdir + fname).c_str(), "rb");
90 if(!f) { addMessage(XLAT("geodesic table missing")); pmodel = mdPerspective; return; }
91 ignore(fread(&PRECX, 4, 1, f));
92 ignore(fread(&PRECY, 4, 1, f));
93 ignore(fread(&PRECZ, 4, 1, f));
94 tab.resize(PRECX * PRECY * PRECZ);
95 ignore(fread(&tab[0], sizeof(compressed_point) * PRECX * PRECY * PRECZ, 1, f));
96 fclose(f);
97 loaded = true;
98 }
99
get(ld ix,ld iy,ld iz,bool lazy)100 hyperpoint tabled_inverses::get(ld ix, ld iy, ld iz, bool lazy) {
101 ix *= PRECX-1;
102 iy *= PRECY-1;
103 iz *= PRECZ-1;
104
105 hyperpoint res;
106
107 if(lazy) {
108 return decompress(get_int(int(ix+.5), int(iy+.5), int(iz+.5)));
109 }
110
111 else {
112
113 if(ix >= PRECX-1) ix = PRECX-2;
114 if(iy >= PRECX-1) iy = PRECX-2;
115 if(iz >= PRECZ-1) iz = PRECZ-2;
116
117 int ax = ix, bx = ax+1;
118 int ay = iy, by = ay+1;
119 int az = iz, bz = az+1;
120
121 #define S0(x,y,z) get_int(x, y, z)[t]
122 #define S1(x,y) (S0(x,y,az) * (bz-iz) + S0(x,y,bz) * (iz-az))
123 #define S2(x) (S1(x,ay) * (by-iy) + S1(x,by) * (iy-ay))
124
125 for(int t=0; t<3; t++)
126 res[t] = S2(ax) * (bx-ix) + S2(bx) * (ix-ax);
127
128 res[3] = 0;
129 }
130
131 return res;
132 }
133
get_texture_id()134 GLuint tabled_inverses::get_texture_id() {
135 if(!toload) return texture_id;
136
137 load();
138 if(!loaded) return 0;
139
140 println(hlog, "installing table");
141 toload = false;
142
143 if(texture_id == 0) glGenTextures(1, &texture_id);
144
145 glBindTexture( GL_TEXTURE_3D, texture_id);
146
147 glTexParameteri(GL_TEXTURE_3D, GL_TEXTURE_MIN_FILTER, GL_LINEAR);
148 glTexParameteri(GL_TEXTURE_3D, GL_TEXTURE_MAG_FILTER, GL_LINEAR);
149
150 glTexParameteri(GL_TEXTURE_3D, GL_TEXTURE_WRAP_S, GL_CLAMP_TO_EDGE);
151 glTexParameteri(GL_TEXTURE_3D, GL_TEXTURE_WRAP_T, GL_CLAMP_TO_EDGE);
152 glTexParameteri(GL_TEXTURE_3D, GL_TEXTURE_WRAP_R, GL_CLAMP_TO_EDGE);
153
154 auto xbuffer = new glvertex[PRECZ*PRECY*PRECX];
155
156 for(int z=0; z<PRECZ*PRECY*PRECX; z++) {
157 auto& t = tab[z];
158 xbuffer[z] = glhr::makevertex(t[0], t[1], t[2]);
159 }
160
161 #if !ISWEB
162 glTexImage3D(GL_TEXTURE_3D, 0, 34836 /*GL_RGBA32F*/, PRECX, PRECX, PRECZ, 0, GL_RGBA, GL_FLOAT, xbuffer);
163 #else
164 // glTexStorage3D(GL_TEXTURE_3D, 1, 34836 /*GL_RGBA32F*/, PRECX, PRECX, PRECZ);
165 // glTexSubImage3D(GL_TEXTURE_3D, 0, 0, 0, 0, PRECX, PRECY, PRECZ, GL_RGBA, GL_FLOAT, xbuffer);
166 #endif
167 delete[] xbuffer;
168 return texture_id;
169 }
170
x_to_ix(ld u)171 EX ld x_to_ix(ld u) {
172 if(u == 0.) return 0.;
173 ld diag = u*u/2.;
174
175 ld x = diag;
176 ld y = u;
177 ld z = diag+1.;
178
179 x /= (1.+z);
180 y /= (1.+z);
181
182 return 0.5 - atan((0.5-x) / y) / M_PI;
183 }
184
ix_to_x(ld ix)185 EX ld ix_to_x(ld ix) {
186 ld minx = 0;
187 while(x_to_ix(minx) <= ix) minx++;
188 ld maxx = minx; minx--;
189 for(int it=0; it<20; it++) {
190 ld x = (minx + maxx) / 2;
191 if(x_to_ix(x) < ix) minx = x;
192 else maxx = x;
193 }
194 return minx;
195 }
196
z_to_iz(ld z)197 EX ld z_to_iz(ld z) {
198 z = sinh(z) / (1 + cosh(z));
199 if(nih) z = (z+1) / 2;
200 return z;
201 }
202
iz_to_z(ld iz)203 EX ld iz_to_z(ld iz) {
204 ld minz = 0;
205 while(z_to_iz(minz) <= iz) minz++;
206 while(z_to_iz(minz) > iz) minz--;
207 ld maxz = minz + 1;
208 for(int it=0; it<20; it++) {
209 ld z = (minz + maxz) / 2;
210 if(z_to_iz(z) < iz) minz = z;
211 else maxz = z;
212 }
213 return (minz+maxz) / 2;
214 }
215
azeq_to_table(hyperpoint x)216 EX hyperpoint azeq_to_table(hyperpoint x) {
217 // azimuthal equidistant to Poincare
218 ld r = hypot_d(3, x);
219 if(r == 0) return point3(0,0,0);
220 ld make_r = sinh(r) / (1 + cosh(r));
221 ld d = make_r / r;
222 return x * d;
223 }
224
table_to_azeq(hyperpoint x)225 EX hyperpoint table_to_azeq(hyperpoint x) {
226 // Poincare to azimuthal equidistant
227 ld hr = sqhypot_d(3, x);
228 if(hr < 1e-5) return x * 2;
229 if(hr >= 1) return x * 60;
230 ld hz = (1 + hr) / (1 - hr);
231 ld d = (hz+1) * acosh(hz) / sinh(acosh(hz));
232 return x * d;
233 }
234
235
236 struct hrmap_solnih : hrmap {
237 hrmap *binary_map;
238 hrmap *ternary_map; /* nih only */
239 map<pair<heptagon*, heptagon*>, heptagon*> at;
240 map<heptagon*, pair<heptagon*, heptagon*>> coords;
241
242 heptagon *origin;
243
getOriginhr::sn::hrmap_solnih244 heptagon *getOrigin() override { return origin; }
245
get_athr::sn::hrmap_solnih246 heptagon *get_at(heptagon *x, heptagon *y) {
247 auto& h = at[make_pair(x, y)];
248 if(h) return h;
249 h = init_heptagon(S7);
250 h->c7 = newCell(S7, h);
251 coords[h] = make_pair(x, y);
252 h->distance = x->distance;
253 h->zebraval = x->emeraldval;
254 h->emeraldval = y->emeraldval;
255 return h;
256 }
257
hrmap_solnihhr::sn::hrmap_solnih258 hrmap_solnih() {
259
260 heptagon *alt;
261 heptagon *alt3;
262
263 if(true) {
264 dynamicval<eGeometry> g(geometry, gBinary4);
265 alt = init_heptagon(S7);
266 alt->s = hsOrigin;
267 alt->alt = alt;
268 binary_map = bt::new_alt_map(alt);
269 }
270
271 if(nih) {
272 dynamicval<eGeometry> g(geometry, gTernary);
273 alt3 = init_heptagon(S7);
274 alt3->s = hsOrigin;
275 alt3->alt = alt3;
276 ternary_map = bt::new_alt_map(alt3);
277 }
278 else {
279 alt3 = alt;
280 ternary_map = nullptr;
281 }
282
283 origin = get_at(alt, alt3);
284 }
285
altstephr::sn::hrmap_solnih286 heptagon *altstep(heptagon *h, int d) {
287 dynamicval<eGeometry> g(geometry, gBinary4);
288 dynamicval<hrmap*> cm(currentmap, binary_map);
289 return h->cmove(d);
290 }
291
altstep3hr::sn::hrmap_solnih292 heptagon *altstep3(heptagon *h, int d) {
293 dynamicval<eGeometry> g(geometry, gTernary);
294 dynamicval<hrmap*> cm(currentmap, ternary_map);
295 return h->cmove(d);
296 }
297
create_stephr::sn::hrmap_solnih298 heptagon *create_step(heptagon *parent, int d) override {
299 auto p = coords[parent];
300 auto pf = p.first, ps = p.second;
301 auto rule = [&] (heptagon *c1, heptagon *c2, int d1) {
302 auto g = get_at(c1, c2);
303 parent->c.connect(d, g, d1, false);
304 return g;
305 };
306
307 switch(geometry){
308
309 case gSol: switch(d) {
310 case 0: // right
311 return rule(altstep(pf, 2), ps, 4);
312 case 1: // up
313 return rule(pf, altstep(ps, 2), 5);
314 case 2: // front left
315 return rule(altstep(pf, 0), altstep(ps, 3), ps->zebraval ? 7 : 6);
316 case 3: // front right
317 return rule(altstep(pf, 1), altstep(ps, 3), ps->zebraval ? 7 : 6);
318 case 4: // left
319 return rule(altstep(pf, 4), ps, 0);
320 case 5: // down
321 return rule(pf, altstep(ps, 4), 1);
322 case 6: // back down
323 return rule(altstep(pf, 3), altstep(ps, 0), pf->zebraval ? 3 : 2);
324 case 7: // back up
325 return rule(altstep(pf, 3), altstep(ps, 1), pf->zebraval ? 3 : 2);
326 default:
327 return NULL;
328 }
329
330 case gNIH: switch(d) {
331 case 0: // right
332 return rule(altstep(pf, 2), ps, 2);
333 case 1: // up
334 return rule(pf, altstep3(ps, 3), 3);
335 case 2: // left
336 return rule(altstep(pf, 4), ps, 0);
337 case 3: // down
338 return rule(pf, altstep3(ps, 5), 1);
339 case 4: // back
340 return rule(altstep(pf, 3), altstep3(ps, 4), 5 + pf->zebraval + 2 * ps->zebraval);
341 default:
342 return rule(altstep(pf, (d-5) % 2), altstep3(ps, (d-5)/2), 4);
343 }
344
345 case gSolN: switch(d) {
346 case 0: // right
347 return rule(altstep(pf, 2), ps, 2);
348 case 1: // up
349 return rule(pf, altstep3(ps, 3), 3);
350 case 2: // left
351 return rule(altstep(pf, 4), ps, 0);
352 case 3: // down
353 return rule(pf, altstep3(ps, 5), 1);
354 case 4: case 5:
355 return rule(altstep(pf, d-4), altstep3(ps, 4), ps->zebraval + 6);
356 case 6: case 7: case 8:
357 return rule(altstep(pf, 3), altstep3(ps, d-6), pf->zebraval + 4);
358 default:
359 return NULL;
360 }
361
362 default: throw hr_exception("not solnihv");
363 }
364 }
365
~hrmap_solnihhr::sn::hrmap_solnih366 ~hrmap_solnih() {
367 delete binary_map;
368 if(ternary_map) delete ternary_map;
369 for(auto& p: at) clear_heptagon(p.second);
370 }
371
adjmatrixhr::sn::hrmap_solnih372 transmatrix adjmatrix(int i, int j) {
373 switch(geometry) {
374 case gSol: {
375 ld z = log(2);
376 ld bw = vid.binary_width * z;
377 switch(i) {
378 case 0: return xpush(+bw);
379 case 1: return ypush(+bw);
380 case 2: case 3:
381 return ypush(bw*(6.5-j)) * zpush(+z) * xpush(bw*(i-2.5));
382 case 4: return xpush(-bw);
383 case 5: return ypush(-bw);
384 case 6: case 7:
385 return xpush(bw*(2.5-j)) * zpush(-z) * ypush(bw*(i-6.5));
386 default:return Id;
387 }
388 }
389 case gNIH: {
390 ld bw = vid.binary_width;
391 switch(i) {
392 case 0: return xpush(+bw);
393 case 1: return ypush(+bw);
394 case 2: return xpush(-bw);
395 case 3: return ypush(-bw);
396 case 4: return xpush(-((j-5)%2-.5)*bw) * ypush(-((j-5)/2-1)*bw) * zpush(1);
397 default:
398 return zpush(-1) * xpush(((i-5)%2-.5)*bw) * ypush(((i-5)/2-1)*bw);
399 }
400 }
401 case gSolN: {
402 ld bw = vid.binary_width;
403 switch(i) {
404 case 0: return xpush(+bw);
405 case 1: return ypush(+bw);
406 case 2: return xpush(-bw);
407 case 3: return ypush(-bw);
408 case 4:
409 case 5:
410 return ypush(bw*(7-j)) * zpush(+1) * xpush(bw*(i-4.5));
411 case 6:
412 case 7:
413 case 8:
414 return xpush(bw*(4.5-j)) * zpush(-1) * ypush(bw*(i-7));
415 default:
416 throw hr_exception("wrong dir");
417 }
418 }
419
420 default: throw hr_exception("wrong geometry");
421 }
422 }
423
adjhr::sn::hrmap_solnih424 transmatrix adj(heptagon *h, int d) override {
425 h->cmove(d); return adjmatrix(d, h->c.spin(d));
426 }
427
relative_matrixhhr::sn::hrmap_solnih428 transmatrix relative_matrixh(heptagon *h2, heptagon *h1, const hyperpoint& hint) override {
429 for(int i=0; i<h1->type; i++) if(h1->move(i) == h2) return adjmatrix(i, h1->c.spin(i));
430 if(gmatrix0.count(h2->c7) && gmatrix0.count(h1->c7))
431 return inverse_shift(gmatrix0[h1->c7], gmatrix0[h2->c7]);
432
433 transmatrix front = Id, back = Id;
434
435 int up, down;
436
437 switch(geometry) {
438 case gSol: up = 2; down = 6; break;
439 case gSolN: up = 4; down = 7; break;
440 case gNIH: up = 4; down = 4; break;
441 default: throw hr_exception("not nihsolv");
442 }
443
444 while(h1->distance > h2->distance) front = front * adj(h1, down), h1 = h1->cmove(down);
445 while(h1->distance < h2->distance) back = iadj(h2, down) * back, h2 = h2->cmove(down);
446 while(coords[h1].first != coords[h2].first) front = front * adj(h1, down), back = iadj(h2, down) * back, h1 = h1->cmove(down), h2 = h2->cmove(down);
447 while(coords[h1].second != coords[h2].second) front = front * adj(h1, up), back = iadj(h2, up) * back, h1 = h1->cmove(up), h2 = h2->cmove(up);
448 return front * back;
449 }
450 };
451
getcoord(heptagon * h)452 EX pair<heptagon*,heptagon*> getcoord(heptagon *h) {
453 return ((hrmap_solnih*)currentmap)->coords[h];
454 }
455
get_at(heptagon * h1,heptagon * h2,bool gen)456 EX heptagon *get_at(heptagon *h1, heptagon *h2, bool gen) {
457 auto m = ((hrmap_solnih*)currentmap);
458 if(!gen && !m->at.count(make_pair(h1, h2))) return nullptr;
459 return m->get_at(h1, h2);
460 }
461
462 EX string common =
463 "uniform mediump sampler3D tInvExpTable;"
464 "uniform mediump float PRECX, PRECY, PRECZ;"
465
466 "float x_to_ix(float u) {"
467 " if(u < 1e-6) return 0.;"
468 " float diag = u*u/2.;"
469
470 " float x = diag;"
471 " float y = u;"
472 " float z = diag+1.;"
473
474 " x /= (1.+z);"
475 " y /= (1.+z);"
476
477 " return 0.5 - atan((0.5-x) / y) / 3.1415926535897932384626433832795;"
478 " }"
479
480 "float z_to_iz_s(float z) {"
481 "return sinh(z) / (1. + cosh(z));"
482 "}"
483
484 "float z_to_iz_ns(float z) {"
485 "z = sinh(z) / (1. + cosh(z));"
486 "return (z+1.)/2.;"
487 "}";
488
christoffel(const hyperpoint at,const hyperpoint velocity,const hyperpoint transported)489 hyperpoint christoffel(const hyperpoint at, const hyperpoint velocity, const hyperpoint transported) {
490 const ld l2 = log(2);
491 const ld l3 = log(3);
492 switch(geom()) {
493 case gSolN:
494 return hpxyz3(
495 -(velocity[2] * transported[0] + velocity[0] * transported[2]) * l2,
496 (velocity[2] * transported[1] + velocity[1] * transported[2]) * l3,
497 velocity[0] * transported[0] * exp(2*l2*at[2]) * l2 - velocity[1] * transported[1] * exp(-2*l3*at[2]) * l3,
498 0
499 );
500 case gSol:
501 return hpxyz3(
502 -velocity[2] * transported[0] - velocity[0] * transported[2],
503 velocity[2] * transported[1] + velocity[1] * transported[2],
504 velocity[0] * transported[0] * exp(2*at[2]) - velocity[1] * transported[1] * exp(-2*at[2]),
505 0
506 );
507 case gNIH:
508 return hpxyz3(
509 (velocity[2] * transported[0] + velocity[0] * transported[2]) * l2,
510 (velocity[2] * transported[1] + velocity[1] * transported[2]) * l3,
511 -(velocity[0] * transported[0] * exp(-2*l2*at[2]) * l2 + velocity[1] * transported[1] * exp(-2*l3*at[2]) * l3),
512 0
513 );
514 default:
515 throw hr_exception("christoffel not in solnihv");
516 }
517 }
518
get_inverse_exp_symsol(hyperpoint h,flagtype flags)519 EX hyperpoint get_inverse_exp_symsol(hyperpoint h, flagtype flags) {
520 auto& s = get_tabled();
521 s.load();
522
523 ld ix = h[0] >= 0. ? sn::x_to_ix(h[0]) : sn::x_to_ix(-h[0]);
524 ld iy = h[1] >= 0. ? sn::x_to_ix(h[1]) : sn::x_to_ix(-h[1]);
525 ld iz = sn::z_to_iz(h[2]);
526
527 if(h[2] < 0.) { iz = -iz; swap(ix, iy); }
528
529 hyperpoint res = s.get(ix, iy, iz, flags & pfNO_INTERPOLATION);
530
531 if(h[2] < 0.) { swap(res[0], res[1]); res[2] = -res[2]; }
532 if(h[0] < 0.) res[0] = -res[0];
533 if(h[1] < 0.) res[1] = -res[1];
534
535 if(flags & pfNO_DISTANCE) return res;
536 return table_to_azeq(res);
537 }
538
get_inverse_exp_nsym(hyperpoint h,flagtype flags)539 EX hyperpoint get_inverse_exp_nsym(hyperpoint h, flagtype flags) {
540 auto& s = get_tabled();
541 s.load();
542
543 ld ix = h[0] >= 0. ? sn::x_to_ix(h[0]) : sn::x_to_ix(-h[0]);
544 ld iy = h[1] >= 0. ? sn::x_to_ix(h[1]) : sn::x_to_ix(-h[1]);
545 ld iz = sn::z_to_iz(h[2]);
546
547 hyperpoint res = s.get(ix, iy, iz, flags & pfNO_INTERPOLATION);
548
549 if(h[0] < 0.) res[0] = -res[0];
550 if(h[1] < 0.) res[1] = -res[1];
551
552 if(flags & pfNO_DISTANCE) return res;
553 return table_to_azeq(res);
554 }
555
556 EX string shader_symsol = sn::common +
557
558 "vec4 inverse_exp(vec4 h) {"
559
560 "float ix = h[0] >= 0. ? x_to_ix(h[0]) : x_to_ix(-h[0]);"
561 "float iy = h[1] >= 0. ? x_to_ix(h[1]) : x_to_ix(-h[1]);"
562 "float iz = z_to_iz_s(h[2]);"
563
564 "if(h[2] < 1e-6) { iz = -iz; float s = ix; ix = iy; iy = s; }"
565 "if(iz < 0.) iz = 0.;"
566
567 "vec4 res;"
568
569 "float cx = ix*(1.-1./PRECX) + .5/PRECX;"
570 "float cy = iy*(1.-1./PRECY) + .5/PRECY;"
571 "float cz = iz*(1.-1./PRECZ) + .5/PRECZ;"
572
573 // "if(ix > .5 && iy > .6 && ix < iy + .05 && iz < .2 && iz < (iy - 0.5) * 0.6)"
574 "\n#ifndef SOLV_ALL\n"
575
576 "bool ok = true;"
577
578 // hard to tell which triangles fall on the other sides
579 "if(iz < .03 && ix > .65 && iy > .65) ok = false;"
580 "if(iz < .013 && ix > .55 && iy > .55) ok = false;"
581 "if(iz < .0075 && ix > .45 && iy > .45) ok = false;"
582 "if(iz > 0.004 && ix > 0.4 && iy > 0.4 && ix < .6 && iy < .6) ok = true;"
583 "if(iz > 0.000004 && ix > 0.4 && ix < 0.7 && iy > 0.4 && iy < 0.7) ok = true;"
584 "if(iz < 0.04 && ix > 0.70 && ix < 0.8 && iy > 0.5 && iy < 0.7) ok = false;"
585 "if(iz < 0.05 && ix > .45 && iy > .75 && ix < .55 && iy < .95) ok = false;"
586 "if(iz < 0.05 && ix > .85 && iy > .45 && iy < .75) ok = false;"
587 "if(iz < 0.025 && ix > .65 && iy > .65 && ix < .8 && iy < .8) ok = false;"
588
589 "if(!ok) res = vec4(0./0.,0./0.,0./0.,1);"
590 "else "
591
592 "\n#endif\n"
593
594 "res = texture3D(tInvExpTable, vec3(cx, cy, cz));"
595
596 "if(h[2] < 1e-6) { res.xy = res.yx; res[2] = -res[2]; }"
597 "if(h[0] < 0.) res[0] = -res[0];"
598 "if(h[1] < 0.) res[1] = -res[1];"
599
600 "return res;"
601 "}";
602
603 EX string shader_nsymsol = sn::common + R"*(
604
605 vec4 inverse_exp(vec4 h) {
606
607 float ix = h[0] >= 0. ? x_to_ix(h[0]) : x_to_ix(-h[0]);
608 float iy = h[1] >= 0. ? x_to_ix(h[1]) : x_to_ix(-h[1]);
609 float iz = z_to_iz_ns(h[2]);
610
611 vec4 res;
612
613 float cx = ix*(1.-1./PRECX) + .5/PRECX;
614 float cy = iy*(1.-1./PRECY) + .5/PRECY;
615 float cz = iz*(1.-1./PRECZ) + .5/PRECZ;
616
617 if(ix > .65 && iy > .5 && iz > .45 && iz < .55)
618 res = vec4(0.,0.,0.,1.);
619 else if(ix > .55 && iy > .75 && ix < .7 && iz > .45 && iz < .55)
620 res = vec4(0.,0.,0.,1.);
621 else if(ix > .45 && iy > .75 && ix < .7 && iz > .4 && iz < .5)
622 res = vec4(0.,0.,0.,1.);
623 else if(ix > .85 && iy > .5 && iz > .55 && iz < .75)
624 res = vec4(0.,0.,0.,1.);
625 else if(ix > .7 && iy > .55 && iz > .42 && iz < .58)
626 res = vec4(0.,0.,0.,1.);
627 else if(iz > 0.45 && ix > 0.8 && iy > 0.3 && iy < 0.6)
628 res = vec4(0.,0.,0.,1.);
629 else if(iz > 0.45 && ix > 0.8 && iy > 0.3 && iy < 0.6)
630 res = vec4(0.,0.,0.,1.);
631 else if(iz > .4 && iz < .55 && ix > .7 && iy > .36 && iy < .5 && ix < .8 && ix+iy > 1.2)
632 res = vec4(0.,0.,0.,1.);
633 else res = texture3D(tInvExpTable, vec3(cx, cy, cz));
634
635 if(h[0] < 0.) res[0] = -res[0];
636 if(h[1] < 0.) res[1] = -res[1];
637
638 return res;
639 })*";
640
641 EX string shader_nsym = sn::common +
642
643 "vec4 inverse_exp(vec4 h) {"
644
645 "float ix = h[0] >= 0. ? x_to_ix(h[0]) : x_to_ix(-h[0]);"
646 "float iy = h[1] >= 0. ? x_to_ix(h[1]) : x_to_ix(-h[1]);"
647 "float iz = z_to_iz_ns(h[2]);"
648
649 "vec4 res;"
650
651 "float cx = ix*(1.-1./PRECX) + .5/PRECX;"
652 "float cy = iy*(1.-1./PRECY) + .5/PRECY;"
653 "float cz = iz*(1.-1./PRECZ) + .5/PRECZ;"
654
655 "res = texture3D(tInvExpTable, vec3(cx, cy, cz));"
656
657 "if(h[0] < 0.) res[0] = -res[0];"
658 "if(h[1] < 0.) res[1] = -res[1];"
659
660 "return res;"
661 "}";
662
663 EX ld solrange_xy = 15;
664 EX ld solrange_z = 4;
665
in_table_range(hyperpoint h)666 EX bool in_table_range(hyperpoint h) {
667 return abs(h[0]) < solrange_xy && abs(h[1]) < solrange_xy && abs(h[2]) < solrange_z;
668 }
669
670 EX tabled_inverses solt = sn::tabled_inverses("solv-geodesics.dat");
671 EX tabled_inverses niht = sn::tabled_inverses("shyp-geodesics.dat");
672 EX tabled_inverses sont = sn::tabled_inverses("ssol-geodesics.dat");
673
get_tabled()674 EX tabled_inverses& get_tabled() {
675 switch(geom()) {
676 case gSol: return solt;
677 case gNIH: return niht;
678 case gSolN: return sont;
679 default: throw hr_exception("not solnih");
680 }
681 }
682
approx_distance(heptagon * h1,heptagon * h2)683 EX int approx_distance(heptagon *h1, heptagon *h2) {
684 auto m = (sn::hrmap_solnih*) currentmap;
685 dynamicval<eGeometry> g(geometry, gBinary4);
686 dynamicval<hrmap*> cm(currentmap, m->binary_map);
687 int d1 = bt::celldistance3_approx(m->coords[h1].first, m->coords[h2].first);
688 int d2 = bt::celldistance3_approx(m->coords[h1].second, m->coords[h2].second);
689 return d1 + d2 - abs(h1->distance - h2->distance);
690 }
691
create_faces()692 EX void create_faces() {
693 if(geometry == gSol) {
694 ld zstep = -log(2) / 2;
695 ld bwh = vid.binary_width * zstep;
696 auto pt = [&] (int x, int y, int z) { return xpush(bwh*x) * ypush(bwh*y) * zpush(zstep*z) * C0; };
697 add_wall(0, {pt(-1,-1,-1), pt(-1,-1,+1), pt(-1,00,+1), pt(-1,+1,+1), pt(-1,+1,-1)});
698 add_wall(1, {pt(-1,-1,-1), pt(00,-1,-1), pt(+1,-1,-1), pt(+1,-1,+1), pt(-1,-1,+1)});
699 add_wall(2, {pt(+1,+1,-1), pt(+1,-1,-1), pt(00,-1,-1), pt(00,+1,-1)});
700 add_wall(3, {pt(00,+1,-1), pt(00,-1,-1), pt(-1,-1,-1), pt(-1,+1,-1)});
701 add_wall(4, {pt(+1,-1,-1), pt(+1,-1,+1), pt(+1,00,+1), pt(+1,+1,+1), pt(+1,+1,-1)});
702 add_wall(5, {pt(-1,+1,-1), pt(00,+1,-1), pt(+1,+1,-1), pt(+1,+1,+1), pt(-1,+1,+1)});
703 add_wall(6, {pt(-1,+1,+1), pt(+1,+1,+1), pt(+1,00,+1), pt(-1,00,+1)});
704 add_wall(7, {pt(-1,00,+1), pt(+1,00,+1), pt(+1,-1,+1), pt(-1,-1,+1)});
705 }
706
707 if(geometry == gNIH) {
708 ld zstep = .5;
709 ld bwh = vid.binary_width / 6;
710 auto pt = [&] (int x, int y, int z) { return xpush(bwh*x) * ypush(bwh*y) * zpush(zstep*z) * C0; };
711 add_wall(0, {pt(+3,-3,-1), pt(+3,-3,+1), pt(+3,+3,+1), pt(+3,+3,-1), pt(+3,+1,-1), pt(+3,-1,-1) });
712 add_wall(1, {pt(-3,+3,-1), pt(-3,+3,+1), pt(+3,+3,+1), pt(+3,+3,-1), pt(+0,+3,-1) });
713 add_wall(2, {pt(-3,-3,-1), pt(-3,-3,+1), pt(-3,+3,+1), pt(-3,+3,-1), pt(-3,+1,-1), pt(-3,-1,-1) });
714 add_wall(3, {pt(-3,-3,-1), pt(-3,-3,+1), pt(+3,-3,+1), pt(+3,-3,-1), pt(+0,-3,-1)});
715
716 add_wall(4, {pt(-3,-3,+1), pt(-3,+3,+1), pt(+3,+3,+1), pt(+3,-3,+1)});
717
718 for(int i=0; i<6; i++) {
719 int x = -3 + (i%2) * 3;
720 int y = -3 + (i/2) * 2;
721 add_wall(5+i, {pt(x,y,-1), pt(x+3,y,-1), pt(x+3,y+2,-1), pt(x,y+2,-1)});
722 }
723 }
724
725 if(geometry == gSolN) {
726 ld zstep = -.5;
727 ld bwh = vid.binary_width / 6;
728 auto pt = [&] (int x, int y, int z) { return xpush(bwh*x) * ypush(bwh*y) * zpush(zstep*z) * C0; };
729 add_wall(0, {pt(+3,-3,-1), pt(+3,-3,+1), pt(+3,-1,+1), pt(+3,+1,+1), pt(+3,+3,+1), pt(+3,+3,-1)});
730 add_wall(1, {pt(-3,+3,-1), pt(00,+3,-1), pt(+3,+3,-1), pt(+3,+3,+1), pt(-3,+3,+1)});
731 add_wall(2, {pt(-3,-3,-1), pt(-3,-3,+1), pt(-3,-1,+1), pt(-3,+1,+1), pt(-3,+3,+1), pt(-3,+3,-1)});
732 add_wall(3, {pt(-3,-3,-1), pt(00,-3,-1), pt(+3,-3,-1), pt(+3,-3,+1), pt(-3,-3,+1)});
733 add_wall(4, {pt(-3,+3,-1), pt(-3,-3,-1), pt(00,-3,-1), pt(00,+3,-1)});
734 add_wall(5, {pt(00,+3,-1), pt(00,-3,-1), pt(+3,-3,-1), pt(+3,+3,-1)});
735 add_wall(6, {pt(-3,-3,+1), pt(+3,-3,+1), pt(+3,-1,+1), pt(-3,-1,+1)});
736 add_wall(7, {pt(-3,-1,+1), pt(+3,-1,+1), pt(+3,+1,+1), pt(-3,+1,+1)});
737 add_wall(8, {pt(-3,+1,+1), pt(+3,+1,+1), pt(+3,+3,+1), pt(-3,+3,+1)});
738 }
739
740 get_hsh().compute_hept();
741 }
742 EX }
743 #endif
744
745 EX namespace nilv {
746
christoffel(const hyperpoint Position,const hyperpoint Velocity,const hyperpoint Transported)747 hyperpoint christoffel(const hyperpoint Position, const hyperpoint Velocity, const hyperpoint Transported) {
748 ld x = Position[0];
749 return point3(
750 x * Velocity[1] * Transported[1] - 0.5 * (Velocity[1] * Transported[2] + Velocity[2] * Transported[1]),
751 -.5 * x * (Velocity[1] * Transported[0] + Velocity[0] * Transported[1]) + .5 * (Velocity[2] * Transported[0] + Velocity[0] * Transported[2]),
752 -.5 * (x*x-1) * (Velocity[1] * Transported[0] + Velocity[0] * Transported[1]) + .5 * x * (Velocity[2] * Transported[0] + Velocity[0] * Transported[2])
753 );
754 }
755
formula_exp(hyperpoint v)756 EX hyperpoint formula_exp(hyperpoint v) {
757 // copying Modelling Nil-geometry in Euclidean Space with Software Presentation
758 // v[0] = c cos alpha
759 // v[1] = c sin alpha
760 // v[2] = w
761
762 if(v[0] == 0 && v[1] == 0) return point31(v[0], v[1], v[2]);
763
764 if(v[2] == 0) return point31(v[0], v[1], v[0] * v[1] / 2);
765
766 ld alpha = atan2(v[1], v[0]);
767 ld w = v[2];
768 ld c = hypot(v[0], v[1]) / v[2];
769
770 return point31(
771 2 * c * sin(w/2) * cos(w/2 + alpha),
772 2 * c * sin(w/2) * sin(w/2 + alpha),
773 w * (1 + (c*c/2) * ((1 - sin(w)/w) + (1-cos(w))/w * sin(w + 2 * alpha)))
774 );
775 }
776
get_inverse_exp(hyperpoint h,flagtype prec IS (pNORMAL))777 EX hyperpoint get_inverse_exp(hyperpoint h, flagtype prec IS(pNORMAL)) {
778 ld wmin, wmax;
779
780 ld side = h[2] - h[0] * h[1] / 2;
781
782 if(hypot_d(2, h) < 1e-6) return point3(h[0], h[1], h[2]);
783 else if(side > 1e-6) {
784 wmin = 0, wmax = 2 * M_PI;
785 }
786 else if(side < -1e-6) {
787 wmin = - 2 * M_PI, wmax = 0;
788 }
789 else return point3(h[0], h[1], 0);
790
791 ld alpha_total = h[0] ? atan(h[1] / h[0]) : M_PI/2;
792
793 ld b;
794 if(abs(h[0]) > abs(h[1]))
795 b = h[0] / 2 / cos(alpha_total);
796 else
797 b = h[1] / 2 / sin(alpha_total);
798
799 ld s = sin(2 * alpha_total);
800
801 int max_iter = (prec & pfLOW_BS_ITER) ? 5 : 20;
802
803 for(int it=0;; it++) {
804 ld w = (wmin + wmax) / 2;
805 ld z = b * b * (s + (sin(w) - w)/(cos(w) - 1)) + w;
806
807 if(it == max_iter) {
808 ld alpha = alpha_total - w/2;
809 ld c = b / sin(w/2);
810 return point3(c * w * cos(alpha), c * w * sin(alpha), w);
811 }
812 if(h[2] > z) wmin = w;
813 else wmax = w;
814 }
815 }
816
817 EX string nilshader =
818 "vec4 inverse_exp(vec4 h) {"
819 "float wmin, wmax;"
820 "float side = h[2] - h[0] * h[1] / 2.;"
821 "if(h[0]*h[0] + h[1]*h[1] < 1e-12) return vec4(h[0], h[1], h[2], 1);"
822 "if(side > 1e-6) { wmin = 0.; wmax = 2.*PI; }"
823 "else if(side < -1e-6) { wmin = -2.*PI; wmax = 0.; }"
824 "else return vec4(h[0], h[1], 0., 1.);"
825 "float at = h[0] != 0. ? atan(h[1] / h[0]) : PI/2.;"
826 "float b = abs(h[0]) > abs(h[1]) ? h[0] / 2. / cos(at) : h[1] / 2. / sin(at);"
827 "float s = sin(2. * at);"
828
829 "for(int it=0; it<50; it++) {"
830 "float w = (wmin + wmax) / 2.;"
831 // the formula after ':' produces visible numerical artifacts for w~0
832 "float z = b * b * (s + (abs(w) < .1 ? w/3. + w*w*w/90. + w*w*w*w*w/2520.: (sin(w) - w)/(cos(w) - 1.))) + w;"
833 "if(h[2] > z) wmin = w;"
834 "else wmax = w;"
835 "}"
836
837 "float w = (wmin + wmax) / 2.;"
838 "float alpha = at - w/2.;"
839 "float c = b / sin(w/2.);"
840 "return vec4(c*w*cos(alpha), c*w*sin(alpha), w, 1.);"
841 "}";
842
843 #if HDR
844 struct mvec : array<int, 3> {
845
mvechr::nilv::mvec846 mvec() { }
847
mvechr::nilv::mvec848 mvec(int x, int y, int z) {
849 auto& a = *this;
850 a[0] = x; a[1] = y; a[2] = z;
851 }
inversehr::nilv::mvec852 mvec inverse() {
853 auto& a = *this;
854 return mvec(-a[0], -a[1], -a[2]+a[1] * a[0]);
855 }
operator *hr::nilv::mvec856 mvec operator * (const mvec b) {
857 auto& a = *this;
858 return mvec(a[0] + b[0], a[1] + b[1], a[2] + b[2] + a[0] * b[1]);
859 }
860 };
861 #endif
862
863 static const mvec mvec_zero = mvec(0, 0, 0);
864
865 EX ld nilwidth = 1;
866
mvec_to_point(mvec m)867 hyperpoint mvec_to_point(mvec m) { return hpxy3(m[0] * nilwidth, m[1] * nilwidth, m[2] * nilwidth * nilwidth); }
868
869 #if HDR
870 struct nilstructure {
871 vector<mvec> movevectors;
872 vector<vector<hyperpoint>> facevertices;
873 };
874 #endif
875
876 nilstructure ns6 = {
877 {{ mvec(-1,0,0), mvec(0,-1,0), mvec(0,0,-1), mvec(1,0,0), mvec(0,1,0), mvec(0,0,1) }},
878
879 {{
880 { point31(-0.5,-0.5,-0.25), point31(-0.5,-0.5,0.75), point31(-0.5,0.5,0.25), point31(-0.5,0.5,-0.75), },
881 { point31(0.5,-0.5,-0.5), point31(0.5,-0.5,0.5), point31(-0.5,-0.5,0.5), point31(-0.5,-0.5,-0.5), },
882 { point31(0,0,-0.5), point31(-0.5,0.5,-0.75), point31(-0.5,-0.5,-0.25), point31(0,0,-0.5), point31(-0.5,-0.5,-0.25), point31(-0.5,-0.5,-0.5), point31(0,0,-0.5), point31(-0.5,-0.5,-0.5), point31(0.5,-0.5,-0.5), point31(0,0,-0.5), point31(0.5,-0.5,-0.5), point31(0.5,-0.5,-0.75), point31(0,0,-0.5), point31(0.5,-0.5,-0.75), point31(0.5,0.5,-0.25), point31(0,0,-0.5), point31(0.5,0.5,-0.25), point31(0.5,0.5,-0.5), point31(0,0,-0.5), point31(0.5,0.5,-0.5), point31(-0.5,0.5,-0.5), point31(0,0,-0.5), point31(-0.5,0.5,-0.5), point31(-0.5,0.5,-0.75), },
883 { point31(0.5,0.5,-0.25), point31(0.5,0.5,0.75), point31(0.5,-0.5,0.25), point31(0.5,-0.5,-0.75), },
884 { point31(-0.5,0.5,-0.5), point31(-0.5,0.5,0.5), point31(0.5,0.5,0.5), point31(0.5,0.5,-0.5), },
885 { point31(0,0,0.5), point31(-0.5,0.5,0.25), point31(-0.5,-0.5,0.75), point31(0,0,0.5), point31(-0.5,-0.5,0.75), point31(-0.5,-0.5,0.5), point31(0,0,0.5), point31(-0.5,-0.5,0.5), point31(0.5,-0.5,0.5), point31(0,0,0.5), point31(0.5,-0.5,0.5), point31(0.5,-0.5,0.25), point31(0,0,0.5), point31(0.5,-0.5,0.25), point31(0.5,0.5,0.75), point31(0,0,0.5), point31(0.5,0.5,0.75), point31(0.5,0.5,0.5), point31(0,0,0.5), point31(0.5,0.5,0.5), point31(-0.5,0.5,0.5), point31(0,0,0.5), point31(-0.5,0.5,0.5), point31(-0.5,0.5,0.25), },
886 }}
887 };
888
889 nilstructure ns8 = {
890 {{ mvec(-1,0,0), mvec(-1,0,1), mvec(0,-1,0), mvec(0,0,-1), mvec(1,0,0), mvec(1,0,-1), mvec(0,1,0), mvec(0,0,1) }},
891
892 {{
893 { point31(-0.5,-0.5,-0.25), point31(-0.5,-0.5,0.75), point31(-0.5,0.5,-0.25), },
894 { point31(-0.5,-0.5,0.75), point31(-0.5,0.5,0.75), point31(-0.5,0.5,-0.25), },
895 { point31(-0.5,-0.5,-0.25), point31(-0.5,-0.5,0.75), point31(0.5,-0.5,0.25), point31(0.5,-0.5,-0.75), },
896 { point31(-0.5,-0.5,-0.25), point31(-0.5,0.5,-0.25), point31(0.5,0.5,-0.75), point31(0.5,-0.5,-0.75), },
897 { point31(0.5,0.5,0.25), point31(0.5,-0.5,0.25), point31(0.5,-0.5,-0.75), },
898 { point31(0.5,0.5,-0.75), point31(0.5,0.5,0.25), point31(0.5,-0.5,-0.75), },
899 { point31(-0.5,0.5,0.75), point31(-0.5,0.5,-0.25), point31(0.5,0.5,-0.75), point31(0.5,0.5,0.25), },
900 { point31(-0.5,-0.5,0.75), point31(-0.5,0.5,0.75), point31(0.5,0.5,0.25), point31(0.5,-0.5,0.25), },
901 }}
902 };
903
current_ns()904 EX nilstructure& current_ns() { return S7 == 6 ? ns6 : ns8; }
905
906 EX array<int,3> nilperiod, nilperiod_edit;
907 int S7_edit;
908
adjmatrix(int i)909 EX transmatrix adjmatrix(int i) {
910 return nisot::translate(mvec_to_point(current_ns().movevectors[i]));
911 }
912
913 struct hrmap_nil : hrmap {
914 map<mvec, heptagon*> at;
915 map<heptagon*, mvec> coords;
916
getOriginhr::nilv::hrmap_nil917 heptagon *getOrigin() override { return get_at(mvec_zero); }
918
~hrmap_nilhr::nilv::hrmap_nil919 ~hrmap_nil() {
920 for(auto& p: at) clear_heptagon(p.second);
921 }
922
get_athr::nilv::hrmap_nil923 heptagon *get_at(mvec c) {
924 auto& h = at[c];
925 if(h) return h;
926 h = init_heptagon(S7);
927 h->c7 = newCell(S7, h);
928 coords[h] = c;
929 h->zebraval = c[0];
930 h->emeraldval = c[1];
931 h->fieldval = c[2];
932 return h;
933 }
934
create_stephr::nilv::hrmap_nil935 heptagon *create_step(heptagon *parent, int d) override {
936 auto p = coords[parent];
937 auto q = p * current_ns().movevectors[d];
938 for(int a=0; a<3; a++) q[a] = zgmod(q[a], nilperiod[a]);
939 auto child = get_at(q);
940 parent->c.connect(d, child, (d + S7/2) % S7, false);
941 return child;
942 }
943
adjhr::nilv::hrmap_nil944 transmatrix adj(heptagon *h, int i) override { return adjmatrix(i); }
945
relative_matrixhhr::nilv::hrmap_nil946 transmatrix relative_matrixh(heptagon *h2, heptagon *h1, const hyperpoint& hint) override {
947 for(int a=0; a<S7; a++) if(h2 == h1->move(a)) return adjmatrix(a);
948 auto p = coords[h1].inverse() * coords[h2];
949 for(int a=0; a<3; a++) p[a] = szgmod(p[a], nilperiod[a]);
950 return nisot::translate(mvec_to_point(p));
951 }
952 };
953
get_coord(heptagon * h)954 EX mvec get_coord(heptagon *h) { return ((hrmap_nil*)currentmap)->coords[h]; }
955
get_heptagon_at(mvec m)956 EX heptagon *get_heptagon_at(mvec m) { return ((hrmap_nil*)currentmap)->get_at(m); }
957
set_flags()958 EX void set_flags() {
959 int coords = 0;
960 for(int a=0; a<3; a++) if(nilperiod[a]) coords++;
961 set_flag(ginf[gNil].flags, qANYQ, coords);
962 set_flag(ginf[gNil].flags, qBOUNDED, coords == 3);
963 set_flag(ginf[gNil].flags, qSMALL, coords == 3 && nilperiod[0] * nilperiod[1] * nilperiod[2] <= 4096);
964 }
965
on_geodesic(hyperpoint s0,hyperpoint s1,ld x)966 EX hyperpoint on_geodesic(hyperpoint s0, hyperpoint s1, ld x) {
967 hyperpoint local = nisot::translate(s0, -1) * s1;
968 hyperpoint h = get_inverse_exp(local);
969 return nisot::translate(s0) * formula_exp(h * x);
970 }
971
colorize(cell * c,char whichCanvas)972 EX color_t colorize(cell *c, char whichCanvas) {
973 mvec at = ((hrmap_nil*)currentmap)->coords[c->master];
974 color_t res = 0;
975
976 auto setres = [&] (int z, color_t which) {
977 if(zgmod(at[2] - z, nilperiod[2]) == 0) res = which;
978 if(zgmod(at[2] - z-1, nilperiod[2]) == 0) res = which;
979 };
980
981 if(at[1] == 0 && at[0] >=0 && at[0] < 4)
982 setres(-at[0], gradient(0x1FF0000, 0x10000FF, 0, at[0], 4));
983 else if(at[0] == 4 && at[1] >= 0 && at[1] < 4)
984 setres(at[1]*3-4, gradient(0x10000FF, 0x100FF00, 0, at[1], 4));
985 else if(at[1] == 4 && at[0] >= 0 && at[0] <= 4)
986 setres(4+at[0], gradient(0x100FF00, 0x1FFFF00, 4, at[0], 0));
987 else if(at[0] == 0 && at[1] >= 0 && at[1] <= 4)
988 setres(at[1], gradient(0x1FFFF00, 0x1FF0000, 4, at[1], 0));
989
990 return res;
991 }
992
prepare_niltorus3()993 EX void prepare_niltorus3() {
994 nilperiod_edit = nilperiod;
995 S7_edit = ginf[gNil].sides;
996 }
997
show_niltorus3()998 EX void show_niltorus3() {
999 cmode = sm::SIDE | sm::MAYDARK;
1000 gamescreen(1);
1001 dialog::init(XLAT("Nil quotient spaces"));
1002 for(int a=0; a<3; a++) {
1003 string title = XLAT("%1 period", s0+char('X'+a));
1004 dialog::addSelItem(title, its(nilperiod_edit[a]), 'x');
1005 dialog::add_action([=] {
1006 dialog::editNumber(nilperiod_edit[a], 0, 60, 1, 0, title,
1007 XLAT("Set to 0 to make it non-periodic.")
1008 );
1009 dialog::bound_low(0);
1010 });
1011 }
1012 dialog::addSelItem(XLAT("honeycomb"), its(S7_edit), 'h');
1013 dialog::add_action([] { S7_edit = S7_edit ^ 6 ^ 8; });
1014
1015 bool ok = (!nilperiod_edit[1]) || (nilperiod_edit[2] && nilperiod_edit[1] % nilperiod_edit[2] == 0);
1016
1017 dialog::addBreak(50);
1018
1019 if(ok) {
1020 dialog::addItem(XLAT("activate"), 'a');
1021 dialog::add_action([] {
1022 stop_game();
1023 nilperiod = nilperiod_edit;
1024 ginf[gNil].sides = S7_edit;
1025 set_flags();
1026 geometry = gNil;
1027 start_game();
1028 });
1029 }
1030 else dialog::addInfo(XLAT("Y period must be divisible by Z period"));
1031
1032 dialog::addBreak(50);
1033 dialog::addBack();
1034 dialog::display();
1035 }
1036
create_faces()1037 EX void create_faces() {
1038 for(int i=0; i<S7; i++) {
1039 vector<hyperpoint> fvs = nilv::current_ns().facevertices[i];
1040 using nilv::nilwidth;
1041 for(auto& h: fvs) h[0] *= nilwidth, h[1] *= nilwidth, h[2] *= nilwidth * nilwidth;
1042 add_wall(i, fvs);
1043 }
1044 get_hsh().compute_hept();
1045 }
1046
1047 EX }
1048
in_s2xe()1049 EX bool in_s2xe() { return prod && hybrid::under_class() == gcSphere; }
in_h2xe()1050 EX bool in_h2xe() { return prod && hybrid::under_class() == gcHyperbolic; }
in_e2xe()1051 EX bool in_e2xe() { return prod && hybrid::under_class() == gcEuclid; }
1052
1053 EX namespace hybrid {
1054
1055 EX eGeometry underlying;
1056 EX geometry_information *underlying_cgip;
1057
under_class()1058 EX eGeometryClass under_class() { return ginf[hybrid::underlying].cclass; }
1059
1060 EX int csteps;
1061
1062 EX int disc_quotient = 0;
1063
1064 EX map<heptagon*, short> altmap_heights;
1065
configure(eGeometry g)1066 EX void configure(eGeometry g) {
1067 if(WDIM == 3) return;
1068 ray::reset_raycaster();
1069 check_cgi();
1070 cgi.require_basics();
1071 underlying = geometry;
1072 underlying_cgip = cgip;
1073 bool sph = sphere;
1074 auto keep = ginf[g].menu_displayed_name;
1075 ginf[g] = ginf[underlying];
1076 ginf[g].menu_displayed_name = keep;
1077 if(g == gRotSpace) {
1078 ginf[g].g = sph ? giSphere3 : giSL2;
1079 ginf[g].tiling_name = "Iso(" + ginf[g].tiling_name + ")";
1080 string& qn = ginf[g].quotient_name;
1081 if(csteps && csteps != (sph ? cgi.psl_steps*2 : 0)) {
1082 string qplus;
1083 if(csteps == cgi.psl_steps)
1084 qplus = sph ? "elliptic" : "PSL";
1085 else if(csteps == 2 * cgi.psl_steps && !sph)
1086 qplus = "SL";
1087 else qplus = its(csteps);
1088 if(qn == "none") qn = qplus;
1089 else qn = qn + "/" + qplus;
1090 }
1091 if(sph) ginf[g].flags |= qELLIPTIC;
1092 if(csteps && csteps != cgi.psl_steps && csteps != 2*cgi.psl_steps)
1093 ginf[g].flags |= qANYQ;
1094 }
1095 else {
1096 ginf[g].cclass = g == gRotSpace ? gcSL2 : gcProduct;
1097 ginf[g].g.gameplay_dimension++;
1098 ginf[g].g.graphical_dimension++;
1099 ginf[g].tiling_name += "xZ";
1100 if(csteps) ginf[g].flags |= qANYQ, ginf[g].tiling_name += its(csteps);
1101 }
1102 ginf[g].flags |= qHYBRID;
1103 }
1104
reconfigure()1105 EX void reconfigure() {
1106 if(!hybri) return;
1107 stop_game();
1108 auto g = geometry;
1109 geometry = underlying;
1110 configure(g);
1111 geometry = g;
1112 }
1113
1114 EX hrmap *pmap;
1115 EX geometry_information *pcgip;
1116 EX eGeometry actual_geometry;
1117
in_actual(const T & t)1118 template<class T> auto in_actual(const T& t) -> decltype(t()) {
1119 dynamicval<eGeometry> g(geometry, actual_geometry);
1120 dynamicval<geometry_information*> gc(cgip, pcgip);
1121 dynamicval<hrmap*> gu(currentmap, pmap);
1122 dynamicval<hrmap*> gup(pmap, NULL);
1123 return t();
1124 }
1125
1126 struct hrmap_hybrid : hrmap {
1127
1128 hrmap *underlying_map;
1129
1130 bool twisted;
1131 map<cell*, pair<cellwalker, cellwalker>> spins;
1132
1133 map<pair<cell*, int>, cell*> at;
1134 map<cell*, pair<cell*, int>> where;
1135
getOriginhr::hybrid::hrmap_hybrid1136 heptagon *getOrigin() override { return underlying_map->getOrigin(); }
1137
1138 const int SHIFT_UNKNOWN = 30000;
1139 map<cell*, vector<int>> shifts;
1140
make_shifthr::hybrid::hrmap_hybrid1141 EX vector<int>& make_shift(cell *c) {
1142 auto& res = shifts[c];
1143 if(res.empty()) res = vector<int> (c->type+1, SHIFT_UNKNOWN);
1144 return res;
1145 }
1146
get_shift_currenthr::hybrid::hrmap_hybrid1147 EX int& get_shift_current(cellwalker cw) {
1148 return make_shift(cw.at)[cw.spin];
1149 }
1150
have_shifthr::hybrid::hrmap_hybrid1151 EX bool have_shift(cellwalker cw) {
1152 return shifts.count(cw.at) && get_shift_current(cw) != SHIFT_UNKNOWN;
1153 }
1154
get_shifthr::hybrid::hrmap_hybrid1155 EX int get_shift(cellwalker cw0) {
1156 if(S3 >= OINF) return 0;
1157 auto& v = get_shift_current(cw0);
1158 if(v != SHIFT_UNKNOWN) return v;
1159
1160 vector<int> candidates;
1161
1162 for(int a: {1, -1}) {
1163 cellwalker cw = cw0;
1164 cw += wstep; cw += a;
1165 int s = 0;
1166 while(cw != cw0) {
1167 if(!have_shift(cw)) goto next;
1168 s += shifts[cw.at][cw.spin];
1169 cw += wstep;
1170 cw += a;
1171 }
1172 candidates.push_back(-a * cgi.single_step * (sphere ? -1 : 1) - s);
1173 next: ;
1174 }
1175
1176 if(candidates.size() == 2 && candidates[0] != candidates[1]) {
1177 int val = candidates[0] - candidates[1];
1178 if(disc_quotient == 0) disc_quotient = val;
1179 disc_quotient = gcd(val, disc_quotient);
1180 if(disc_quotient < 0) disc_quotient = -disc_quotient;
1181 }
1182
1183 int val = 0;
1184
1185 auto cw1 = cw0+wstep;
1186 if(1) {
1187 /* the value from PSL, helps to draw the underlying space correctly */
1188 auto ps = cgi.psl_steps;
1189 val = cw0.spin*ps / cw0.at->type - cw1.spin*ps / cw1.at->type + ps/2;
1190 }
1191 if(!candidates.empty()) val = candidates[0];
1192
1193 v = val;
1194 get_shift_current(cw1) = -val;
1195
1196 return val;
1197 }
1198
ensure_shiftshr::hybrid::hrmap_hybrid1199 EX void ensure_shifts(cell *c) {
1200 if(S3 >= OINF) return;
1201 if(!make_shift(c)[c->type]) return;
1202 forCellEx(c1, c)
1203 for(int a=0; a<c->type; a++) {
1204 cellwalker cw0(c, a);
1205 cellwalker cw = cw0;
1206 while(cw != cw0) {
1207 get_shift(cw);
1208 cw += wstep;
1209 cw += a;
1210 }
1211 }
1212 make_shift(c)[c->type] = 0;
1213 }
1214
cycle_discrepancyhr::hybrid::hrmap_hybrid1215 EX int cycle_discrepancy(cellwalker cw0) {
1216 int total = 0;
1217 auto cw = cw0;
1218 do {
1219 total += get_shift(cw);
1220 cw += wstep;
1221 cw++;
1222 }
1223 while(cw != cw0);
1224 return total + cgi.single_step * (sphere ? -1 : 1);
1225 }
1226
fix_bounded_cycleshr::hybrid::hrmap_hybrid1227 EX void fix_bounded_cycles() {
1228 if(!rotspace) return;
1229 if(!bounded) return;
1230 in_underlying([&] {
1231 cellwalker final(currentmap->gamestart(), 0);
1232 auto& ac = currentmap->allcells();
1233 for(cell *c: ac) for(int i=0; i<c->type; i++) {
1234 cellwalker cw(c, i);
1235 int cd = cycle_discrepancy(cw);
1236 if(!cd) continue;
1237 while(cw != final) {
1238 if(celldist(cw.peek()) < celldist(cw.at)) {
1239 cw += wstep;
1240 cw++;
1241 }
1242 else {
1243 get_shift_current(cw) -= cd;
1244 get_shift_current(cw+wstep) += cd;
1245 cw++;
1246 }
1247 }
1248 }
1249
1250 disc_quotient = abs(cycle_discrepancy(final));
1251
1252 if(debugflags & DF_GEOM) for(cell *c: ac) for(int i=0; i<c->type; i++) {
1253 cellwalker cw(c, i);
1254 if(cycle_discrepancy(cw)) println(hlog, cw, cycle_discrepancy(cw));
1255 }
1256 });
1257 }
1258
in_underlyinghr::hybrid::hrmap_hybrid1259 template<class T> auto in_underlying(const T& t) -> decltype(t()) {
1260 pcgip = cgip;
1261 dynamicval<hrmap*> gpm(pmap, this);
1262 dynamicval<eGeometry> gag(actual_geometry, geometry);
1263 dynamicval<eGeometry> g(geometry, underlying);
1264 dynamicval<int> gss(underlying_cgip->single_step, cgi.single_step);
1265 dynamicval<int> gsp(underlying_cgip->psl_steps, cgi.psl_steps);
1266 dynamicval<geometry_information*> gc(cgip, underlying_cgip);
1267 dynamicval<hrmap*> gu(currentmap, underlying_map);
1268 return t();
1269 }
1270
getCellhr::hybrid::hrmap_hybrid1271 cell *getCell(cell *u, int h) {
1272 if(twisted) {
1273 if(!spins.count(u))
1274 println(hlog, "link missing: ", u);
1275 else {
1276 while(h >= csteps) h -= csteps, u = spins[u].first.at;
1277 while(h < 0) h += csteps, u = spins[u].second.at;
1278 }
1279 }
1280 h = zgmod(h, csteps);
1281 cell*& c = at[make_pair(u, h)];
1282 if(!c) { c = newCell(u->type+2, u->master); where[c] = {u, h}; }
1283 return c;
1284 }
1285
gamestarthr::hybrid::hrmap_hybrid1286 cell* gamestart() override { return getCell(underlying_map->gamestart(), 0); }
1287
hrmap_hybridhr::hybrid::hrmap_hybrid1288 hrmap_hybrid() {
1289 twisted = false;
1290 disc_quotient = 0;
1291 in_underlying([this] { initcells(); underlying_map = currentmap; });
1292 for(hrmap*& m: allmaps) if(m == underlying_map) m = NULL;
1293 fix_bounded_cycles();
1294 }
1295
~hrmap_hybridhr::hybrid::hrmap_hybrid1296 ~hrmap_hybrid() {
1297 in_underlying([] { delete currentmap; });
1298 for(auto& p: at) destroy_cell(p.second);
1299 }
1300
find_cell_connectionhr::hybrid::hrmap_hybrid1301 void find_cell_connection(cell *c, int d) override {
1302 hybrid::find_cell_connection(c, d);
1303 }
1304
shvidhr::hybrid::hrmap_hybrid1305 int shvid(cell *c) override {
1306 cell *c1 = hybrid::get_where(c).first;
1307 return PIU( hr::shvid(c1) );
1308 }
1309
__anon0a4104a40c02null1310 virtual transmatrix spin_to(cell *c, int d, ld bonus) override { if(d >= c->type-2) return Id; c = get_where(c).first; return in_underlying([&] { return currentmap->spin_to(c, d, bonus); }); }
__anon0a4104a40d02null1311 virtual transmatrix spin_from(cell *c, int d, ld bonus) override { if(d >= c->type-2) return Id; c = get_where(c).first; return in_underlying([&] { return currentmap->spin_from(c, d, bonus); }); }
1312
get_cellshapehr::hybrid::hrmap_hybrid1313 subcellshape& get_cellshape(cell *c) override {
1314 int id = full_shvid(c);
1315 return generate_subcellshape_if_needed(c, id);
1316 }
1317 };
1318
hmap()1319 hrmap_hybrid* hmap() { return (hrmap_hybrid*) currentmap; }
1320
get_at(cell * base,int level)1321 EX cell *get_at(cell *base, int level) {
1322 return hmap()->getCell(base, level);
1323 }
1324
get_where(cell * c)1325 EX pair<cell*, int> get_where(cell *c) { return hmap()->where[c]; }
1326
find_cell_connection(cell * c,int d)1327 EX void find_cell_connection(cell *c, int d) {
1328 auto m = hmap();
1329 if(d >= c->type - 2) {
1330 int s = cgi.single_step;
1331 int lev = m->where[c].second + (d == c->type-1 ? s : -s);
1332 cell *c1 = get_at(m->where[c].first, lev);
1333 c->c.connect(d, c1, c1->type - 3 + c->type - d, false);
1334 }
1335 else {
1336 auto cu = m->where[c].first;
1337 auto cu1 = m->in_underlying([&] { return cu->cmove(d); });
1338 int d1 = cu->c.spin(d);
1339 int s = 0;
1340 if(geometry == gRotSpace) {
1341 auto cm = (hrmap_hybrid*)currentmap;
1342 m->in_underlying([&] { cm->ensure_shifts(cu); });
1343 s = ((hrmap_hybrid*)currentmap)->get_shift(cellwalker(cu, d));
1344 }
1345 cell *c1 = get_at(cu1, m->where[c].second + s);
1346 c->c.connect(d, c1, d1, cu->c.mirror(d));
1347 }
1348 }
1349
get_umap()1350 EX hrmap* get_umap() { if(!dynamic_cast<hrmap_hybrid*>(currentmap)) return nullptr; else return ((hrmap_hybrid*)currentmap)->underlying_map; }
1351
1352 #if HDR
in_underlying_geometry(const T & f)1353 template<class T> auto in_underlying_geometry(const T& f) -> decltype(f()) {
1354 if(!hybri) return f();
1355 dynamicval<eGeometry> g(geometry, underlying);
1356 dynamicval<eGeometry> gag(actual_geometry, geometry);
1357 dynamicval<int> gss(underlying_cgip->single_step, cgi.single_step);
1358 dynamicval<int> gsp(underlying_cgip->psl_steps, cgi.psl_steps);
1359 dynamicval<geometry_information*> gc(cgip, underlying_cgip);
1360 dynamicval<hrmap*> gpm(pmap, currentmap);
1361 dynamicval<hrmap*> gm(currentmap, get_umap());
1362 return f();
1363 }
1364
1365 #define PIU(x) hr::hybrid::in_underlying_geometry([&] { return (x); })
1366 #endif
1367
get_corner(cell * c,int i,int next,ld z)1368 EX hyperpoint get_corner(cell *c, int i, int next, ld z) {
1369 ld lev = cgi.plevel * z / 2;
1370 if(WDIM == 2) {
1371 ld zz = lerp(cgi.FLOOR, cgi.WALL, (1+z) / 2);
1372 hyperpoint h = zshift(get_corner_position(c, i+next), zz);
1373 return h;
1374 }
1375 if(prod) {
1376 dynamicval<eGeometry> g(geometry, hybrid::underlying);
1377 dynamicval<geometry_information*> gc(cgip, hybrid::underlying_cgip);
1378 dynamicval<hrmap*> gm(currentmap, ((hrmap_hybrid*)currentmap)->underlying_map);
1379 return mscale(get_corner_position(c, i+next), exp(lev));
1380 }
1381 else {
1382 #if MAXMDIM >= 4
1383 ld tf, he, alpha;
1384 in_underlying_geometry([&] {
1385 hyperpoint h1 = get_corner_position(c, i);
1386 hyperpoint h2 = get_corner_position(c, i+1);
1387 hyperpoint hm = mid(h1, h2);
1388 tf = hdist0(hm)/2;
1389 he = hdist(hm, h2)/2;
1390 alpha = atan2(hm[1], hm[0]);
1391 });
1392 return spin(alpha) * rots::uxpush(tf) * rots::uypush(next?he:-he) * rots::uzpush(lev) * C0;
1393 #else
1394 throw hr_exception();
1395 #endif
1396 }
1397 }
1398
__anon0a4104a41102() 1399 auto clear_samples = addHook(hooks_clearmemory, 40, [] () {
1400 for(auto& c: cgis) for(auto& v: c.second.walloffsets)
1401 v.second = nullptr;
1402 altmap_heights.clear();
1403 });
1404
gen_sample_list()1405 EX vector<pair<int, cell*>> gen_sample_list() {
1406 if(!hybri && WDIM != 2 && PURE)
1407 return {make_pair(0, centerover), make_pair(centerover->type, nullptr)};
1408 vector<pair<int, cell*>> result;
1409 for(auto& v: cgi.walloffsets) if(v.first >= 0) result.push_back(v);
1410 sort(result.begin(), result.end());
1411 result.emplace_back(isize(cgi.wallstart), nullptr);
1412 return result;
1413 }
1414
1415 vector<cell*> to_link;
1416
will_link(cell * c)1417 EX void will_link(cell *c) { if(pmap && ((hrmap_hybrid*) pmap)->twisted) to_link.push_back(c); }
1418
1419 EX bool in_link = false;
1420
link()1421 EX void link() {
1422 if(in_link) return;
1423 dynamicval<bool> b(in_link, true);
1424 auto pm = (hrmap_hybrid*) pmap;
1425 if(!pm) return;
1426 auto& ss = pm->spins;
1427 int success = -1;
1428 while(success) {
1429 vector<cell*> xlink = std::move(to_link);
1430 success = 0;
1431 for(cell *c: xlink) {
1432 bool success_here = ss.count(c);
1433 if(!success_here) forCellIdEx(c2, i, c) if(ss.count(c2)) {
1434 ss[c].first = ss[c2].first + c->c.spin(i) + wstep - i;
1435 ss[c].second = ss[c2].second + c->c.spin(i) + wstep - i;
1436 success++;
1437 success_here = true;
1438 break;
1439 }
1440 if(!success_here) to_link.push_back(c);
1441 }
1442 }
1443 }
1444
celldistance(cell * c1,cell * c2)1445 EX int celldistance(cell *c1, cell *c2) {
1446 if(sl2) {
1447 auto w1 = hybrid::get_where(c1), w2 = hybrid::get_where(c2);
1448 return PIU (hr::celldistance(w1.first, w2.first));
1449 }
1450 else if(csteps == 0) {
1451 auto w1 = hybrid::get_where(c1), w2 = hybrid::get_where(c2);
1452 return PIU (hr::celldistance(w1.first, w2.first)) + abs(w1.second - w2.second);
1453 }
1454 else {
1455 int s = 0;
1456 int a = 999999, b = -999999;
1457 auto c = c1;
1458 do {
1459 auto w1 = hybrid::get_where(c), w2 = hybrid::get_where(c2);
1460 if(w1.second == w2.second) {
1461 int d = PIU(hr::celldistance(w1.first, w2.first));
1462 a = min(s+d, a);
1463 b = max(s-d, b);
1464 }
1465 c = c->cmove(c1->type-1); s++;
1466 }
1467 while(c != c1);
1468 return min(a, s-b);
1469 }
1470 }
1471
configure_period()1472 EX void configure_period() {
1473 static int s;
1474 s = csteps / cgi.single_step;
1475 string str = "";
1476 if(rotspace)
1477 str = XLAT(
1478 "If the 2D underlying manifold is bounded, the period should be a divisor of the 'rotation space' "
1479 "value (PSL(2,R)) times the Euler characteristics of the underlying manifold. "
1480 "For unbounded underlying manifold, any value should work (theoretically, "
1481 "the current implementation in HyperRogue is not perfect).\n\n"
1482 "We won't stop you from trying illegal numbers, but they won't work correctly.");
1483 dialog::editNumber(s, 0, 16, 1, 0, XLAT("%1 period", "Z"), str);
1484 dialog::bound_low(0);
1485 auto set_s = [] (int s1, bool ret) {
1486 return [s1,ret] {
1487 if(ret) popScreen();
1488 if(csteps == s1) return;
1489 stop_game();
1490 csteps = s1 * cgi.single_step;
1491 hybrid::reconfigure();
1492 start_game();
1493 };
1494 };
1495 dialog::extra_options = [=] () {
1496 if(rotspace) {
1497 int e_steps = cgi.psl_steps / gcd(cgi.single_step, cgi.psl_steps);
1498 bool ubounded = PIU(bounded);
1499 dialog::addSelItem( sphere ? XLAT("elliptic") : XLAT("PSL(2,R)"), its(e_steps), 'P');
1500 dialog::add_action(set_s(e_steps, true));
1501 dialog::addSelItem( sphere ? XLAT("sphere") : XLAT("SL(2,R)"), its(2*e_steps), 'P');
1502 dialog::add_action(set_s(2*e_steps, true));
1503 if(sl2 && !ubounded) {
1504 dialog::addSelItem( XLAT("universal cover"), its(0), 'P');
1505 dialog::add_action(set_s(0, true));
1506 }
1507 dialog::addSelItem(ubounded ? XLAT("maximum") : XLAT("works correctly so far"), its(disc_quotient), 'Q');
1508 dialog::add_action(set_s(disc_quotient, true));
1509 }
1510 else {
1511 dialog::addSelItem( XLAT("non-periodic"), its(0), 'N');
1512 dialog::add_action(set_s(0, true));
1513 }
1514 dialog::reaction_final = set_s(s, false);
1515 };
1516 }
1517
1518 EX }
1519
1520 EX namespace product {
1521
1522 int z0;
1523
1524 struct hrmap_product : hybrid::hrmap_hybrid {
relative_matrixchr::product::hrmap_product1525 transmatrix relative_matrixc(cell *c2, cell *c1, const hyperpoint& hint) override {
1526 return in_underlying([&] { return calc_relative_matrix(where[c2].first, where[c1].first, hint); }) * mscale(Id, cgi.plevel * szgmod(where[c2].second - where[c1].second, hybrid::csteps));
1527 }
1528
adjhr::product::hrmap_product1529 transmatrix adj(cell *c, int i) override {
1530 if(twisted && i == c->type-1 && where[c].second == hybrid::csteps-1) {
1531 auto b = spins[where[c].first].first;
1532 transmatrix T = mscale(Id, cgi.plevel);
1533 T = T * spin(2 * M_PI * b.spin / b.at->type);
1534 if(b.mirrored) T = T * Mirror;
1535 return T;
1536 }
1537 if(twisted && i == c->type-2 && where[c].second == 0) {
1538 auto b = spins[where[c].first].second;
1539 transmatrix T = mscale(Id, -cgi.plevel);
1540 T = T * spin(2 * M_PI * b.spin / b.at->type);
1541 if(b.mirrored) T = T * Mirror;
1542 return T;
1543 }
1544 if(i == c->type-1) return mscale(Id, cgi.plevel);
1545 else if(i == c->type-2) return mscale(Id, -cgi.plevel);
1546 c = where[c].first;
1547 return PIU(currentmap->adj(c, i));
1548 }
1549
hrmap_producthr::product::hrmap_product1550 hrmap_product() {
1551 current_spin_invalid = false;
1552 using hybrid::csteps;
1553 if((cspin || cmirror) && csteps) {
1554 in_underlying([&] {
1555 twisted = validate_spin();
1556 if(!twisted) { current_spin_invalid = true; return; }
1557 auto ugs = currentmap->gamestart();
1558 spins[ugs] = make_pair(
1559 cellwalker(ugs, gmod(+cspin, ugs->type), cmirror),
1560 cellwalker(ugs, gmod(-cspin, ugs->type), cmirror)
1561 );
1562 manual_celllister cl;
1563 cl.add(ugs);
1564 for(int i=0; i<isize(cl.lst); i++) {
1565 cell *c = cl.lst[i];
1566 hybrid::will_link(c);
1567 forCellEx(c2, c) cl.add(c2);
1568 }
1569 hybrid::link();
1570 });
1571 }
1572 }
1573
ray_iadjhr::product::hrmap_product1574 virtual transmatrix ray_iadj(cell *c, int i) override {
1575 if(i == c->type-2) return (mscale(Id, +cgi.plevel));
1576 if(i == c->type-1) return (mscale(Id, -cgi.plevel));
1577 transmatrix T;
1578 cell *cw = hybrid::get_where(c).first;
1579 hybrid::in_underlying_geometry([&] {
1580 T = currentmap->ray_iadj(cw, i);
1581 });
1582 return T;
1583 }
1584 };
1585
1586 EX bool current_spin_invalid, cmirror;
1587 EX int cspin;
1588
1589 /* todo might need a shiftpoint version */
inverse_exp(hyperpoint h)1590 EX hyperpoint inverse_exp(hyperpoint h) {
1591 hyperpoint res;
1592 res[2] = zlevel(h);
1593 h = zshift(h, -res[2]);
1594 ld r = hypot_d(2, h);
1595 if(hybrid::under_class() == gcEuclid) {
1596 res[0] = h[0];
1597 res[1] = h[1];
1598 }
1599 else if(r < 1e-6) {
1600 res[0] = h[0];
1601 res[1] = h[1];
1602 }
1603 else {
1604 auto c = acos_auto_clamp(h[2]);
1605 r = c / r;
1606 res[0] = h[0] * r;
1607 res[1] = h[1] * r;
1608 }
1609 return res;
1610 }
1611
direct_exp(hyperpoint h)1612 EX hyperpoint direct_exp(hyperpoint h) {
1613 hyperpoint res;
1614 ld d = hypot_d(2, h);
1615 ld cd = d == 0 ? 0 : sin_auto(d) / d;
1616 res[0] = h[0] * cd;
1617 res[1] = h[1] * cd;
1618 res[2] = cos_auto(d);
1619 return zshift(res, h[2]);
1620 }
1621
validate_spin()1622 EX bool validate_spin() {
1623 if(prod) return hybrid::in_underlying_geometry(validate_spin);
1624 if(kite::in()) return false;
1625 if(!quotient && !arcm::in()) return true;
1626 map<cell*, cellwalker> cws;
1627 manual_celllister cl;
1628 cell *start = currentmap->gamestart();
1629 cl.add(start);
1630 cws[start] = cellwalker(start, gmod(cspin, start->type), cmirror);
1631 for(int i=0; i<isize(cl.lst); i++) {
1632 cell *c = cl.lst[i];
1633 cellwalker cwc = cws.at(c);
1634 forCellIdEx(c2, j, c) {
1635 cellwalker cwc2 = cwc + j + wstep - c->c.spin(j);
1636 if(!cws.count(c2)) cws[c2] = cwc2;
1637 else if(cws[c2] != cwc2) return false;
1638 cl.add(c2);
1639 }
1640 }
1641 return true;
1642 }
1643
show_config()1644 EX void show_config() {
1645 cmode = sm::SIDE | sm::MAYDARK;
1646 gamescreen(1);
1647 dialog::init(XLAT("quotient product spaces"));
1648 dialog::addSelItem(XLAT("%1 period", "Z"), its(hybrid::csteps), 'z');
1649 dialog::add_action(hybrid::configure_period);
1650 dialog::addSelItem(XLAT("rotation"), its(cspin), 'r');
1651 dialog::add_action([] {
1652 static int s;
1653 dialog::editNumber(s, 0, 16, 1, 0, XLAT("rotation", "Z"),
1654 XLAT("Works if the underlying space is symmetric.")
1655 );
1656 dialog::reaction_final = [] {
1657 if(s == cspin) return;
1658 stop_game();
1659 cspin = s;
1660 start_game();
1661 };
1662 });
1663 dialog::addBoolItem(XLAT("reflect"), cmirror, 'f');
1664 dialog::add_action([]{
1665 stop_game();
1666 cmirror = !cmirror;
1667 start_game();
1668 });
1669 if(current_spin_invalid)
1670 dialog::addInfo("the current rotation is invalid");
1671 else
1672 dialog::addBreak(100);
1673
1674 dialog::addBreak(50);
1675 dialog::addBack();
1676 dialog::display();
1677 }
1678
1679 EX }
1680
1681 EX namespace slr {
1682
1683 EX ld range_xy = 2;
1684 EX int steps = 15;
1685
translate(hyperpoint h)1686 EX transmatrix translate(hyperpoint h) {
1687 return matrix4(
1688 h[3], -h[2], h[1], h[0],
1689 h[2], h[3], -h[0], h[1],
1690 h[1], -h[0], h[3], h[2],
1691 h[0], h[1], -h[2], h[3]
1692 );
1693 }
1694
polar(ld r,ld theta,ld phi)1695 EX hyperpoint polar(ld r, ld theta, ld phi) {
1696 return hyperpoint(sinh(r) * cos(theta-phi), sinh(r) * sin(theta-phi), cosh(r) * sin(phi), cosh(r) * cos(phi));
1697 }
1698
xyz_point(ld x,ld y,ld z)1699 EX hyperpoint xyz_point(ld x, ld y, ld z) {
1700 ld r = hypot(x, y);
1701 ld f = r ? sinh(r) / r : 1;
1702 return hyperpoint(x * f * cos(z) + y * f * sin(z), y * f * cos(z) - x * f * sin(z), cosh(r) * sin(z), cosh(r) * cos(z));
1703 }
1704
get_inverse_exp(shiftpoint h)1705 EX hyperpoint get_inverse_exp(shiftpoint h) {
1706 ld xy = hypot_d(2, h.h);
1707 ld phi = atan2(h[2], h[3]) + h.shift;
1708
1709 if(xy < 1e-6) return point31(0.,0.,phi);
1710
1711 bool flipped = phi > 0;
1712 if(flipped) phi = -phi;
1713
1714 ld SV = stretch::not_squared();
1715 ld K = -1;
1716
1717 ld alpha = flipped ? atan2(h[1], h[0]) - h.shift : atan2(h[1], -h[0]) + h.shift;
1718
1719 hyperpoint res;
1720
1721 ld fiber_barrier = atan(1/SV);
1722
1723 ld flip_barrier = atan( 1 / tanh(asinh(xy)) / SV);
1724
1725 // test the side of the flip barrier
1726
1727 int part = -1;
1728
1729 if(1) {
1730 ld kk = flip_barrier;
1731
1732 ld x_part = cos(kk);
1733 ld z_part = sin(kk);
1734
1735 ld rparam = x_part / z_part / SV;
1736
1737 ld r = atanh(rparam);
1738
1739 ld cr = cosh(r);
1740 ld sr = sinh(r);
1741
1742 // sinh(r) = xy
1743 // r = tanh(sinh(xy))
1744
1745
1746 ld z = cr * (K - 1/SV/SV);
1747
1748 ld k = M_PI/2;
1749 ld a = k / K;
1750 ld zw = xy * cr / sr;
1751 ld u = z * a;
1752
1753 ld phi1 = atan2(zw, cos(k)) - u;
1754
1755 if(phi < phi1) part = 2;
1756 }
1757
1758 if(part == -1) {
1759 ld zw = xy;
1760
1761 ld u = xy * (K - 1/SV/SV) / K;
1762 ld phi1 = atan2(zw, 1) - u;
1763
1764 if(phi > phi1) part = 0; else part = 1;
1765 }
1766
1767 if(part == 2) {
1768 ld min_k = fiber_barrier;
1769 ld max_k = flip_barrier;
1770
1771 for(int it=0; it<30; it++) {
1772 ld kk = (min_k + max_k) / 2;
1773
1774 ld x_part = cos(kk);
1775 ld z_part = sin(kk);
1776
1777 ld rparam = x_part / z_part / SV;
1778
1779 assert(rparam <= 1);
1780
1781 ld r = atanh(rparam);
1782 ld cr = cosh(r);
1783 ld sr = sinh(r);
1784
1785 ld z = cr * (K - 1/SV/SV);
1786
1787 ld k = M_PI - asin(xy / sr);
1788 ld a = k / K;
1789 ld len = a * hypot(sr, cr/SV);
1790 ld zw = xy * cr / sr;
1791 ld u = z * a;
1792
1793 ld phi1 = atan2(zw, cos(k)) - u;
1794
1795 if(phi < phi1) max_k = kk;
1796 else min_k = kk;
1797
1798 ld r_angle = alpha + u;
1799 res = point3(cos(r_angle) * x_part * len, -sin(r_angle) * x_part * len, z_part * len);
1800 }
1801 }
1802
1803 if(part == 0) {
1804 ld min_k = 0;
1805 ld max_k = fiber_barrier;
1806
1807 for(int it=0; it<30; it++) {
1808 ld kk = (min_k + max_k) / 2;
1809
1810 ld x_part = cos(kk);
1811 ld z_part = sin(kk);
1812
1813 ld rparam = x_part / z_part / SV;
1814
1815 ld cr = 1 / sqrt(rparam*rparam - 1);
1816 ld sr = rparam * cr;
1817
1818 ld z = cr * (K - 1/SV/SV);
1819
1820 ld k = asinh(xy / sr);
1821 ld a = k / K;
1822 ld len = a * hypot(sr, cr/SV);
1823
1824 ld zw = xy * cr / sr;
1825
1826 ld u = z * a;
1827 ld phi1 = atan2(zw, cosh(k)) - u;
1828
1829 if(phi > phi1) max_k = kk; else min_k = kk;
1830
1831 ld r_angle = alpha + u;
1832 res = point3(cos(r_angle) * x_part * len, -sin(r_angle) * x_part * len, z_part * len);
1833 }
1834 }
1835
1836 if(part == 1) {
1837 ld min_k = fiber_barrier;
1838 ld max_k = flip_barrier;
1839
1840 for(int it=0; it<30; it++) {
1841 ld kk = (min_k + max_k) / 2;
1842
1843 ld x_part = cos(kk);
1844 ld z_part = sin(kk);
1845
1846 ld rparam = x_part / z_part / SV;
1847
1848 ld r = atanh(rparam);
1849 ld cr = cosh(r);
1850 ld sr = sinh(r);
1851
1852 ld z = cr * (K - 1/SV/SV);
1853
1854 ld k = asin(xy / sr);
1855 ld a = k / K;
1856 ld len = a * hypot(sr, cr/SV);
1857 ld zw = xy * cr / sr;
1858 ld u = z * a;
1859
1860 ld phi1 = atan2(zw, cos(k)) - u;
1861
1862 if(isnan(phi1)) max_k = kk;
1863 else if(phi > phi1) max_k = kk;
1864 else min_k = kk;
1865
1866 ld r_angle = alpha + u;
1867 res = point3(cos(r_angle) * x_part * len, -sin(r_angle) * x_part * len, z_part * len);
1868 }
1869 }
1870
1871 if(flipped) res[0] *= -1, res[2] *= -1;
1872
1873 return res;
1874 }
1875
1876 #if ISWEB
1877 #define ITERATE " for(int it=0; it<50; it++) { if(it >= uIterations) break; "
1878 #else
1879 #define ITERATE " for(int it=0; it<uIterations; it++) {"
1880 #endif
1881
1882 EX string slshader =
1883
1884 "uniform mediump float uIndexSL;"
1885 "uniform mediump int uIterations;"
1886 "uniform mediump float uSV;"
1887
1888 "vec4 inverse_exp(vec4 h) {"
1889
1890 "float xy = length(h.xy);"
1891 "float phi = atan2(h[2], h[3]) + uIndexSL;"
1892
1893 "if(xy < 1e-6) return vec4(0.,0.,phi,1.);"
1894
1895 "vec4 res = vec4(sqrt(-1.),sqrt(-1.),sqrt(-1.),sqrt(-1.));"
1896
1897 "bool flipped = phi > 0.;"
1898
1899 "if(flipped) phi = -phi;"
1900
1901 "float alpha = flipped ? atan2(h[1], h[0]) - uIndexSL : atan2(h[1], -h[0]) + uIndexSL;"
1902
1903 "float fiber_barrier = atan(1./uSV);"
1904
1905 "float flip_barrier = atan(1. / tanh(asinh(xy)) / uSV);"
1906
1907 "int part = 0;"
1908
1909 "if(true) {"
1910 "float x_part = cos(flip_barrier);"
1911 "float z_part = sin(flip_barrier);"
1912 "float rparam = x_part / z_part / uSV;"
1913 "float r = atanh(rparam);"
1914 "float cr = cosh(r);"
1915 "float sr = sinh(r);"
1916 "float z = cr * (-1.-1./uSV/uSV);"
1917 "float k = PI/2.;"
1918 "float a = -k;"
1919 "float zw = xy * cr / sr;"
1920 "float u = z * a;"
1921 "float phi1 = atan2(zw, cos(k)) - u;"
1922 "if(phi < phi1) part = 2;"
1923 "}\n"
1924
1925 "if(part == 0) {"
1926 "float zw = xy;"
1927 "float u = xy * (1. + 1./uSV/uSV);"
1928 "float phi1 = atan2(zw, 1.) - u;"
1929 "if(phi > phi1) part = 0; else part = 1;"
1930 "}\n"
1931
1932 "if(part == 2) {"
1933 "float min_k = fiber_barrier;"
1934 "float max_k = flip_barrier;"
1935
1936 ITERATE
1937 "float kk = (min_k + max_k) / 2.;"
1938 "float x_part = cos(kk);"
1939 "float z_part = sin(kk);"
1940 "float rparam = x_part / z_part / uSV;"
1941 "float r = atanh(rparam);"
1942 "float cr = cosh(r);"
1943 "float sr = sinh(r);"
1944
1945 "float z = cr * (-1. - 1./uSV/uSV);"
1946 "float k = PI - asin(xy / sr);"
1947 "float a = -k;"
1948 "float len = a * length(vec2(sr, cr/uSV));"
1949 "float zw = xy * cr / sr;"
1950 "float u = z * a;"
1951 "float phi1 = atan2(zw, cos(k)) - u;"
1952 "if(phi < phi1) max_k = kk; else min_k = kk;"
1953 "float r_angle = alpha + u;"
1954 "res = vec4(cos(r_angle) * x_part * len, -sin(r_angle) * x_part * len, z_part * len, 1);"
1955 "}"
1956 "}\n"
1957
1958 "if(part == 0) {"
1959 "float min_k = 0.;"
1960 "float max_k = fiber_barrier;"
1961
1962 ITERATE
1963 "float kk = (min_k + max_k) / 2.;"
1964 "float x_part = cos(kk);"
1965 "float z_part = sin(kk);"
1966 "float rparam = x_part / z_part / uSV;"
1967 "float cr = 1. / sqrt(rparam*rparam - 1.);"
1968 "float sr = rparam * cr;"
1969 "float z = cr * (-1. - 1./uSV/uSV);"
1970 "float k = asinh(xy / sr);"
1971 "float a = -k;"
1972 "float len = a * length(vec2(sr, cr/uSV));"
1973 "float zw = xy * cr / sr;"
1974 "float u = z * a;"
1975 "float phi1 = atan2(zw, cosh(k)) - u;"
1976
1977 "if(phi > phi1) max_k = kk; else min_k = kk;"
1978
1979 "float r_angle = alpha + u;"
1980 "res = vec4(cos(r_angle) * x_part * len, -sin(r_angle) * x_part * len, z_part * len, 1);"
1981 "}"
1982 "}\n"
1983
1984 "if(part == 1) {"
1985 "float min_k = fiber_barrier;"
1986 "float max_k = flip_barrier;"
1987
1988 ITERATE
1989 "float kk = (min_k + max_k) / 2.;"
1990
1991 "float x_part = cos(kk);"
1992 "float z_part = sin(kk);"
1993
1994 "float rparam = x_part / z_part / uSV;"
1995
1996 "float r = atanh(rparam);"
1997 "float cr = cosh(r);"
1998 "float sr = sinh(r);"
1999
2000 "float z = cr * (-1. - 1./uSV/uSV);"
2001
2002 "float k = asin(xy / sr);"
2003 "float a = -k;"
2004 "float len = a * length(vec2(sr, cr/uSV));"
2005 "float zw = xy * cr / sr;"
2006 "float u = z * a;"
2007
2008 "float phi1 = atan2(zw, cos(k)) - u;"
2009
2010 "if(phi > phi1) max_k = kk;"
2011 "else min_k = kk;"
2012
2013 "float r_angle = alpha + u;"
2014 "res = vec4(cos(r_angle) * x_part * len, -sin(r_angle) * x_part * len, z_part * len, 1);"
2015 "}"
2016 "}\n"
2017
2018 "if(flipped) res[0] *= -1., res[2] *= -1.;"
2019
2020 "return res;"
2021 "}";
2022
2023 EX }
2024
2025 EX namespace rots {
2026 EX ld underlying_scale = 0;
2027
2028 #if MAXMDIM >= 4
uxpush(ld x)2029 EX transmatrix uxpush(ld x) {
2030 if(sl2) return xpush(x);
2031 return cspin(1, 3, x) * cspin(0, 2, x);
2032 }
2033
uypush(ld y)2034 EX transmatrix uypush(ld y) {
2035 if(sl2) return ypush(y);
2036 return cspin(0, 3, -y) * cspin(1, 2, y);
2037 }
2038
uzpush(ld z)2039 EX transmatrix uzpush(ld z) {
2040 if(sl2) return zpush(z);
2041 return cspin(3, 2, -z) * cspin(0, 1, -z);
2042 }
2043
lift_matrix(const transmatrix & T)2044 EX transmatrix lift_matrix(const transmatrix& T) {
2045 hyperpoint d;
2046 ld alpha, beta, distance;
2047 transmatrix Spin;
2048 hybrid::in_underlying_geometry([&] {
2049 hyperpoint h = tC0(T);
2050 Spin = iso_inverse(gpushxto0(h) * T);
2051 d = hr::inverse_exp(shiftless(h));
2052 alpha = atan2(Spin[0][1], Spin[0][0]);
2053 distance = hdist0(h);
2054 beta = atan2(h[1], h[0]);
2055 });
2056 for(int k=0; k<3; k++) Spin[3][k] = Spin[k][3] = 0; Spin[3][3] = 1;
2057 return spin(beta) * uxpush(distance/2) * spin(-beta+alpha);
2058 }
2059
2060 EX std::map<int, transmatrix> saved_matrices_ray;
2061
2062 struct hrmap_rotation_space : hybrid::hrmap_hybrid {
2063
2064 std::map<int, transmatrix> saved_matrices;
2065
adjhr::rots::hrmap_rotation_space2066 transmatrix adj(cell *c1, int i) override {
2067 if(i == c1->type-2) return uzpush(-cgi.plevel) * spin(-2*cgi.plevel);
2068 if(i == c1->type-1) return uzpush(+cgi.plevel) * spin(+2*cgi.plevel);
2069 cell *c2 = c1->cmove(i);
2070 #if CAP_ARCM
2071 int id1 = hybrid::underlying == gArchimedean ? arcm::id_of(c1->master) + 20 * arcm::parent_index_of(c1->master) : shvid(c1);
2072 int id2 = hybrid::underlying == gArchimedean ? arcm::id_of(c2->master) + 20 * arcm::parent_index_of(c2->master) : shvid(c2);
2073 #else
2074 int id1 = shvid(c1);
2075 int id2 = shvid(c2);
2076 #endif
2077 int j = c1->c.spin(i);
2078 int id = id1 + (id2 << 10) + (i << 20) + (j << 26);
2079 auto &M = saved_matrices[id];
2080 if(M[3][3]) return M;
2081
2082 cell *cw = where[c1].first;
2083 return M = lift_matrix(PIU(currentmap->adj(cw, i)));
2084 }
2085
relative_matrixchr::rots::hrmap_rotation_space2086 transmatrix relative_matrixc(cell *c2, cell *c1, const hyperpoint& hint) override {
2087 if(c1 == c2) return Id;
2088 if(gmatrix0.count(c2) && gmatrix0.count(c1))
2089 return inverse_shift(gmatrix0[c1], gmatrix0[c2]);
2090 for(int i=0; i<c1->type; i++) if(c1->move(i) == c2) return adj(c1, i);
2091 return Id; // not implemented yet
2092 }
2093
ray_iadjhr::rots::hrmap_rotation_space2094 transmatrix ray_iadj(cell *c1, int i) override {
2095 if(i == c1->type-1) return uzpush(-cgi.plevel) * spin(-2*cgi.plevel);
2096 if(i == c1->type-2) return uzpush(+cgi.plevel) * spin(+2*cgi.plevel);
2097 cell *c2 = c1->cmove(i);
2098 #if CAP_ARCM
2099 int id1 = hybrid::underlying == gArchimedean ? arcm::id_of(c1->master) + 20 * arcm::parent_index_of(c1->master) : shvid(c1);
2100 int id2 = hybrid::underlying == gArchimedean ? arcm::id_of(c2->master) + 20 * arcm::parent_index_of(c2->master) : shvid(c2);
2101 #else
2102 int id1 = shvid(c1);
2103 int id2 = shvid(c2);
2104 #endif
2105 int j = c1->c.spin(i);
2106 int id = id1 + (id2 << 10) + (i << 20) + (j << 26);
2107 auto &M = saved_matrices_ray[id];
2108 if(M[3][3]) return M;
2109
2110 cell *cw = hybrid::get_where(c1).first;
2111
2112 transmatrix T;
2113 hybrid::in_underlying_geometry([&] {
2114 hyperpoint h0 = get_corner_position(cw, i);
2115 hyperpoint h1 = get_corner_position(cw, (i+1));
2116 T = to_other_side(h0, h1);
2117 });
2118
2119 return M = lift_matrix(T);
2120 }
2121 };
2122
2123 /** reinterpret the given point of rotspace as a rotation matrix in the underlying geometry (note: this is the inverse) */
qtm(hyperpoint h)2124 EX transmatrix qtm(hyperpoint h) {
2125
2126 ld& x = h[0];
2127 ld& y = h[1];
2128 ld& z = h[2];
2129 ld& w = h[3];
2130
2131 ld xx = x*x;
2132 ld yy = y*y;
2133 ld zz = z*z;
2134 ld ww = w*w;
2135
2136 ld xy = x*y;
2137 ld xz = x*z;
2138 ld xw = x*w;
2139 ld yz = y*z;
2140 ld yw = y*w;
2141 ld zw = z*w;
2142
2143 transmatrix M;
2144
2145 M[0][0] = +xx - yy - zz + ww;
2146 M[1][1] = -xx + yy - zz + ww;
2147 M[2][2] = -xx - yy + zz + ww;
2148
2149 M[0][1] = -2 * (xy + zw);
2150 M[1][0] = -2 * (xy - zw);
2151
2152 M[0][2] = 2 * (xz - yw);
2153 M[2][0] = 2 * (xz + yw);
2154
2155 M[1][2] = -2 * (yz + xw);
2156 M[2][1] = -2 * (yz - xw);
2157
2158 if(hyperbolic) {
2159 swap(M[0][2], M[1][2]);
2160 swap(M[2][0], M[2][1]);
2161 M[1][2] *= -1;
2162 M[2][0] *= -1;
2163 M[2][2] = xx + yy + zz + ww;
2164 return M;
2165 }
2166
2167
2168 return M;
2169 }
2170
2171 EX bool drawing_underlying = false;
2172
draw_underlying(bool cornermode)2173 EX void draw_underlying(bool cornermode) {
2174 if(underlying_scale <= 0) return;
2175 ld d = hybrid::get_where(centerover).second;
2176 d *= cgi.plevel;
2177 transmatrix T = rots::uzpush(-d) * spin(-2*d);
2178
2179 if(det(T) < 0) T = centralsym * T;
2180
2181 if(prod) d = 0;
2182
2183 hyperpoint h = inverse(View * spin(master_to_c7_angle()) * T) * C0;
2184
2185 auto g = std::move(gmatrix);
2186 auto g0 = std::move(gmatrix0);
2187
2188 ld alpha = atan2(ortho_inverse(NLP) * point3(1, 0, 0));
2189
2190 bool inprod = prod;
2191 transmatrix pView = View;
2192 if(inprod) {
2193 pView = spin(alpha) * View;
2194 ld z = zlevel(tC0(View));
2195 for(int a=0; a<3; a++) pView[a] *= exp(-z);
2196 }
2197
2198 cell *co = hybrid::get_where(centerover).first;
2199
2200 hybrid::in_underlying_geometry([&] {
2201 cgi.require_shapes();
2202 dynamicval<int> pcc(corner_centering, cornermode ? 1 : 2);
2203 dynamicval<bool> pf(playerfound, true);
2204 dynamicval<cell*> m5(centerover, co);
2205 dynamicval<transmatrix> m2(View, inprod ? pView : ypush(0) * qtm(h));
2206 if(PURE && !inprod) View = View * pispin;
2207 View = inverse(stretch::mstretch_matrix) * spin(2*d) * View;
2208 dynamicval<shiftmatrix> m3(playerV, shiftless(Id));
2209 dynamicval<transmatrix> m4(actual_view_transform, Id);
2210 dynamicval<shiftmatrix> m6(cwtV, shiftless(Id));
2211 dynamicval<eModel> pm(pmodel, mdDisk);
2212 dynamicval<ld> pss(pconf.scale, (sphere ? 10 : euclid ? .4 : 1) * underlying_scale);
2213 dynamicval<ld> psa(pconf.alpha, sphere ? 10 : 1);
2214 dynamicval<hrmap*> p(hybrid::pmap, NULL);
2215 dynamicval<int> psr(sightrange_bonus, 0);
2216
2217 dynamicval<int> psx(vid.use_smart_range, 2);
2218 dynamicval<ld> psy(vid.smart_range_detail, 1);
2219 dynamicval<bool> pdu(drawing_underlying, true);
2220
2221 calcparam();
2222 reset_projection(); current_display->set_all(0, 0);
2223 ptds.clear();
2224 drawthemap();
2225 drawqueue();
2226 displaychr(current_display->xcenter, current_display->ycenter, 0, 24, '+', 0xFFFFFFFF);
2227 glflush();
2228 });
2229 gmatrix = std::move(g);
2230 gmatrix0 = std::move(g0);
2231 calcparam();
2232 reset_projection(); current_display->set_all(0, 0);
2233 }
2234
2235 /** @brief exponential function for both slr and Berger sphere */
2236
formula_exp(hyperpoint vel)2237 EX hyperpoint formula_exp(hyperpoint vel) {
2238 bool sp = sphere;
2239 ld K = sp ? 1 : -1;
2240
2241 ld len = hypot_d(3, vel);
2242
2243 if(vel[2] < 0) len = -len;
2244
2245 ld z_part = vel[2]/len;
2246 ld x_part = sqrt(max<ld>(1 - z_part * z_part, 0));
2247
2248 ld SV = stretch::not_squared();
2249
2250 ld rparam = x_part / z_part / SV;
2251
2252 ld beta = atan2(vel[1], vel[0]);
2253 if(len < 0) beta += M_PI;
2254
2255 if(sl2 && rparam > 1) {
2256 ld cr = 1 / sqrt(rparam*rparam - 1); // *i
2257 ld sr = rparam * cr; // *i
2258
2259 if(z_part == 0) cr = 0, sr = 1;
2260
2261 ld z = cr * (K - 1/SV/SV); // *i
2262
2263 ld a = len / hypot(sr, cr/SV); // /i
2264
2265 ld k = K*a; // /i
2266 ld u = z*a;
2267
2268 ld xy = sr * sinh(k);
2269 ld zw = cr * sinh(k);
2270
2271 return hyperpoint(K*xy * cos(u+beta), K*xy * sin(u+beta), zw * cos(u) - cosh(k) * sin(u), zw * sin(u) + cosh(k)*cos(u));
2272 }
2273
2274 else {
2275 ld r = atan_auto(rparam);
2276 ld cr = cos_auto(r);
2277 ld sr = sin_auto(r);
2278
2279 ld z = cr * (K - 1/SV/SV);
2280
2281 ld a = len / hypot(sr, cr/SV);
2282
2283 ld k = K*a;
2284 ld u = z*a;
2285
2286 ld xy = sr * sin(k);
2287 ld zw = cr * sin(k);
2288
2289 return hyperpoint(K*xy * cos(u+beta), K*xy * sin(u+beta), zw * cos(u) - cos(k) * sin(u), zw * sin(u) + cos(k)*cos(u));
2290 }
2291 }
2292
2293 #endif
2294 EX }
2295
2296 /** stretched rotation space (S3 or SLR) */
2297 EX namespace stretch {
2298
2299 EX ld factor;
2300
2301 EX bool mstretch;
2302
2303 EX transmatrix m_itoa, m_atoi, m_pd;
2304 EX ld ms_christoffel[3][3][3];
2305
2306 EX transmatrix mstretch_matrix;
2307
enable_mstretch()2308 EX void enable_mstretch() {
2309 mstretch = true;
2310
2311 for(int a=0; a<4; a++)
2312 for(int b=0; b<4; b++)
2313 if(a==3 || b==3) m_atoi[a][b] = (a==b);
2314
2315 m_itoa = inverse3(m_atoi);
2316
2317 for(int a=0; a<4; a++)
2318 for(int b=0; b<4; b++)
2319 if(a==3 || b==3)
2320 m_itoa[a][b] = m_atoi[a][b] = 0;
2321
2322 for(int j=0; j<3; j++)
2323 for(int k=0; k<3; k++) {
2324 m_pd[j][k] = 0;
2325 for(int i=0; i<3; i++)
2326 m_pd[j][k] += m_atoi[i][j] * m_atoi[i][k];
2327 }
2328
2329 auto& c = ms_christoffel;
2330
2331 ld A00 = m_pd[0][0];
2332 ld A11 = m_pd[1][1];
2333 ld A22 = m_pd[2][2];
2334 ld A01 = m_pd[0][1] + m_pd[1][0];
2335 ld A02 = m_pd[0][2] + m_pd[2][0];
2336 ld A12 = m_pd[2][1] + m_pd[1][2];
2337 ld B01 = A01 * A01;
2338 ld B02 = A02 * A02;
2339 ld B12 = A12 * A12;
2340 ld B00 = A00 * A00;
2341 ld B11 = A11 * A11;
2342 ld B22 = A22 * A22;
2343
2344 ld den = (-4*A00*A11*A22 + A00*B12 + B01*A22 - A01*A02*A12 + B02*A11);
2345
2346 if(sl2) {
2347 c[ 0 ][ 0 ][ 0 ] = (A01*(A01*A12 - 2*A02*A11) - A02*(2*A01*A22 - A02*A12))/den;
2348 c[ 0 ][ 0 ][ 1 ] = (A00*A01*A12 - 2*A00*A02*A11 - A01*A11*A12 + A01*A12*A22 + 2*A02*B11 + 2*A02*A11*A22 - A02*B12)/-den ;
2349 c[ 0 ][ 0 ][ 2 ] = (-A01*(4*A11*A22 - B12)/2 + A12*(A01*A12 - 2*A02*A11)/2 - (A00 + A22)*(2*A01*A22 - A02*A12))/den;
2350 c[ 0 ][ 1 ][ 0 ] = (A00*A01*A12 - 2*A00*A02*A11 - A01*A11*A12 + A01*A12*A22 + 2*A02*B11 + 2*A02*A11*A22 - A02*B12)/-den ;
2351 c[ 0 ][ 1 ][ 1 ] = -(A01*(A01*A12 - 2*A02*A11) + A12*(4*A11*A22 - B12))/den;
2352 c[ 0 ][ 1 ][ 2 ] = (B01*A22 - B02*A11 + 4*B11*A22 - A11*B12 + 4*A11*B22 - B12*A22)/-den ;
2353 c[ 0 ][ 2 ][ 0 ] = (-A01*(4*A11*A22 - B12)/2 + A12*(A01*A12 - 2*A02*A11)/2 - (A00 + A22)*(2*A01*A22 - A02*A12))/den;
2354 c[ 0 ][ 2 ][ 1 ] = (B01*A22 - B02*A11 + 4*B11*A22 - A11*B12 + 4*A11*B22 - B12*A22)/-den ;
2355 c[ 0 ][ 2 ][ 2 ] = -(A02*(2*A01*A22 - A02*A12) + A12*(4*A11*A22 - B12))/den;
2356 c[ 1 ][ 0 ][ 0 ] = (-A01*(2*A00*A12 - A01*A02) + A02*(4*A00*A22 - B02))/den;
2357 c[ 1 ][ 0 ][ 1 ] = (A02*(2*A01*A22 - A02*A12)/2 + A12*(4*A00*A22 - B02)/2 + (A00 - A11)*(2*A00*A12 - A01*A02))/den;
2358 c[ 1 ][ 0 ][ 2 ] = (-4*B00*A22 + A00*B02 + A00*B12 - 4*A00*B22 - B01*A22 + B02*A22)/-den ;
2359 c[ 1 ][ 1 ][ 0 ] = (A02*(2*A01*A22 - A02*A12)/2 + A12*(4*A00*A22 - B02)/2 + (A00 - A11)*(2*A00*A12 - A01*A02))/den;
2360 c[ 1 ][ 1 ][ 1 ] = (A01*(2*A00*A12 - A01*A02) + A12*(2*A01*A22 - A02*A12))/den;
2361 c[ 1 ][ 1 ][ 2 ] = (A01*(4*A00*A22 - B02)/2 + A02*(2*A00*A12 - A01*A02)/2 + (A11 + A22)*(2*A01*A22 - A02*A12))/den;
2362 c[ 1 ][ 2 ][ 0 ] = (-4*B00*A22 + A00*B02 + A00*B12 - 4*A00*B22 - B01*A22 + B02*A22)/-den ;
2363 c[ 1 ][ 2 ][ 1 ] = (A01*(4*A00*A22 - B02)/2 + A02*(2*A00*A12 - A01*A02)/2 + (A11 + A22)*(2*A01*A22 - A02*A12))/den;
2364 c[ 1 ][ 2 ][ 2 ] = (A02*(4*A00*A22 - B02) + A12*(2*A01*A22 - A02*A12))/den;
2365 c[ 2 ][ 0 ][ 0 ] = (A01*(4*A00*A11 - B01) - A02*(2*A00*A12 - A01*A02))/den;
2366 c[ 2 ][ 0 ][ 1 ] = (4*B00*A11 - A00*B01 - 4*A00*B11 + A00*B12 + B01*A11 - B02*A11)/-den ;
2367 c[ 2 ][ 0 ][ 2 ] = (-A01*(A01*A12 - 2*A02*A11)/2 + A12*(4*A00*A11 - B01)/2 - (A00 + A22)*(2*A00*A12 - A01*A02))/den;
2368 c[ 2 ][ 1 ][ 0 ] = (4*B00*A11 - A00*B01 - 4*A00*B11 + A00*B12 + B01*A11 - B02*A11)/-den ;
2369 c[ 2 ][ 1 ][ 1 ] = -(A01*(4*A00*A11 - B01) + A12*(A01*A12 - 2*A02*A11))/den;
2370 c[ 2 ][ 1 ][ 2 ] = (A00*A01*A12 + 2*A00*A02*A11 - B01*A02 + A01*A11*A12 + A01*A12*A22 - 2*A02*B11 - 2*A02*A11*A22)/-den ;
2371 c[ 2 ][ 2 ][ 0 ] = (-A01*(A01*A12 - 2*A02*A11)/2 + A12*(4*A00*A11 - B01)/2 - (A00 + A22)*(2*A00*A12 - A01*A02))/den;
2372 c[ 2 ][ 2 ][ 1 ] = (A00*A01*A12 + 2*A00*A02*A11 - B01*A02 + A01*A11*A12 + A01*A12*A22 - 2*A02*B11 - 2*A02*A11*A22)/-den ;
2373 c[ 2 ][ 2 ][ 2 ] = -(A02*(2*A00*A12 - A01*A02) + A12*(A01*A12 - 2*A02*A11))/den;
2374 }
2375 else {
2376 c[ 0 ][ 0 ][ 0 ] = (A01*(A01*A12 - 2*A02*A11) + A02*(2*A01*A22 - A02*A12))/den ;
2377 c[ 0 ][ 0 ][ 1 ] = (A02*(4*A11*A22 - B12)/2 + A12*(2*A01*A22 - A02*A12)/2 - (A00 - A11)*(A01*A12 - 2*A02*A11))/den ;
2378 c[ 0 ][ 0 ][ 2 ] = (-A01*(4*A11*A22 - B12)/2 + A12*(A01*A12 - 2*A02*A11)/2 - (A00 - A22)*(2*A01*A22 - A02*A12))/den ;
2379 c[ 0 ][ 1 ][ 0 ] = (A02*(4*A11*A22 - B12)/2 + A12*(2*A01*A22 - A02*A12)/2 - (A00 - A11)*(A01*A12 - 2*A02*A11))/den ;
2380 c[ 0 ][ 1 ][ 1 ] = (-A01*(A01*A12 - 2*A02*A11) + A12*(4*A11*A22 - B12))/den ;
2381 c[ 0 ][ 1 ][ 2 ] = (B01*A22 - B02*A11 + 4*B11*A22 - A11*B12 - 4*A11*B22 + B12*A22)/(4*A00*A11*A22 - A00*B12 - B01*A22 + A01*A02*A12 - B02*A11) ;
2382 c[ 0 ][ 2 ][ 0 ] = (-A01*(4*A11*A22 - B12)/2 + A12*(A01*A12 - 2*A02*A11)/2 - (A00 - A22)*(2*A01*A22 - A02*A12))/den ;
2383 c[ 0 ][ 2 ][ 1 ] = (B01*A22 - B02*A11 + 4*B11*A22 - A11*B12 - 4*A11*B22 + B12*A22)/(4*A00*A11*A22 - A00*B12 - B01*A22 + A01*A02*A12 - B02*A11) ;
2384 c[ 0 ][ 2 ][ 2 ] = -(A02*(2*A01*A22 - A02*A12) + A12*(4*A11*A22 - B12))/den ;
2385 c[ 1 ][ 0 ][ 0 ] = -(A01*(2*A00*A12 - A01*A02) + A02*(4*A00*A22 - B02))/den ;
2386 c[ 1 ][ 0 ][ 1 ] = (-A02*(2*A01*A22 - A02*A12)/2 - A12*(4*A00*A22 - B02)/2 + (A00 - A11)*(2*A00*A12 - A01*A02))/den ;
2387 c[ 1 ][ 0 ][ 2 ] = (-4*B00*A22 + A00*B02 + A00*B12 + 4*A00*B22 - B01*A22 - B02*A22)/(4*A00*A11*A22 - A00*B12 - B01*A22 + A01*A02*A12 - B02*A11) ;
2388 c[ 1 ][ 1 ][ 0 ] = (-A02*(2*A01*A22 - A02*A12)/2 - A12*(4*A00*A22 - B02)/2 + (A00 - A11)*(2*A00*A12 - A01*A02))/den ;
2389 c[ 1 ][ 1 ][ 1 ] = (A01*(2*A00*A12 - A01*A02) - A12*(2*A01*A22 - A02*A12))/den ;
2390 c[ 1 ][ 1 ][ 2 ] = (A01*(4*A00*A22 - B02)/2 + A02*(2*A00*A12 - A01*A02)/2 + (A11 - A22)*(2*A01*A22 - A02*A12))/den ;
2391 c[ 1 ][ 2 ][ 0 ] = (-4*B00*A22 + A00*B02 + A00*B12 + 4*A00*B22 - B01*A22 - B02*A22)/(4*A00*A11*A22 - A00*B12 - B01*A22 + A01*A02*A12 - B02*A11) ;
2392 c[ 1 ][ 2 ][ 1 ] = (A01*(4*A00*A22 - B02)/2 + A02*(2*A00*A12 - A01*A02)/2 + (A11 - A22)*(2*A01*A22 - A02*A12))/den ;
2393 c[ 1 ][ 2 ][ 2 ] = (A02*(4*A00*A22 - B02) + A12*(2*A01*A22 - A02*A12))/den ;
2394 c[ 2 ][ 0 ][ 0 ] = (A01*(4*A00*A11 - B01) + A02*(2*A00*A12 - A01*A02))/den ;
2395 c[ 2 ][ 0 ][ 1 ] = (4*B00*A11 - A00*B01 - 4*A00*B11 - A00*B12 + B01*A11 + B02*A11)/(4*A00*A11*A22 - A00*B12 - B01*A22 + A01*A02*A12 - B02*A11) ;
2396 c[ 2 ][ 0 ][ 2 ] = (-A01*(A01*A12 - 2*A02*A11)/2 + A12*(4*A00*A11 - B01)/2 - (A00 - A22)*(2*A00*A12 - A01*A02))/den ;
2397 c[ 2 ][ 1 ][ 0 ] = (4*B00*A11 - A00*B01 - 4*A00*B11 - A00*B12 + B01*A11 + B02*A11)/(4*A00*A11*A22 - A00*B12 - B01*A22 + A01*A02*A12 - B02*A11) ;
2398 c[ 2 ][ 1 ][ 1 ] = (-A01*(4*A00*A11 - B01) + A12*(A01*A12 - 2*A02*A11))/den ;
2399 c[ 2 ][ 1 ][ 2 ] = (A00*A01*A12 + 2*A00*A02*A11 - B01*A02 + A01*A11*A12 - A01*A12*A22 - 2*A02*B11 + 2*A02*A11*A22)/(4*A00*A11*A22 - A00*B12 - B01*A22 + A01*A02*A12 - B02*A11) ;
2400 c[ 2 ][ 2 ][ 0 ] = (-A01*(A01*A12 - 2*A02*A11)/2 + A12*(4*A00*A11 - B01)/2 - (A00 - A22)*(2*A00*A12 - A01*A02))/den ;
2401 c[ 2 ][ 2 ][ 1 ] = (A00*A01*A12 + 2*A00*A02*A11 - B01*A02 + A01*A11*A12 - A01*A12*A22 - 2*A02*B11 + 2*A02*A11*A22)/(4*A00*A11*A22 - A00*B12 - B01*A22 + A01*A02*A12 - B02*A11) ;
2402 c[ 2 ][ 2 ][ 2 ] = -(A02*(2*A00*A12 - A01*A02) + A12*(A01*A12 - 2*A02*A11))/den ;
2403 }
2404
2405 for(int i=0; i<3; i++)
2406 for(int j=0; j<3; j++)
2407 for(int k=0; k<3; k++)
2408 if(c[i][j][k])
2409 println(hlog, tie(i,j,k), " : ", c[i][j][k]);
2410
2411
2412 println(hlog, "ATOI = ", m_atoi);
2413 println(hlog, "ITOA = ", m_itoa, " vs ", 1/not_squared());
2414 println(hlog, "PD = ", m_pd, " vs ", factor);
2415
2416 ray::reset_raycaster();
2417 }
2418
applicable()2419 EX bool applicable() {
2420 return rotspace || (cgflags & qSTRETCHABLE);
2421 }
2422
in()2423 EX bool in() {
2424 return (factor || mstretch) && applicable();
2425 }
2426
translate(hyperpoint h)2427 EX transmatrix translate(hyperpoint h) {
2428 if(!sphere) return slr::translate(h);
2429 return matrix4(
2430 h[3], -h[2], h[1], h[0],
2431 h[2], h[3], -h[0], h[1],
2432 -h[1], h[0], h[3], h[2],
2433 -h[0], -h[1], -h[2], h[3]
2434 );
2435 }
2436
itranslate(hyperpoint h)2437 EX transmatrix itranslate(hyperpoint h) {
2438 h[0] = -h[0];
2439 h[1] = -h[1];
2440 h[2] = -h[2];
2441 if(!sphere) return slr::translate(h);
2442 return translate(h);
2443 }
2444
mulz(const hyperpoint at,const hyperpoint velocity,ld zf)2445 hyperpoint mulz(const hyperpoint at, const hyperpoint velocity, ld zf) {
2446 auto vel = itranslate(at) * velocity;
2447 vel[2] *= zf;
2448 return translate(at) * vel;
2449 }
2450
squared()2451 EX ld squared() {
2452 return abs(1 + factor);
2453 }
2454
not_squared()2455 EX ld not_squared() {
2456 return sqrt(squared());
2457 }
2458
isometric_to_actual(const hyperpoint at,const hyperpoint velocity)2459 EX hyperpoint isometric_to_actual(const hyperpoint at, const hyperpoint velocity) {
2460 if(mstretch)
2461 return translate(at) * m_itoa * itranslate(at) * velocity;
2462 else
2463 return mulz(at, velocity, 1/not_squared());
2464 }
2465
actual_to_isometric(const hyperpoint at,const hyperpoint velocity)2466 EX hyperpoint actual_to_isometric(const hyperpoint at, const hyperpoint velocity) {
2467 if(mstretch)
2468 return translate(at) * m_atoi * itranslate(at) * velocity;
2469 else
2470 return mulz(at, velocity, not_squared());
2471 }
2472
christoffel(const hyperpoint at,const hyperpoint velocity,const hyperpoint transported)2473 EX hyperpoint christoffel(const hyperpoint at, const hyperpoint velocity, const hyperpoint transported) {
2474
2475 auto vel = itranslate(at) * velocity;
2476 auto tra = itranslate(at) * transported;
2477
2478 hyperpoint c;
2479
2480 if(mstretch) {
2481 c = Hypc;
2482 for(int i=0; i<3; i++)
2483 for(int j=0; j<3; j++)
2484 for(int k=0; k<3; k++)
2485 c[i] += vel[j] * tra[k] * ms_christoffel[i][j][k];
2486 }
2487
2488 else {
2489 auto K = factor;
2490 c[0] = (sphere ? -K : K+2) * (vel[1] * tra[2] + vel[2] * tra[1]);
2491 c[1] = (sphere ? K : -(K+2)) * (vel[0] * tra[2] + vel[2] * tra[0]);
2492 c[2] = 0;
2493 c[3] = 0;
2494 }
2495
2496 return translate(at) * c;
2497 }
2498
sqnorm(hyperpoint at,hyperpoint h)2499 EX ld sqnorm(hyperpoint at, hyperpoint h) {
2500 if(sphere)
2501 return sqhypot_d(4, h);
2502 h = itranslate(at) * h;
2503 return h[0] * h[0] + h[1] * h[1] + h[2] * h[2];
2504 }
2505
inverse_exp_all(hyperpoint h,int generations)2506 EX vector<hyperpoint> inverse_exp_all(hyperpoint h, int generations) {
2507
2508 vector<hyperpoint> res;
2509
2510 ld SV = stretch::not_squared();
2511
2512 if(stretch::factor == 0) {
2513 ld d = hypot_d(3, h);
2514 if(h[3] >= 1 || h[3] <= -1|| d == 0) return res;
2515 ld a = acos(h[3]);
2516
2517 res.push_back(point31(h[0] * a / d, h[1] * a / d, h[2] * a / d));
2518
2519 a = a - 2 * M_PI;
2520
2521 res.push_back(point31(h[0] * a / d, h[1] * a / d, h[2] * a / d));
2522
2523 return res;
2524 }
2525
2526 if(h[0] == 0 && h[1] == 0) {
2527 ld a = atan2(h[2], h[3]);
2528
2529 for(int it=-generations; it<generations; it++) {
2530 res.push_back(point31(0, 0, (a + 2 * M_PI * it) * SV));
2531 }
2532
2533 return res;
2534 }
2535
2536 ld xy = hypot_d(2, h);
2537
2538 ld base_min_a = asin(xy);
2539 ld base_max_a = M_PI - base_min_a;
2540
2541 ld seek = M_PI/2-atan2(h[3], h[2]);
2542
2543 auto ang = [&] (ld a) {
2544 ld rp = xy / sin(a);
2545 ld co = abs(rp) >= 1 ? 0 : sqrt(1-rp*rp);
2546
2547 return atan2(co * sin(a), cos(a)) - co * (1 - 1/SV/SV) * a;
2548
2549 // while(a0 > M_PI) a0 -= 2 * M_PI;
2550 // while(a0 < -M_PI) a0 += 2 * M_PI;
2551 };
2552
2553 for(int shift=-generations; shift<generations; shift++) {
2554 ld min_a = base_min_a + M_PI * shift;
2555 ld max_a = base_max_a + M_PI * shift;
2556
2557 ld ang_min = ang(min_a);
2558 ld ang_max = ang(max_a);
2559
2560 for(int mi=0; mi<2; mi++) {
2561 // 0 : minimum, 1 : maximum
2562 ld tl = min_a, tr = max_a;
2563 for(int it=0; it<20; it++) {
2564 ld t1 = tl * .51 + tr * .49;
2565 ld t2 = tl * .49 + tr * .51;
2566 if((ang(t1) < ang(t2)) == mi)
2567 tr = t1;
2568 else
2569 tl = t2;
2570 }
2571 ld extreme = (tl + tr) / 2;
2572 ld ang_extreme = ang(extreme);
2573 for(int t=0; t<2; t++) {
2574 ld mmin = t == 0 ? min_a : extreme;
2575 ld mmax = t == 0 ? extreme : max_a;
2576 ld vmin = t == 0 ? ang_min : ang_extreme;
2577 ld vmax = t == 0 ? ang_extreme : ang_max;
2578
2579 // make it increasing
2580 if(t != mi) swap(mmin, mmax), swap(vmin, vmax);
2581
2582 // println(hlog, "*** ", mi, t, " ** ", tie(min_a, ang_min), tie(extreme, ang_extreme), tie(max_a, ang_max), " -> ", vmin, " to ", vmax);
2583
2584 int cmin = ceil((vmin - seek) / 2 / M_PI);
2585 int cmax = floor((vmax - seek) / 2 / M_PI);
2586 for(int c = cmin; c <= cmax; c++) {
2587 ld cseek = seek + c * 2 * M_PI;
2588
2589 for(int it=0; it<40; it++) {
2590
2591 ld a = (mmin + mmax) / 2;
2592
2593 ld cros = ang(a);
2594 if(cros > cseek) mmax = a; else mmin = a;
2595 }
2596
2597 ld a = (mmin + mmax) / 2;
2598
2599 ld r = asin_clamp( xy / sin(a) );
2600
2601 ld z_part = 1;
2602 ld x_part = SV * tan(r);
2603
2604 ld db = hypot(x_part, z_part);
2605 x_part /= db;
2606 z_part /= db;
2607
2608 ld alpha = atan2(-h[1], h[0]);
2609
2610 ld z = cos(r) * (1 - 1/SV/SV);
2611 ld u = z * a;
2612
2613 ld r_angle = alpha + u;
2614
2615 ld len = a * hypot(sin_auto(r), cos_auto(r)/SV);
2616
2617 auto answer = point3(cos(r_angle) * x_part * len, -sin(r_angle) * x_part * len, z_part * len);
2618
2619 // int id = (shift << 10) + (mi << 9) + (t << 8) + c;
2620
2621 /*
2622 auto f = formula_exp(answer);
2623
2624 ld err = sqhypot_d(4, f - h);
2625
2626 println(hlog, "************************* ", answer, ": error = ", err, " id = ", id, " params = ", tie(shift, mi, t, c));
2627 */
2628
2629 res.emplace_back(answer);
2630 }
2631 }
2632 }
2633 }
2634
2635 return res;
2636 }
2637
2638
2639 EX }
2640
2641 EX namespace nisot {
2642
christoffel(const hyperpoint at,const hyperpoint velocity,const hyperpoint transported)2643 EX hyperpoint christoffel(const hyperpoint at, const hyperpoint velocity, const hyperpoint transported) {
2644 if(nil) return nilv::christoffel(at, velocity, transported);
2645 #if CAP_SOLV
2646 else if(sn::in()) return sn::christoffel(at, velocity, transported);
2647 #endif
2648 else if(stretch::in() || sl2) return stretch::christoffel(at, velocity, transported);
2649 else return point3(0, 0, 0);
2650 }
2651
in_table_range(hyperpoint h)2652 EX bool in_table_range(hyperpoint h) {
2653 #if CAP_SOLV
2654 if(sol) return sn::in_table_range(h);
2655 #endif
2656 return true;
2657 }
2658
get_acceleration(const hyperpoint & at,const hyperpoint & vel)2659 EX hyperpoint get_acceleration(const hyperpoint& at, const hyperpoint& vel) {
2660 return christoffel(at, vel, vel);
2661 }
2662
geodesic_step(hyperpoint & at,hyperpoint & vel)2663 EX void geodesic_step(hyperpoint& at, hyperpoint& vel) {
2664 /* RK4 method */
2665 auto acc1 = get_acceleration(at, vel);
2666 auto acc2 = get_acceleration(at + vel/2, vel + acc1/2);
2667 auto acc3 = get_acceleration(at + vel/2 + acc1/4, vel + acc2/2);
2668 auto acc4 = get_acceleration(at + vel + acc2/2, vel + acc3);
2669
2670 at += vel + (acc1+acc2+acc3)/6;
2671 vel += (acc1+2*acc2+2*acc3+acc4)/6;
2672 }
2673
2674 EX int rk_steps = 20;
2675
numerical_exp(hyperpoint v)2676 EX hyperpoint numerical_exp(hyperpoint v) {
2677 hyperpoint at = point31(0, 0, 0);
2678 v /= rk_steps;
2679 v[3] = 0;
2680 for(int i=0; i<rk_steps; i++) geodesic_step(at, v);
2681 return at;
2682 }
2683
parallel_transport_bare(transmatrix Pos,hyperpoint h)2684 EX transmatrix parallel_transport_bare(transmatrix Pos, hyperpoint h) {
2685
2686 bool stretch = stretch::in() || sl2;
2687
2688 h[3] = 0;
2689
2690 if(stretch::in() && stretch::mstretch)
2691 Pos = stretch::mstretch_matrix * Pos;
2692
2693 auto tPos = transpose(Pos);
2694
2695 h = Pos * h;
2696
2697 int steps = rk_steps;
2698 h /= steps;
2699
2700 auto& at = tPos[3];
2701 auto& vel = h;
2702
2703 array<ld, 4> ms;
2704
2705 if(stretch) {
2706 for(int i=0; i<3; i++) {
2707 ms[i] = stretch::sqnorm(at, tPos[i]);
2708 tPos[i] = stretch::isometric_to_actual(at, tPos[i]);
2709 }
2710 ms[3] = stretch::sqnorm(at, vel);
2711 if(!ms[3]) return Pos;
2712 vel = stretch::isometric_to_actual(at, vel);
2713 }
2714
2715 for(int i=0; i<steps; i++) {
2716 auto acc1 = get_acceleration(at, vel);
2717 auto at1 = at + vel/2; auto vel1 = vel + acc1/2;
2718 auto acc2 = get_acceleration(at1, vel1);
2719 auto at2 = at1 + acc1/4; auto vel2 = vel + acc2/2;
2720 auto acc3 = get_acceleration(at2, vel2);
2721 auto at3 = at + vel + acc2/2; auto vel3 = vel + acc3;
2722 auto acc4 = get_acceleration(at3, vel3);
2723
2724 for(int j=0; j<3; j++) {
2725 auto& tra = tPos[j];
2726
2727 auto tacc1 = christoffel(at, vel, tra);
2728 auto tacc2 = christoffel(at1, vel1, tra + tacc1/2);
2729 auto tacc3 = christoffel(at2, vel2, tra + tacc2/2);
2730 auto tacc4 = christoffel(at3, vel3, tra + tacc3);
2731
2732 tra += (tacc1+tacc2*2+tacc3*2+tacc4) / 6;
2733 }
2734
2735 at += vel + (acc1+acc2+acc3)/6;
2736 vel += (acc1+2*acc2+2*acc3+acc4)/6;
2737
2738 if(stretch) {
2739 at = normalize(at);
2740
2741 auto fix = [&] (hyperpoint& h, ld& m) {
2742 h = stretch::itranslate(at) * h;
2743 h[3] = 0;
2744 ld m1;
2745 if(stretch::mstretch) {
2746 m1 = 0;
2747 for(int i=0; i<3; i++) for(int j=0; j<3; j++)
2748 m1 += h[i] * stretch::m_pd[i][j] * h[j];
2749 }
2750 else
2751 m1 = h[0] * h[0] + h[1] * h[1] + h[2] * h[2] * stretch::squared();
2752 h /= sqrt(m1/m);
2753 h = stretch::translate(at) * h;
2754 };
2755
2756 for(int i=0; i<3; i++) fix(tPos[i], ms[i]);
2757 fix(vel, ms[3]);
2758 }
2759
2760 }
2761
2762 if(stretch) {
2763 vel = stretch::actual_to_isometric(at, vel);
2764 for(int i=0; i<3; i++) tPos[i] = stretch::actual_to_isometric(at, tPos[i]);
2765 }
2766
2767 Pos = transpose(tPos);
2768
2769 if(stretch::in() && stretch::mstretch)
2770 Pos = inverse(stretch::mstretch_matrix) * Pos;
2771
2772 return Pos;
2773 }
2774
fixmatrix(transmatrix & T)2775 EX void fixmatrix(transmatrix& T) {
2776 if(sphere) return hr::fixmatrix(T);
2777 transmatrix push = eupush( tC0(T) );
2778 transmatrix push_back = eupush(tC0(T), -1);
2779 transmatrix gtl = push_back * T;
2780 fix_rotation(gtl);
2781 T = push * gtl;
2782 }
2783
parallel_transport(const transmatrix Position,const hyperpoint direction)2784 EX transmatrix parallel_transport(const transmatrix Position, const hyperpoint direction) {
2785 auto P = Position;
2786 nisot::fixmatrix(P);
2787 if(!geodesic_movement) return eupush(Position * translate(-direction) * inverse(Position) * C0, -1) * Position;
2788 return parallel_transport_bare(P, direction);
2789 }
2790
spin_towards(const transmatrix Position,const hyperpoint goal,flagtype prec IS (pNORMAL))2791 EX transmatrix spin_towards(const transmatrix Position, const hyperpoint goal, flagtype prec IS(pNORMAL)) {
2792
2793 hyperpoint at = tC0(Position);
2794 transmatrix push_back = translate(at, -1);
2795 hyperpoint back_goal = push_back * goal;
2796 back_goal = inverse_exp(shiftless(back_goal), prec);
2797
2798 transmatrix back_Position = push_back * Position;
2799
2800 return rspintox(inverse(back_Position) * back_goal);
2801 }
2802
new_map()2803 EX hrmap *new_map() {
2804 #if CAP_SOLV
2805 if(sn::in()) return new sn::hrmap_solnih;
2806 #endif
2807 if(prod) return new product::hrmap_product;
2808 #if MAXMDIM >= 4
2809 if(nil) return new nilv::hrmap_nil;
2810 if(hybri) return new rots::hrmap_rotation_space;
2811 #endif
2812 return NULL;
2813 }
2814
2815 #if CAP_COMMANDLINE
__anon0a4104a42002() 2816 auto config = addHook(hooks_args, 0, [] () {
2817 using namespace arg;
2818 #if CAP_SOLV
2819 if(argis("-solrange")) {
2820 shift_arg_formula(sn::solrange_xy);
2821 shift_arg_formula(sn::solrange_z);
2822 return 0;
2823 }
2824 #endif
2825 if(argis("-slrange")) {
2826 shift_arg_formula(slr::range_xy);
2827 return 0;
2828 }
2829 #if CAP_SOLV
2830 else if(argis("-fsol")) {
2831 shift(); sn::solt.fname = args();
2832 return 0;
2833 }
2834 else if(argis("-nihsol")) {
2835 shift(); sn::niht.fname = args();
2836 return 0;
2837 }
2838 #endif
2839 else if(argis("-solgeo")) {
2840 geodesic_movement = true;
2841 pmodel = mdGeodesic;
2842 return 0;
2843 }
2844 else if(argis("-solnogeo")) {
2845 geodesic_movement = false;
2846 pmodel = mdPerspective;
2847 return 0;
2848 }
2849 else if(argis("-product")) {
2850 PHASEFROM(2);
2851 set_geometry(gProduct);
2852 return 0;
2853 }
2854 else if(argis("-s2xe")) {
2855 PHASEFROM(2);
2856 shift(); s2xe::qrings = argi();
2857 return 0;
2858 }
2859 else if(argis("-rotspace")) {
2860 PHASEFROM(2);
2861 set_geometry(gRotSpace);
2862 return 0;
2863 }
2864 else if(argis("-rot_uscale")) {
2865 PHASEFROM(2);
2866 shift_arg_formula(rots::underlying_scale);
2867 return 0;
2868 }
2869 else if(argis("-nilperiod")) {
2870 PHASEFROM(2);
2871 if(nil) stop_game();
2872 for(int a=0; a<3; a++) { shift(); nilv::nilperiod[a] = argi(); }
2873 nilv::set_flags();
2874 return 0;
2875 }
2876 else if(argis("-nilwidth")) {
2877 PHASEFROM(2);
2878 shift_arg_formula(nilv::nilwidth);
2879 return 0;
2880 }
2881 else if(argis("-nilh")) {
2882 PHASEFROM(2);
2883 stop_game();
2884 shift(); ginf[gNil].sides = argi();
2885 nilv::set_flags();
2886 start_game();
2887 }
2888 else if(argis("-rk-steps")) {
2889 PHASEFROM(2);
2890 shift(); rk_steps = argi();
2891 return 0;
2892 }
2893 else if(argis("-nilv")) {
2894 PHASEFROM(2);
2895 if(nil) stop_game();
2896 shift();
2897 ginf[gNil].sides = argi();
2898 return 0;
2899 }
2900 #if CAP_SOLV
2901 else if(argis("-catperiod")) {
2902 PHASEFROM(2);
2903 if(sol) stop_game();
2904 shift(); asonov::period_xy = argi();
2905 shift(); asonov::period_z = argi();
2906 asonov::set_flags();
2907 return 0;
2908 }
2909 #endif
2910 else if(argis("-prodperiod")) {
2911 PHASEFROM(2);
2912 if(prod) stop_game();
2913 shift(); hybrid::csteps = argi();
2914 hybrid::reconfigure();
2915 return 0;
2916 }
2917 else if(argis("-rot-stretch")) {
2918 PHASEFROM(2);
2919 shift_arg_formula(stretch::factor, ray::reset_raycaster);
2920 return 0;
2921 }
2922 else if(argis("-mstretch")) {
2923 PHASEFROM(2);
2924 auto& M = stretch::m_atoi;
2925 M = Id;
2926 stretch::enable_mstretch();
2927 while(true) {
2928 shift();
2929 string s = args();
2930 if(isize(s) == 2 && among(s[0], 'a', 'b','c') && among(s[1], 'a', 'b', 'c'))
2931 shift_arg_formula(M[s[0]-'a'][s[1]-'a'], stretch::enable_mstretch);
2932 else break;
2933 }
2934 // shift_arg_formula(stretch::yfactor, ray::reset_raycaster);
2935 return 0;
2936 }
2937 else if(argis("-mstretch1")) {
2938 PHASEFROM(2);
2939 auto& M = stretch::m_atoi;
2940 M = Id;
2941 M[2][2] = stretch::not_squared();
2942 stretch::enable_mstretch();
2943 // shift_arg_formula(stretch::yfactor, ray::reset_raycaster);
2944 return 0;
2945 }
2946 else if(argis("-prodturn")) {
2947 PHASEFROM(2);
2948 if(prod) stop_game();
2949 shift(); product::cspin = argi();
2950 shift(); product::cmirror = argi();
2951 return 0;
2952 }
2953 return 1;
2954 });
2955 #endif
2956 }
2957
2958 }
2959