1 // Hyperbolic Rogue -- basic computations in non-Euclidean geometry
2 // Copyright (C) 2011-2019 Zeno Rogue, see 'hyper.cpp' for details
3
4 /** \file hyperpoint.cpp
5 * \brief basic computations in non-Euclidean geometry
6 *
7 * This implements hyperpoint (a point in non-Euclidean space), transmatrix (a transformation matrix),
8 * and various basic routines related to them: rotations, translations, inverses and determinants, etc.
9 * For nonisotropic geometries, it rather refers to nonisotropic.cpp.
10 */
11
12 #include "hyper.h"
13 namespace hr {
14
15 #if HDR
16 static const ld full_circle = 2 * M_PI;
17 static const ld quarter_circle = M_PI / 2;
18 static const ld degree = M_PI / 180;
19 static const ld golden_phi = (sqrt(5)+1)/2;
20 static const ld log_golden_phi = log(golden_phi);
21 #endif
22
23 eGeometry geometry;
24 eVariation variation;
25
26
27 #if HDR
28 /** \brief A point in our continuous space
29 *
30 * Originally used for representing points in the hyperbolic plane.
31 * Currently used for all kinds of supported spaces, as well as
32 * for all vector spaces (up to 4 dimensions). We are using
33 * the normalized homogeneous coordinates, which allows us to work with most
34 * geometries in HyperRogue in a uniform way.
35
36 * In the hyperbolic plane, this is the Minkowski hyperboloid model:
37 * (x,y,z) such that x*x+y*y-z*z == -1 and z > 0.
38 *
39 * In spherical geometry, we have x*x+y*y+z*z == 1.
40 *
41 * In Euclidean geometry, we have z = 1.
42 *
43 * In isotropic 3D geometries an extra coordinate is added.
44 *
45 * In nonisotropic coordinates h[3] == 1.
46 *
47 * In product geometries the 'z' coordinate is modelled by multiplying all
48 * three coordinates with exp(z).
49 *
50 */
51
52 struct hyperpoint : array<ld, MAXMDIM> {
hyperpointhr::hyperpoint53 hyperpoint() {}
54
55 #if MAXMDIM == 4
hyperpointhr::hyperpoint56 constexpr hyperpoint(ld x, ld y, ld z, ld w) : array<ld, MAXMDIM> {{x,y,z,w}} {}
57 #else
hyperpointhr::hyperpoint58 constexpr hyperpoint(ld x, ld y, ld z, ld w) : array<ld, MAXMDIM> {{x,y,z}} {}
59 #endif
60
operator *=hr::hyperpoint61 inline hyperpoint& operator *= (ld d) {
62 for(int i=0; i<MXDIM; i++) self[i] *= d;
63 return self;
64 }
65
operator /=hr::hyperpoint66 inline hyperpoint& operator /= (ld d) {
67 for(int i=0; i<MXDIM; i++) self[i] /= d;
68 return self;
69 }
70
operator +=hr::hyperpoint71 inline hyperpoint& operator += (const hyperpoint h2) {
72 for(int i=0; i<MXDIM; i++) self[i] += h2[i];
73 return self;
74 }
75
operator -=hr::hyperpoint76 inline hyperpoint& operator -= (const hyperpoint h2) {
77 for(int i=0; i<MXDIM; i++) self[i] -= h2[i];
78 return self;
79 }
80
operator *(ld d,hyperpoint h)81 inline friend hyperpoint operator * (ld d, hyperpoint h) { return h *= d; }
operator *(hyperpoint h,ld d)82 inline friend hyperpoint operator * (hyperpoint h, ld d) { return h *= d; }
operator /(hyperpoint h,ld d)83 inline friend hyperpoint operator / (hyperpoint h, ld d) { return h /= d; }
operator +(hyperpoint h,hyperpoint h2)84 inline friend hyperpoint operator + (hyperpoint h, hyperpoint h2) { return h += h2; }
operator -(hyperpoint h,hyperpoint h2)85 inline friend hyperpoint operator - (hyperpoint h, hyperpoint h2) { return h -= h2; }
86
operator -(hyperpoint h)87 inline friend hyperpoint operator - (hyperpoint h) { return h * -1; }
88
89 // cross product
operator ^(hyperpoint h1,hyperpoint h2)90 inline friend hyperpoint operator ^ (hyperpoint h1, hyperpoint h2) {
91 return hyperpoint(
92 h1[1] * h2[2] - h1[2] * h2[1],
93 h1[2] * h2[0] - h1[0] * h2[2],
94 h1[0] * h2[1] - h1[1] * h2[0],
95 0
96 );
97 }
98
dot_d(int c,hyperpoint h1,hyperpoint h2)99 friend ld dot_d(int c, hyperpoint h1, hyperpoint h2) {
100 ld sum = 0;
101 for(int i=0; i<c; i++) sum += h1[i] * h2[i];
102 return sum;
103 }
104
105 // Euclidean inner product
operator |(hyperpoint h1,hyperpoint h2)106 inline friend ld operator | (hyperpoint h1, hyperpoint h2) {
107 return dot_d(MXDIM, h1, h2);
108 }
109 };
110
111 /** \brief A matrix acting on hr::hyperpoint
112 *
113 * Since we are using homogeneous coordinates for hr::hyperpoint,
114 * rotations and translations can be represented
115 * as matrix multiplications. Other applications of matrices in HyperRogue
116 * (in dimension up to 4) are also implemented using transmatrix.
117 */
118 struct transmatrix {
119 ld tab[MAXMDIM][MAXMDIM];
operator []hr::transmatrix120 hyperpoint& operator [] (int i) { return (hyperpoint&)tab[i][0]; }
operator []hr::transmatrix121 const hyperpoint& operator [] (int i) const { return (const hyperpoint&)tab[i]; }
122
operator *(const transmatrix & T,const hyperpoint & H)123 inline friend hyperpoint operator * (const transmatrix& T, const hyperpoint& H) {
124 hyperpoint z;
125 for(int i=0; i<MXDIM; i++) {
126 z[i] = 0;
127 for(int j=0; j<MXDIM; j++) z[i] += T[i][j] * H[j];
128 }
129 return z;
130 }
131
operator *(const transmatrix & T,const transmatrix & U)132 inline friend transmatrix operator * (const transmatrix& T, const transmatrix& U) {
133 transmatrix R;
134 for(int i=0; i<MXDIM; i++) for(int j=0; j<MXDIM; j++) {
135 R[i][j] = 0;
136 for(int k=0; k<MXDIM; k++)
137 R[i][j] += T[i][k] * U[k][j];
138 }
139 return R;
140 }
141 };
142
143 /** @brief hyperpoint with shift
144 * shift has two uses:
145 * (1) in the 'universal cover of SL' geometry, shift is used for the extra angular coordinate
146 * (2) in band models, shift is used to draw faraway points correctly
147 */
148 struct shiftpoint {
149 hyperpoint h;
150 ld shift;
operator []hr::shiftpoint151 ld& operator [] (int i) { return h[i]; }
operator []hr::shiftpoint152 const ld& operator [] (int i) const { return h[i]; }
operator +(const shiftpoint & h,const hyperpoint & h2)153 inline friend shiftpoint operator + (const shiftpoint& h, const hyperpoint& h2) {
154 return shiftpoint{h.h+h2, h.shift};
155 }
operator -(const shiftpoint & h,const hyperpoint & h2)156 inline friend shiftpoint operator - (const shiftpoint& h, const hyperpoint& h2) {
157 return shiftpoint{h.h-h2, h.shift};
158 }
159 };
160
shiftless(const hyperpoint & h,ld shift=0)161 inline shiftpoint shiftless(const hyperpoint& h, ld shift = 0) {
162 shiftpoint res; res.h = h; res.shift = shift; return res;
163 }
164
165 struct shiftmatrix {
166 transmatrix T;
167 ld shift;
operator []hr::shiftmatrix168 hyperpoint& operator [] (int i) { return T[i]; }
operator []hr::shiftmatrix169 const hyperpoint& operator [] (int i) const { return T[i]; }
operator *(const shiftmatrix & T,const hyperpoint & h)170 inline friend shiftpoint operator * (const shiftmatrix& T, const hyperpoint& h) {
171 return shiftpoint{T.T*h, T.shift};
172 }
operator *(const shiftmatrix & T,const transmatrix & U)173 inline friend shiftmatrix operator * (const shiftmatrix& T, const transmatrix& U) {
174 return shiftmatrix{T.T*U, T.shift};
175 }
176 };
177
shiftless(const transmatrix & T,ld shift=0)178 inline shiftmatrix shiftless(const transmatrix& T, ld shift = 0) {
179 shiftmatrix res; res.T = T; res.shift = shift; return res;
180 }
181
182 /** returns a diagonal matrix */
diag(ld a,ld b,ld c,ld d)183 constexpr transmatrix diag(ld a, ld b, ld c, ld d) {
184 #if MAXMDIM==3
185 return transmatrix{{{a,0,0}, {0,b,0}, {0,0,c}}};
186 #else
187 return transmatrix{{{a,0,0,0}, {0,b,0,0}, {0,0,c,0}, {0,0,0,d}}};
188 #endif
189 }
190
191 constexpr hyperpoint Hypc = hyperpoint(0, 0, 0, 0);
192
193 /** identity matrix */
194 constexpr transmatrix Id = diag(1,1,1,1);
195
196 /** zero matrix */
197 constexpr transmatrix Zero = diag(0,0,0,0);
198
199 /** mirror image */
200 constexpr transmatrix Mirror = diag(1,-1,1,1);
201
202 /** mirror image: flip in the Y coordinate */
203 constexpr transmatrix MirrorY = diag(1,-1,1,1);
204
205 /** mirror image: flip in the X coordinate */
206 constexpr transmatrix MirrorX = diag(-1,1,1,1);
207
208 /** mirror image: flip in the Z coordinate */
209 constexpr transmatrix MirrorZ = diag(1,1,-1,1);
210
211 /** rotate by PI in the XY plane */
212 constexpr transmatrix pispin = diag(-1,-1,1,1);
213
214 /** central symmetry matrix */
215 constexpr transmatrix centralsym = diag(-1,-1,-1,-1);
216
hpxyz(ld x,ld y,ld z)217 inline hyperpoint hpxyz(ld x, ld y, ld z) { return MDIM == 3 ? hyperpoint(x,y,z,0) : hyperpoint(x,y,0,z); }
hpxyz3(ld x,ld y,ld z,ld w)218 inline hyperpoint hpxyz3(ld x, ld y, ld z, ld w) { return MDIM == 3 ? hyperpoint(x,y,w,0) : hyperpoint(x,y,z,w); }
point3(ld x,ld y,ld z)219 constexpr hyperpoint point3(ld x, ld y, ld z) { return hyperpoint(x,y,z,0); }
point30(ld x,ld y,ld z)220 constexpr hyperpoint point30(ld x, ld y, ld z) { return hyperpoint(x,y,z,0); }
point31(ld x,ld y,ld z)221 constexpr hyperpoint point31(ld x, ld y, ld z) { return hyperpoint(x,y,z,1); }
point2(ld x,ld y)222 constexpr hyperpoint point2(ld x, ld y) { return hyperpoint(x,y,0,0); }
223
224 constexpr hyperpoint C02 = hyperpoint(0,0,1,0);
225 constexpr hyperpoint C03 = hyperpoint(0,0,0,1);
226
227 /** C0 is the origin in our space */
228 #define C0 (MDIM == 3 ? C02 : C03)
229 #endif
230
231 // basic functions and types
232 //===========================
233
234 #ifndef M_PI
235 #define M_PI 3.14159265358979
236 #endif
237
squar(ld x)238 EX ld squar(ld x) { return x*x; }
239
sig(int z)240 EX int sig(int z) { return ginf[geometry].g.sig[z]; }
241
curvature()242 EX int curvature() {
243 switch(cgclass) {
244 case gcEuclid: return 0;
245 case gcHyperbolic: return -1;
246 case gcSphere: return 1;
247 case gcProduct: return PIU(curvature());
248 default: return 0;
249 }
250 }
251
sin_auto(ld x)252 EX ld sin_auto(ld x) {
253 switch(cgclass) {
254 case gcEuclid: return x;
255 case gcHyperbolic: return sinh(x);
256 case gcSphere: return sin(x);
257 case gcProduct: return PIU(sin_auto(x));
258 case gcSL2: return sinh(x);
259 default: return x;
260 }
261 }
262
asin_auto(ld x)263 EX ld asin_auto(ld x) {
264 switch(cgclass) {
265 case gcEuclid: return x;
266 case gcHyperbolic: return asinh(x);
267 case gcSphere: return asin(x);
268 case gcProduct: return PIU(asin_auto(x));
269 case gcSL2: return asinh(x);
270 default: return x;
271 }
272 }
273
acos_auto(ld x)274 EX ld acos_auto(ld x) {
275 switch(cgclass) {
276 case gcHyperbolic: return acosh(x);
277 case gcSphere: return acos(x);
278 case gcProduct: return PIU(acos_auto(x));
279 case gcSL2: return acosh(x);
280 default: return x;
281 }
282 }
283
284 /** \brief volume of a three-dimensional ball of radius r in the current isotropic geometry */
volume_auto(ld r)285 EX ld volume_auto(ld r) {
286 switch(cgclass) {
287 case gcEuclid: return 4 * r * r * r / 3 * M_PI;
288 case gcHyperbolic: return M_PI * (sinh(2*r) - 2 * r);
289 case gcSphere: return M_PI * (2 * r - sin(2*r));
290 default: return 0;
291 }
292 }
293
294 /** \brief area of a circle of radius r in the current isotropic geometry */
area_auto(ld r)295 EX ld area_auto(ld r) {
296 switch(cgclass) {
297 case gcEuclid: return r * r * M_PI;
298 case gcHyperbolic: return 2 * M_PI * (cosh(r) - 1);
299 case gcSphere: return 2 * M_PI * (1 - cos(r));
300 default: return 0;
301 }
302 }
303
304 /** \brief volume in 3D, area in 2D */
wvolarea_auto(ld r)305 EX ld wvolarea_auto(ld r) {
306 if(WDIM == 3) return volume_auto(r);
307 else return area_auto(r);
308 }
309
asin_clamp(ld x)310 EX ld asin_clamp(ld x) { return x>1 ? M_PI/2 : x<-1 ? -M_PI/2 : std::isnan(x) ? 0 : asin(x); }
311
acos_clamp(ld x)312 EX ld acos_clamp(ld x) { return x>1 ? 0 : x<-1 ? M_PI : std::isnan(x) ? 0 : acos(x); }
313
asin_auto_clamp(ld x)314 EX ld asin_auto_clamp(ld x) {
315 switch(cgclass) {
316 case gcEuclid: return x;
317 case gcHyperbolic: return asinh(x);
318 case gcSL2: return asinh(x);
319 case gcSphere: return asin_clamp(x);
320 default: return x;
321 }
322 }
323
acos_auto_clamp(ld x)324 EX ld acos_auto_clamp(ld x) {
325 switch(cgclass) {
326 case gcHyperbolic: return x < 1 ? 0 : acosh(x);
327 case gcSL2: return x < 1 ? 0 : acosh(x);
328 case gcSphere: return acos_clamp(x);
329 case gcProduct: return PIU(acos_auto_clamp(x));
330 default: return x;
331 }
332 }
333
cos_auto(ld x)334 EX ld cos_auto(ld x) {
335 switch(cgclass) {
336 case gcEuclid: return 1;
337 case gcHyperbolic: return cosh(x);
338 case gcSL2: return cosh(x);
339 case gcSphere: return cos(x);
340 case gcProduct: return PIU(cos_auto(x));
341 default: return 1;
342 }
343 }
344
tan_auto(ld x)345 EX ld tan_auto(ld x) {
346 switch(cgclass) {
347 case gcEuclid: return x;
348 case gcHyperbolic: return tanh(x);
349 case gcSphere: return tan(x);
350 case gcProduct: return PIU(tan_auto(x));
351 case gcSL2: return tanh(x);
352 default: return 1;
353 }
354 }
355
atan_auto(ld x)356 EX ld atan_auto(ld x) {
357 switch(cgclass) {
358 case gcEuclid: return x;
359 case gcHyperbolic: return atanh(x);
360 case gcSphere: return atan(x);
361 case gcProduct: return PIU(atan_auto(x));
362 case gcSL2: return atanh(x);
363 default: return x;
364 }
365 }
366
atan2_auto(ld y,ld x)367 EX ld atan2_auto(ld y, ld x) {
368 switch(cgclass) {
369 case gcEuclid: return y/x;
370 case gcHyperbolic: return atanh(y/x);
371 case gcSL2: return atanh(y/x);
372 case gcSphere: return atan2(y, x);
373 case gcProduct: return PIU(atan2_auto(y, x));
374 default: return y/x;
375 }
376 }
377
378 /** This function returns the length of the edge opposite the angle alpha in
379 * a triangle with angles alpha, beta, gamma. This is called the cosine rule,
380 * and of course works only in non-Euclidean geometry. */
edge_of_triangle_with_angles(ld alpha,ld beta,ld gamma)381 EX ld edge_of_triangle_with_angles(ld alpha, ld beta, ld gamma) {
382 return acos_auto((cos(alpha) + cos(beta) * cos(gamma)) / (sin(beta) * sin(gamma)));
383 }
384
hpxy(ld x,ld y)385 EX hyperpoint hpxy(ld x, ld y) {
386 if(sl2) return hyperpoint(x, y, 0, sqrt(1+x*x+y*y));
387 if(rotspace) return hyperpoint(x, y, 0, sqrt(1-x*x-y*y));
388 return PIU(hpxyz(x,y, translatable ? 1 : sphere ? sqrt(1-x*x-y*y) : sqrt(1+x*x+y*y)));
389 }
390
hpxy3(ld x,ld y,ld z)391 EX hyperpoint hpxy3(ld x, ld y, ld z) {
392 return hpxyz3(x,y,z, sl2 ? sqrt(1+x*x+y*y-z*z) :translatable ? 1 : sphere ? sqrt(1-x*x-y*y-z*z) : sqrt(1+x*x+y*y+z*z));
393 }
394
395 #if HDR
396 // a point (I hope this number needs no comments ;) )
397 constexpr hyperpoint Cx12 = hyperpoint(1,0,1.41421356237,0);
398 constexpr hyperpoint Cx13 = hyperpoint(1,0,0,1.41421356237);
399
400 #define Cx1 (GDIM==2?Cx12:Cx13)
401 #endif
402
zero_d(int d,hyperpoint h)403 EX bool zero_d(int d, hyperpoint h) {
404 for(int i=0; i<d; i++) if(h[i]) return false;
405 return true;
406 }
407
408 /** this function returns approximate square of distance between two points
409 * (in the spherical analogy, this would be the distance in the 3D space,
410 * through the interior, not on the surface)
411 * also used to verify whether a point h1 is on the hyperbolic plane by using Hypc for h2
412 */
413
intval(const hyperpoint & h1,const hyperpoint & h2)414 EX ld intval(const hyperpoint &h1, const hyperpoint &h2) {
415 ld res = 0;
416 for(int i=0; i<MDIM; i++) res += squar(h1[i] - h2[i]) * sig(i);
417 if(elliptic) {
418 ld res2 = 0;
419 for(int i=0; i<MDIM; i++) res2 += squar(h1[i] + h2[i]) * sig(i);
420 return min(res, res2);
421 }
422 return res;
423 }
424
quickdist(const hyperpoint & h1,const hyperpoint & h2)425 EX ld quickdist(const hyperpoint &h1, const hyperpoint &h2) {
426 if(prod) return hdist(h1, h2);
427 return intval(h1, h2);
428 }
429
430 /** square Euclidean hypotenuse in the first d dimensions */
sqhypot_d(int d,const hyperpoint & h)431 EX ld sqhypot_d(int d, const hyperpoint& h) {
432 ld sum = 0;
433 for(int i=0; i<d; i++) sum += h[i]*h[i];
434 return sum;
435 }
436
437 /** Euclidean hypotenuse in the first d dimensions */
hypot_d(int d,const hyperpoint & h)438 EX ld hypot_d(int d, const hyperpoint& h) {
439 return sqrt(sqhypot_d(d, h));
440 }
441
442 /** @brief h1 and h2 define a line; to_other_side(h1,h2)*x is x moved orthogonally to this line, by double the distance from C0
443 * (I suppose it could be done better)
444 */
to_other_side(hyperpoint h1,hyperpoint h2)445 EX transmatrix to_other_side(hyperpoint h1, hyperpoint h2) {
446
447 ld d = hdist(h1, h2);
448
449 hyperpoint v;
450 if(euclid)
451 v = (h2 - h1) / d;
452 else
453 v = (h1 * cos_auto(d) - h2) / sin_auto(d);
454
455 ld d1;
456 if(euclid)
457 d1 = -(v|h1) / (v|v);
458 else
459 d1 = atan_auto(-v[LDIM] / h1[LDIM]);
460
461 hyperpoint hm = h1 * cos_auto(d1) + (sphere ? -1 : 1) * v * sin_auto(d1);
462
463 return rspintox(hm) * xpush(-hdist0(hm) * 2) * spintox(hm);
464 }
465
466 /** @brief positive for a material vertex, 0 for ideal vertex, negative for ultra-ideal vertex */
material(const hyperpoint & h)467 EX ld material(const hyperpoint& h) {
468 if(sphere) return intval(h, Hypc);
469 else if(hyperbolic) return -intval(h, Hypc);
470 else if(sl2) return h[2]*h[2] + h[3]*h[3] - h[0]*h[0] - h[1]*h[1];
471 else return h[LDIM];
472 }
473
zlevel(const hyperpoint & h)474 EX ld zlevel(const hyperpoint &h) {
475 if(sl2) return sqrt(-intval(h, Hypc));
476 else if(translatable) return h[LDIM];
477 else if(sphere) return sqrt(intval(h, Hypc));
478 else if(in_e2xe()) return log(h[2]);
479 else if(prod) return log(sqrt(abs(intval(h, Hypc)))); /* abs works with both underlying spherical and hyperbolic */
480 else return (h[LDIM] < 0 ? -1 : 1) * sqrt(-intval(h, Hypc));
481 }
482
hypot_auto(ld x,ld y)483 EX ld hypot_auto(ld x, ld y) {
484 switch(cgclass) {
485 case gcEuclid:
486 return hypot(x, y);
487 case gcHyperbolic:
488 return acosh(cosh(x) * cosh(y));
489 case gcSphere:
490 return acos(cos(x) * cos(y));
491 default:
492 return hypot(x, y);
493 }
494 }
495
496 /** normalize the homogeneous coordinates */
normalize(hyperpoint H)497 EX hyperpoint normalize(hyperpoint H) {
498 if(prod) return H;
499 ld Z = zlevel(H);
500 for(int c=0; c<MXDIM; c++) H[c] /= Z;
501 return H;
502 }
503
504 /** like normalize but makes (ultra)ideal points material */
ultra_normalize(hyperpoint H)505 EX hyperpoint ultra_normalize(hyperpoint H) {
506 if(material(H) <= 0) {
507 H[LDIM] = hypot_d(LDIM, H) + 1e-6;
508 }
509 return normalize(H);
510 }
511
512 /** normalize, and in product geometry, also flatten */
normalize_flat(hyperpoint h)513 EX hyperpoint normalize_flat(hyperpoint h) {
514 if(prod) return product_decompose(h).second;
515 if(sl2) h = slr::translate(h) * zpush0(-atan2(h[2], h[3]));
516 return normalize(h);
517 }
518
519 /** get the center of the line segment from H1 to H2 */
mid(const hyperpoint & H1,const hyperpoint & H2)520 EX hyperpoint mid(const hyperpoint& H1, const hyperpoint& H2) {
521 if(prod) {
522 auto d1 = product_decompose(H1);
523 auto d2 = product_decompose(H2);
524 return zshift(PIU( mid(d1.second, d2.second) ), (d1.first + d2.first) / 2);
525 }
526 return normalize(H1 + H2);
527 }
528
mid(const shiftpoint & H1,const shiftpoint & H2)529 EX shiftpoint mid(const shiftpoint& H1, const shiftpoint& H2) {
530 return shiftless(mid(H1.h, H2.h), (H1.shift + H2.shift)/2);
531 }
532
533 /** like mid, but take 3D into account */
midz(const hyperpoint & H1,const hyperpoint & H2)534 EX hyperpoint midz(const hyperpoint& H1, const hyperpoint& H2) {
535 if(prod) return mid(H1, H2);
536 hyperpoint H3 = H1 + H2;
537
538 ld Z = 2;
539
540 if(!euclid) Z = zlevel(H3) * 2 / (zlevel(H1) + zlevel(H2));
541 for(int c=0; c<MXDIM; c++) H3[c] /= Z;
542
543 return H3;
544 }
545
546 // matrices
547 //==========
548
549 /** rotate by alpha degrees in the coordinates a, b */
cspin(int a,int b,ld alpha)550 EX transmatrix cspin(int a, int b, ld alpha) {
551 transmatrix T = Id;
552 T[a][a] = +cos(alpha); T[a][b] = +sin(alpha);
553 T[b][a] = -sin(alpha); T[b][b] = +cos(alpha);
554 return T;
555 }
556
557 /** rotate by alpha degrees in the XY plane */
spin(ld alpha)558 EX transmatrix spin(ld alpha) { return cspin(0, 1, alpha); }
559
random_spin3()560 EX transmatrix random_spin3() {
561 ld alpha2 = asin(randd() * 2 - 1);
562 ld alpha = randd() * 2 * M_PI;
563 ld alpha3 = randd() * 2 * M_PI;
564 return cspin(0, 1, alpha) * cspin(0, 2, alpha2) * cspin(1, 2, alpha3);
565 }
566
random_spin()567 EX transmatrix random_spin() {
568 if(WDIM == 2) return spin(randd() * 2 * M_PI);
569 else return random_spin3();
570 }
571
eupush(ld x,ld y)572 EX transmatrix eupush(ld x, ld y) {
573 transmatrix T = Id;
574 T[0][LDIM] = x;
575 T[1][LDIM] = y;
576 return T;
577 }
578
euclidean_translate(ld x,ld y,ld z)579 EX transmatrix euclidean_translate(ld x, ld y, ld z) {
580 transmatrix T = Id;
581 T[0][LDIM] = x;
582 T[1][LDIM] = y;
583 T[2][LDIM] = z;
584 return T;
585 }
586
euscale(ld x,ld y)587 EX transmatrix euscale(ld x, ld y) {
588 transmatrix T = Id;
589 T[0][0] = x;
590 T[1][1] = y;
591 return T;
592 }
593
euscale3(ld x,ld y,ld z)594 EX transmatrix euscale3(ld x, ld y, ld z) {
595 transmatrix T = Id;
596 T[0][0] = x;
597 T[1][1] = y;
598 T[2][2] = z;
599 return T;
600 }
601
602 EX transmatrix eupush(hyperpoint h, ld co IS(1)) {
603 if(nonisotropic) return nisot::translate(h, co);
604 transmatrix T = Id;
605 for(int i=0; i<GDIM; i++) T[i][LDIM] = h[i] * co;
606 return T;
607 }
608
eupush3(ld x,ld y,ld z)609 EX transmatrix eupush3(ld x, ld y, ld z) {
610 if(sl2) return slr::translate(slr::xyz_point(x, y, z));
611 return eupush(point3(x, y, z));
612 }
613
euscalezoom(hyperpoint h)614 EX transmatrix euscalezoom(hyperpoint h) {
615 transmatrix T = Id;
616 T[0][0] = h[0];
617 T[0][1] = -h[1];
618 T[1][0] = h[1];
619 T[1][1] = h[0];
620 return T;
621 }
622
euaffine(hyperpoint h)623 EX transmatrix euaffine(hyperpoint h) {
624 transmatrix T = Id;
625 T[0][1] = h[0];
626 T[1][1] = exp(h[1]);
627 return T;
628 }
629
cpush(int cid,ld alpha)630 EX transmatrix cpush(int cid, ld alpha) {
631 transmatrix T = Id;
632 if(prod && cid == 2)
633 return mscale(Id, alpha);
634 if(nonisotropic)
635 return eupush3(cid == 0 ? alpha : 0, cid == 1 ? alpha : 0, cid == 2 ? alpha : 0);
636 T[LDIM][LDIM] = T[cid][cid] = cos_auto(alpha);
637 T[cid][LDIM] = sin_auto(alpha);
638 T[LDIM][cid] = -curvature() * sin_auto(alpha);
639 return T;
640 }
641
cmirror(int cid)642 EX transmatrix cmirror(int cid) {
643 transmatrix T = Id;
644 T[cid][cid] = -1;
645 return T;
646 }
647
648 // push alpha units to the right
xpush(ld alpha)649 EX transmatrix xpush(ld alpha) { return cpush(0, alpha); }
650
651 EX bool eqmatrix(transmatrix A, transmatrix B, ld eps IS(.01)) {
652 for(int i=0; i<MXDIM; i++)
653 for(int j=0; j<MXDIM; j++)
654 if(std::abs(A[i][j] - B[i][j]) > eps)
655 return false;
656 return true;
657 }
658
659 #if MAXMDIM >= 4
660 // in the 3D space, move the point h orthogonally to the (x,y) plane by z units
orthogonal_move(const hyperpoint & h,ld z)661 EX hyperpoint orthogonal_move(const hyperpoint& h, ld z) {
662 if(prod) return zshift(h, z);
663 if(sl2) return slr::translate(h) * cpush0(2, z);
664 if(!hyperbolic) return rgpushxto0(h) * cpush(2, z) * C0;
665 if(nil) return nisot::translate(h) * cpush0(2, z);
666 if(translatable) return hpxy3(h[0], h[1], h[2] + z);
667 ld u = 1;
668 if(h[2]) z += asin_auto(h[2]), u /= cos_auto(asin_auto(h[2]));
669 u *= cos_auto(z);
670 return hpxy3(h[0] * u, h[1] * u, sinh(z));
671 }
672 #endif
673
674 // push alpha units vertically
ypush(ld alpha)675 EX transmatrix ypush(ld alpha) { return cpush(1, alpha); }
676
zpush(ld z)677 EX transmatrix zpush(ld z) { return cpush(2, z); }
678
matrix3(ld a,ld b,ld c,ld d,ld e,ld f,ld g,ld h,ld i)679 EX transmatrix matrix3(ld a, ld b, ld c, ld d, ld e, ld f, ld g, ld h, ld i) {
680 #if MAXMDIM==3
681 return transmatrix {{{a,b,c},{d,e,f},{g,h,i}}};
682 #else
683 if(GDIM == 2)
684 return transmatrix {{{a,b,c,0},{d,e,f,0},{g,h,i,0},{0,0,0,1}}};
685 else
686 return transmatrix {{{a,b,0,c},{d,e,0,f},{0,0,1,0},{g,h,0,i}}};
687 #endif
688 }
689
matrix4(ld a,ld b,ld c,ld d,ld e,ld f,ld g,ld h,ld i,ld j,ld k,ld l,ld m,ld n,ld o,ld p)690 EX transmatrix matrix4(ld a, ld b, ld c, ld d, ld e, ld f, ld g, ld h, ld i, ld j, ld k, ld l, ld m, ld n, ld o, ld p) {
691 #if MAXMDIM==3
692 return transmatrix {{{a,b,d},{e,f,h},{m,n,p}}};
693 #else
694 return transmatrix {{{a,b,c,d},{e,f,g,h},{i,j,k,l},{m,n,o,p}}};
695 #endif
696 }
697
698 #if MAXMDIM >= 4
swapmatrix(transmatrix & T)699 EX void swapmatrix(transmatrix& T) {
700 for(int i=0; i<4; i++) swap(T[i][2], T[i][3]);
701 for(int i=0; i<4; i++) swap(T[2][i], T[3][i]);
702 if(GDIM == 3) {
703 for(int i=0; i<4; i++) T[i][2] = T[2][i] = 0;
704 T[2][2] = 1;
705 }
706 fixmatrix(T);
707 for(int i=0; i<4; i++) for(int j=0; j<4; j++) if(isnan(T[i][j])) T = Id;
708 }
709
swapmatrix(hyperpoint & h)710 EX void swapmatrix(hyperpoint& h) {
711 swap(h[2], h[3]);
712 }
713 #endif
714
parabolic1(ld u)715 EX transmatrix parabolic1(ld u) {
716 if(euclid)
717 return ypush(u);
718 else {
719 ld diag = u*u/2;
720 return matrix3(
721 -diag+1, u, diag,
722 -u, 1, u,
723 -diag, u, diag+1
724 );
725 }
726 }
727
parabolic13(ld u,ld v)728 EX transmatrix parabolic13(ld u, ld v) {
729 if(euclid)
730 return ypush(u);
731 else {
732 ld diag = (u*u+v*v)/2;
733 return matrix4(
734 -diag+1, u, v, diag,
735 -u, 1, 0, u,
736 -v, 0, 1, v,
737 -diag, u, v, diag+1
738 );
739 }
740 }
741
parabolic10(hyperpoint h)742 EX hyperpoint parabolic10(hyperpoint h) {
743 if(euclid) { h[LDIM] = 1; return h; }
744 else if(MDIM == 4) return hyperpoint(sinh(h[0]), h[1]/exp(h[0]), h[2]/exp(h[0]), cosh(h[0]));
745 else return hyperpoint(sinh(h[0]), h[1]/exp(h[0]), cosh(h[0]), 0);
746 }
747
deparabolic10(const hyperpoint h)748 EX hyperpoint deparabolic10(const hyperpoint h) {
749 if(euclid) return h;
750 ld x = -log(h[LDIM] - h[0]);
751 if(MDIM == 3) return hyperpoint(x, h[1] * exp(x), 1, 0);
752 return point31(x, h[1] * exp(x), h[2] * exp(x));
753 }
754
spintoc(const hyperpoint & H,int t,int f)755 EX transmatrix spintoc(const hyperpoint& H, int t, int f) {
756 transmatrix T = Id;
757 ld R = hypot(H[f], H[t]);
758 if(R >= 1e-12) {
759 T[t][t] = +H[t]/R; T[t][f] = +H[f]/R;
760 T[f][t] = -H[f]/R; T[f][f] = +H[t]/R;
761 }
762 return T;
763 }
764
765 /** an Euclidean rotation in the axes (t,f) which rotates
766 * the point H to the positive 't' axis
767 */
768
rspintoc(const hyperpoint & H,int t,int f)769 EX transmatrix rspintoc(const hyperpoint& H, int t, int f) {
770 transmatrix T = Id;
771 ld R = hypot(H[f], H[t]);
772 if(R >= 1e-12) {
773 T[t][t] = +H[t]/R; T[t][f] = -H[f]/R;
774 T[f][t] = +H[f]/R; T[f][f] = +H[t]/R;
775 }
776 return T;
777 }
778
779 /** an isometry which takes the point H to the positive X axis
780 * \see rspintox
781 */
spintox(const hyperpoint & H)782 EX transmatrix spintox(const hyperpoint& H) {
783 if(GDIM == 2 || prod) return spintoc(H, 0, 1);
784 transmatrix T1 = spintoc(H, 0, 1);
785 return spintoc(T1*H, 0, 2) * T1;
786 }
787
788 /** inverse of hr::spintox
789 */
rspintox(const hyperpoint & H)790 EX transmatrix rspintox(const hyperpoint& H) {
791 if(GDIM == 2 || prod) return rspintoc(H, 0, 1);
792 transmatrix T1 = spintoc(H, 0, 1);
793 return rspintoc(H, 0, 1) * rspintoc(T1*H, 0, 2);
794 }
795
796 /** for H on the X axis, this matrix pushes H to C0
797 * \see gpushxto0
798 */
799
pushxto0(const hyperpoint & H)800 EX transmatrix pushxto0(const hyperpoint& H) {
801 transmatrix T = Id;
802 T[0][0] = +H[LDIM]; T[0][LDIM] = -H[0];
803 T[LDIM][0] = curvature() * H[0]; T[LDIM][LDIM] = +H[LDIM];
804 return T;
805 }
806
807 /** set the i-th column of T to H */
set_column(transmatrix & T,int i,const hyperpoint & H)808 EX void set_column(transmatrix& T, int i, const hyperpoint& H) {
809 for(int j=0; j<MXDIM; j++)
810 T[j][i] = H[j];
811 }
812
get_column(transmatrix & T,int i)813 EX hyperpoint get_column(transmatrix& T, int i) {
814 hyperpoint h;
815 for(int j=0; j<MXDIM; j++)
816 h[j] = T[j][i];
817 return h;
818 }
819
820 /** build a matrix using the given vectors as columns */
build_matrix(hyperpoint h1,hyperpoint h2,hyperpoint h3,hyperpoint h4)821 EX transmatrix build_matrix(hyperpoint h1, hyperpoint h2, hyperpoint h3, hyperpoint h4) {
822 transmatrix T;
823 for(int i=0; i<MXDIM; i++) {
824 T[i][0] = h1[i],
825 T[i][1] = h2[i],
826 T[i][2] = h3[i];
827 if(MAXMDIM == 4) T[i][3] = h4[i];
828 }
829 return T;
830 }
831
832 /** for H on the X axis, this matrix pushes C0 to H
833 * \see rgpushxto0
834 */
835
rpushxto0(const hyperpoint & H)836 EX transmatrix rpushxto0(const hyperpoint& H) {
837 transmatrix T = Id;
838 T[0][0] = +H[LDIM]; T[0][LDIM] = H[0];
839 T[LDIM][0] = -curvature() * H[0]; T[LDIM][LDIM] = +H[LDIM];
840 return T;
841 }
842
ggpushxto0(const hyperpoint & H,ld co)843 EX transmatrix ggpushxto0(const hyperpoint& H, ld co) {
844 if(translatable)
845 return eupush(H, co);
846 if(prod) {
847 auto d = product_decompose(H);
848 return mscale(PIU(ggpushxto0(d.second, co)), d.first * co);
849 }
850 transmatrix res = Id;
851 if(sqhypot_d(GDIM, H) < 1e-12) return res;
852 ld fac = -curvature()/(H[LDIM]+1);
853 for(int i=0; i<GDIM; i++)
854 for(int j=0; j<GDIM; j++)
855 res[i][j] += H[i] * H[j] * fac;
856
857 for(int d=0; d<GDIM; d++)
858 res[d][LDIM] = co * H[d],
859 res[LDIM][d] = -curvature() * co * H[d];
860 res[LDIM][LDIM] = H[LDIM];
861
862 return res;
863 }
864
865 /** a translation matrix which takes H to 0 */
gpushxto0(const hyperpoint & H)866 EX transmatrix gpushxto0(const hyperpoint& H) {
867 return ggpushxto0(H, -1);
868 }
869
870 /** a translation matrix which takes 0 to H */
rgpushxto0(const hyperpoint & H)871 EX transmatrix rgpushxto0(const hyperpoint& H) {
872 return ggpushxto0(H, 1);
873 }
874
rgpushxto0(const shiftpoint & H)875 EX shiftmatrix rgpushxto0(const shiftpoint& H) {
876 return shiftless(rgpushxto0(H.h), H.shift);
877 }
878
879 /** \brief Fix the numerical inaccuracies in the isometry T
880 *
881 * The nature of hyperbolic geometry makes the computations numerically unstable.
882 * The numerical errors tend to accumulate, eventually destroying the projection.
883 * This function fixes this problem by replacing T with a 'correct' isometry.
884 */
885
fixmatrix(transmatrix & T)886 EX void fixmatrix(transmatrix& T) {
887 if(nonisotropic) ; // T may be inverse... do not do that
888 else if(cgflags & qAFFINE) ; // affine
889 else if(prod) {
890 auto z = zlevel(tC0(T));
891 T = mscale(T, -z);
892 PIU(fixmatrix(T));
893 T = mscale(T, +z);
894 }
895 else if(euclid)
896 fixmatrix_euclid(T);
897 else
898 orthonormalize(T);
899 }
900
fixmatrix_euclid(transmatrix & T)901 EX void fixmatrix_euclid(transmatrix& T) {
902 for(int x=0; x<GDIM; x++) for(int y=0; y<=x; y++) {
903 ld dp = 0;
904 for(int z=0; z<GDIM; z++) dp += T[z][x] * T[z][y];
905
906 if(y == x) dp = 1 - sqrt(1/dp);
907
908 for(int z=0; z<GDIM; z++) T[z][x] -= dp * T[z][y];
909 }
910 for(int x=0; x<GDIM; x++) T[LDIM][x] = 0;
911 T[LDIM][LDIM] = 1;
912 }
913
orthonormalize(transmatrix & T)914 EX void orthonormalize(transmatrix& T) {
915 for(int x=0; x<MDIM; x++) for(int y=0; y<=x; y++) {
916 ld dp = 0;
917 for(int z=0; z<MXDIM; z++) dp += T[z][x] * T[z][y] * sig(z);
918
919 if(y == x) dp = 1 - sqrt(sig(x)/dp);
920
921 for(int z=0; z<MXDIM; z++) T[z][x] -= dp * T[z][y];
922 }
923 }
924
925 /** fix a 3D rotation matrix */
fix_rotation(transmatrix & rot)926 EX void fix_rotation(transmatrix& rot) {
927 dynamicval<eGeometry> g(geometry, gSphere);
928 fixmatrix(rot);
929 for(int i=0; i<3; i++) rot[i][3] = rot[3][i] = 0;
930 rot[3][3] = 1;
931 }
932
933 /** determinant 2x2 */
det2(const transmatrix & T)934 EX ld det2(const transmatrix& T) {
935 return T[0][0] * T[1][1] - T[0][1] * T[1][0];
936 }
937
938 /** determinant 3x3 */
det3(const transmatrix & T)939 EX ld det3(const transmatrix& T) {
940 ld det = 0;
941 for(int i=0; i<3; i++)
942 det += T[0][i] * T[1][(i+1)%3] * T[2][(i+2)%3];
943 for(int i=0; i<3; i++)
944 det -= T[0][i] * T[1][(i+2)%3] * T[2][(i+1)%3];
945 return det;
946 }
947
948 /** determinant */
det(const transmatrix & T)949 EX ld det(const transmatrix& T) {
950 if(MDIM == 3)
951 return det3(T);
952 else {
953 ld det = 1;
954 transmatrix M = T;
955 for(int a=0; a<MDIM; a++) {
956 int max_at = a;
957 for(int b=a; b<MDIM; b++)
958 if(abs(M[b][a]) > abs(M[max_at][a]))
959 max_at = b;
960
961 if(max_at != a)
962 for(int c=a; c<MDIM; c++) tie(M[max_at][c], M[a][c]) = make_pair(-M[a][c], M[max_at][c]);
963
964 if(!M[a][a]) return 0;
965 for(int b=a+1; b<MDIM; b++) {
966 ld co = -M[b][a] / M[a][a];
967 for(int c=a; c<MDIM; c++) M[b][c] += M[a][c] * co;
968 }
969 det *= M[a][a];
970 }
971 return det;
972 }
973 }
974
975 /** warning about incorrect inverse */
inverse_error(const transmatrix & T)976 void inverse_error(const transmatrix& T) {
977 println(hlog, "Warning: inverting a singular matrix: ", T);
978 }
979
980 /** inverse of a 3x3 matrix */
inverse3(const transmatrix & T)981 EX transmatrix inverse3(const transmatrix& T) {
982 ld d = det(T);
983 transmatrix T2;
984 if(d == 0) {
985 inverse_error(T);
986 return Id;
987 }
988
989 for(int i=0; i<3; i++)
990 for(int j=0; j<3; j++)
991 T2[j][i] = (T[(i+1)%3][(j+1)%3] * T[(i+2)%3][(j+2)%3] - T[(i+1)%3][(j+2)%3] * T[(i+2)%3][(j+1)%3]) / d;
992 return T2;
993 }
994
995 /** \brief inverse of a general matrix */
inverse(const transmatrix & T)996 EX transmatrix inverse(const transmatrix& T) {
997 if(MDIM == 3)
998 return inverse3(T);
999 else {
1000 transmatrix T1 = T;
1001 transmatrix T2 = Id;
1002
1003 for(int a=0; a<MDIM; a++) {
1004 int best = a;
1005
1006 for(int b=a+1; b<MDIM; b++)
1007 if(abs(T1[b][a]) > abs(T1[best][a]))
1008 best = b;
1009
1010 int b = best;
1011
1012 if(b != a)
1013 for(int c=0; c<MDIM; c++)
1014 swap(T1[b][c], T1[a][c]), swap(T2[b][c], T2[a][c]);
1015
1016 if(!T1[a][a]) { inverse_error(T); return Id; }
1017 for(int b=a+1; b<=GDIM; b++) {
1018 ld co = -T1[b][a] / T1[a][a];
1019 for(int c=0; c<MDIM; c++) T1[b][c] += T1[a][c] * co, T2[b][c] += T2[a][c] * co;
1020 }
1021 }
1022
1023 for(int a=MDIM-1; a>=0; a--) {
1024 for(int b=0; b<a; b++) {
1025 ld co = -T1[b][a] / T1[a][a];
1026 for(int c=0; c<MDIM; c++) T1[b][c] += T1[a][c] * co, T2[b][c] += T2[a][c] * co;
1027 }
1028 ld co = 1 / T1[a][a];
1029 for(int c=0; c<MDIM; c++) T1[a][c] *= co, T2[a][c] *= co;
1030 }
1031 return T2;
1032 }
1033 }
1034
1035 /** \brief inverse of an orthogonal matrix, i.e., transposition */
ortho_inverse(transmatrix T)1036 EX transmatrix ortho_inverse(transmatrix T) {
1037 for(int i=1; i<MDIM; i++)
1038 for(int j=0; j<i; j++)
1039 swap(T[i][j], T[j][i]);
1040 return T;
1041 }
1042
1043 /** \brief inverse of an orthogonal matrix in Minkowski space */
pseudo_ortho_inverse(transmatrix T)1044 EX transmatrix pseudo_ortho_inverse(transmatrix T) {
1045 for(int i=1; i<MXDIM; i++)
1046 for(int j=0; j<i; j++)
1047 swap(T[i][j], T[j][i]);
1048 for(int i=0; i<MDIM-1; i++)
1049 T[i][MDIM-1] = -T[i][MDIM-1],
1050 T[MDIM-1][i] = -T[MDIM-1][i];
1051 return T;
1052 }
1053
1054 /** \brief inverse of an isometry -- in most geometries this can be done more efficiently than using inverse */
iso_inverse(const transmatrix & T)1055 EX transmatrix iso_inverse(const transmatrix& T) {
1056 if(hyperbolic)
1057 return pseudo_ortho_inverse(T);
1058 if(sphere)
1059 return ortho_inverse(T);
1060 if(nil) {
1061 transmatrix U = Id;
1062 U[2][LDIM] = T[0][LDIM] * T[1][LDIM] - T[2][LDIM];
1063 U[1][LDIM] = -T[1][LDIM];
1064 U[2][1] = U[0][LDIM] = -T[0][LDIM];
1065 return U;
1066 }
1067 if(euclid && !(cgflags & qAFFINE)) {
1068 transmatrix U = Id;
1069 for(int i=0; i<MDIM-1; i++)
1070 for(int j=0; j<MDIM-1; j++)
1071 U[i][j] = T[j][i];
1072 hyperpoint h = U * tC0(T);
1073 for(int i=0; i<MDIM-1; i++)
1074 U[i][MDIM-1] = -h[i];
1075 return U;
1076 }
1077 return inverse(T);
1078 }
1079
1080 /** \brief T inverse a matrix T = O*S, where O is isometry and S is a scaling matrix (todo optimize) */
z_inverse(const transmatrix & T)1081 EX transmatrix z_inverse(const transmatrix& T) {
1082 return inverse(T);
1083 }
1084
1085 /** \brief T inverse a matrix T = O*P, where O is orthogonal and P is an isometry (todo optimize) */
view_inverse(transmatrix T)1086 EX transmatrix view_inverse(transmatrix T) {
1087 if(nonisotropic) return inverse(T);
1088 if(prod) return z_inverse(T);
1089 return iso_inverse(T);
1090 }
1091
1092 /** \brief T inverse a matrix T = P*O, where O is orthogonal and P is an isometry (todo optimize) */
iview_inverse(transmatrix T)1093 EX transmatrix iview_inverse(transmatrix T) {
1094 if(nonisotropic) return inverse(T);
1095 if(prod) return z_inverse(T);
1096 return iso_inverse(T);
1097 }
1098
product_decompose(hyperpoint h)1099 EX pair<ld, hyperpoint> product_decompose(hyperpoint h) {
1100 ld z = zlevel(h);
1101 return make_pair(z, mscale(h, -z));
1102 }
1103
1104 /** distance from mh and 0 */
hdist0(const hyperpoint & mh)1105 EX ld hdist0(const hyperpoint& mh) {
1106 switch(cgclass) {
1107 case gcHyperbolic:
1108 if(mh[LDIM] < 1) return 0;
1109 return acosh(mh[LDIM]);
1110 case gcEuclid: {
1111 return hypot_d(GDIM, mh);
1112 }
1113 case gcSphere: {
1114 ld res = mh[LDIM] >= 1 ? 0 : mh[LDIM] <= -1 ? M_PI : acos(mh[LDIM]);
1115 return res;
1116 }
1117 case gcProduct: {
1118 auto d1 = product_decompose(mh);
1119 return hypot(PIU(hdist0(d1.second)), d1.first);
1120 }
1121 #if MAXMDIM >= 4
1122 case gcSL2: {
1123 auto cosh_r = hypot(mh[2], mh[3]);
1124 auto phi = atan2(mh[2], mh[3]);
1125 return hypot(cosh_r < 1 ? 0 : acosh(cosh_r), phi);
1126 }
1127 case gcNil: {
1128 ld bz = mh[0] * mh[1] / 2;
1129 return hypot(mh[0], mh[1]) + abs(mh[2] - bz);
1130 }
1131 #endif
1132 default:
1133 return hypot_d(GDIM, mh);
1134 }
1135 }
1136
hdist0(const shiftpoint & mh)1137 EX ld hdist0(const shiftpoint& mh) {
1138 return hdist0(unshift(mh));
1139 }
1140
1141 /** length of a circle of radius r */
circlelength(ld r)1142 EX ld circlelength(ld r) {
1143 switch(cgclass) {
1144 case gcEuclid:
1145 return 2 * M_PI * r;
1146 case gcHyperbolic:
1147 return 2 * M_PI * sinh(r);
1148 case gcSphere:
1149 return 2 * M_PI * sin(r);
1150 default:
1151 return 2 * M_PI * r;
1152 }
1153 }
1154
1155 /* distance between h1 and h2 */
hdist(const hyperpoint & h1,const hyperpoint & h2)1156 EX ld hdist(const hyperpoint& h1, const hyperpoint& h2) {
1157 ld iv = intval(h1, h2);
1158 switch(cgclass) {
1159 case gcEuclid:
1160 if(iv < 0) return 0;
1161 return sqrt(iv);
1162 case gcHyperbolic:
1163 if(iv < 0) return 0;
1164 return 2 * asinh(sqrt(iv) / 2);
1165 case gcSphere:
1166 return 2 * asin_auto_clamp(sqrt(iv) / 2);
1167 case gcProduct: {
1168 auto d1 = product_decompose(h1);
1169 auto d2 = product_decompose(h2);
1170 return hypot(PIU(hdist(d1.second, d2.second)), d1.first - d2.first);
1171 }
1172 case gcSL2:
1173 return hdist0(stretch::itranslate(h1) * h2);
1174 default:
1175 if(iv < 0) return 0;
1176 return sqrt(iv);
1177 }
1178 }
1179
hdist(const shiftpoint & h1,const shiftpoint & h2)1180 EX ld hdist(const shiftpoint& h1, const shiftpoint& h2) {
1181 return hdist(h1.h, unshift(h2, h1.shift));
1182 }
1183
mscale(const hyperpoint & t,double fac)1184 EX hyperpoint mscale(const hyperpoint& t, double fac) {
1185 if(GDIM == 3 && !prod) return cpush(2, fac) * t;
1186 if(prod) fac = exp(fac);
1187 hyperpoint res;
1188 for(int i=0; i<MXDIM; i++)
1189 res[i] = t[i] * fac;
1190 return res;
1191 }
1192
mscale(const shiftpoint & t,double fac)1193 EX shiftpoint mscale(const shiftpoint& t, double fac) {
1194 return shiftless(mscale(t.h, fac), t.shift);
1195 }
1196
mscale(const transmatrix & t,double fac)1197 EX transmatrix mscale(const transmatrix& t, double fac) {
1198 if(GDIM == 3 && !prod) {
1199 // if(pmodel == mdFlatten) { transmatrix u = t; u[2][LDIM] -= fac; return u; }
1200 return t * cpush(2, fac);
1201 }
1202 if(prod) fac = exp(fac);
1203 transmatrix res;
1204 for(int i=0; i<MXDIM; i++) {
1205 for(int j=0; j<MDIM; j++)
1206 res[i][j] = t[i][j] * fac;
1207 for(int j=MDIM; j<MXDIM; j++)
1208 res[i][j] = t[i][j];
1209 }
1210 return res;
1211 }
1212
mscale(const shiftmatrix & t,double fac)1213 EX shiftmatrix mscale(const shiftmatrix& t, double fac) {
1214 return shiftless(mscale(t.T, fac), t.shift);
1215 }
1216
xyscale(const transmatrix & t,double fac)1217 EX transmatrix xyscale(const transmatrix& t, double fac) {
1218 transmatrix res;
1219 for(int i=0; i<MXDIM; i++) {
1220 for(int j=0; j<GDIM; j++)
1221 res[i][j] = t[i][j] * fac;
1222 for(int j=GDIM; j<MXDIM; j++)
1223 res[i][j] = t[i][j];
1224 }
1225 return res;
1226 }
1227
xyzscale(const transmatrix & t,double fac,double facz)1228 EX transmatrix xyzscale(const transmatrix& t, double fac, double facz) {
1229 transmatrix res;
1230 for(int i=0; i<MXDIM; i++) {
1231 for(int j=0; j<GDIM; j++)
1232 res[i][j] = t[i][j] * fac;
1233 res[i][LDIM] =
1234 t[i][LDIM] * facz;
1235 for(int j=LDIM+1; j<MXDIM; j++)
1236 res[i][j] = t[i][j];
1237 }
1238 return res;
1239 }
1240
xyzscale(const shiftmatrix & t,double fac,double facz)1241 EX shiftmatrix xyzscale(const shiftmatrix& t, double fac, double facz) {
1242 return shiftless(xyzscale(t.T, fac, facz), t.shift);
1243 }
1244
mzscale(const transmatrix & t,double fac)1245 EX transmatrix mzscale(const transmatrix& t, double fac) {
1246 if(GDIM == 3) return t * cpush(2, fac);
1247 // take only the spin
1248 transmatrix tcentered = gpushxto0(tC0(t)) * t;
1249 // tcentered = tcentered * spin(downspin_zivory);
1250 fac -= 1;
1251 transmatrix res = t * inverse(tcentered) * ypush(-fac) * tcentered;
1252 fac *= .2;
1253 fac += 1;
1254 for(int i=0; i<MXDIM; i++) for(int j=0; j<MXDIM; j++)
1255 res[i][j] = res[i][j] * fac;
1256 return res;
1257 }
1258
mzscale(const shiftmatrix & t,double fac)1259 EX shiftmatrix mzscale(const shiftmatrix& t, double fac) {
1260 return shiftless(mzscale(t.T, fac), t.shift);
1261 }
1262
mid3(hyperpoint h1,hyperpoint h2,hyperpoint h3)1263 EX hyperpoint mid3(hyperpoint h1, hyperpoint h2, hyperpoint h3) {
1264 return mid(h1+h2+h3, h1+h2+h3);
1265 }
1266
mid_at(hyperpoint h1,hyperpoint h2,ld v)1267 EX hyperpoint mid_at(hyperpoint h1, hyperpoint h2, ld v) {
1268 hyperpoint h = h1 * (1-v) + h2 * v;
1269 return mid(h, h);
1270 }
1271
mid_at_actual(hyperpoint h,ld v)1272 EX hyperpoint mid_at_actual(hyperpoint h, ld v) {
1273 return rspintox(h) * xpush0(hdist0(h) * v);
1274 }
1275
1276 /** in 3D, an orthogonal projection of C0 on the given triangle */
orthogonal_of_C0(hyperpoint h0,hyperpoint h1,hyperpoint h2)1277 EX hyperpoint orthogonal_of_C0(hyperpoint h0, hyperpoint h1, hyperpoint h2) {
1278 h0 /= h0[3];
1279 h1 /= h1[3];
1280 h2 /= h2[3];
1281 hyperpoint w = h0;
1282 hyperpoint d1 = h1 - h0;
1283 hyperpoint d2 = h2 - h0;
1284 ld denom = (d1|d1) * (d2|d2) - (d1|d2) * (d1|d2);
1285 ld a1 = (d2|w) * (d1|d2) - (d1|w) * (d2|d2);
1286 ld a2 = (d1|w) * (d1|d2) - (d2|w) * (d1|d1);
1287 hyperpoint h = w * denom + d1 * a1 + d2 * a2;
1288 return normalize(h);
1289 }
1290
zshift(hyperpoint x,ld z)1291 EX hyperpoint zshift(hyperpoint x, ld z) {
1292 if(GDIM == 3 && WDIM == 2) return rgpushxto0(x) * cpush0(2, z);
1293 else if(sl2) return mscale(x, z);
1294 else if(prod) return mscale(x, z);
1295 else return mscale(x, z);
1296 }
1297
hpxd(ld d,ld x,ld y,ld z)1298 EX hyperpoint hpxd(ld d, ld x, ld y, ld z) {
1299 hyperpoint H = hpxyz(d*x, d*y, z);
1300 H = mid(H, H);
1301 return H;
1302 }
1303
signum(ld x)1304 EX ld signum(ld x) { return x<0?-1:x>0?1:0; }
1305
asign(ld y1,ld y2)1306 EX bool asign(ld y1, ld y2) { return signum(y1) != signum(y2); }
1307
xcross(ld x1,ld y1,ld x2,ld y2)1308 EX ld xcross(ld x1, ld y1, ld x2, ld y2) { return x1 + (x2 - x1) * y1 / (y1 - y2); }
1309
parallel_transport(const transmatrix Position,const transmatrix & ori,const hyperpoint direction)1310 EX transmatrix parallel_transport(const transmatrix Position, const transmatrix& ori, const hyperpoint direction) {
1311 if(nonisotropic) return nisot::parallel_transport(Position, direction);
1312 else if(prod) {
1313 hyperpoint h = product::direct_exp(ori * direction);
1314 return Position * rgpushxto0(h);
1315 }
1316 else return Position * rgpushxto0(direct_exp(direction));
1317 }
1318
apply_parallel_transport(transmatrix & Position,const transmatrix orientation,const hyperpoint direction)1319 EX void apply_parallel_transport(transmatrix& Position, const transmatrix orientation, const hyperpoint direction) {
1320 Position = parallel_transport(Position, orientation, direction);
1321 }
1322
rotate_object(transmatrix & Position,transmatrix & orientation,transmatrix R)1323 EX void rotate_object(transmatrix& Position, transmatrix& orientation, transmatrix R) {
1324 if(prod) orientation = orientation * R;
1325 else Position = Position * R;
1326 }
1327
spin_towards(const transmatrix Position,transmatrix & ori,const hyperpoint goal,int dir,int back)1328 EX transmatrix spin_towards(const transmatrix Position, transmatrix& ori, const hyperpoint goal, int dir, int back) {
1329 transmatrix T;
1330 ld alpha = 0;
1331 if(nonisotropic && nisot::geodesic_movement)
1332 T = nisot::spin_towards(Position, goal);
1333 else {
1334 hyperpoint U = inverse(Position) * goal;
1335 if(prod) {
1336 hyperpoint h = product::inverse_exp(U);
1337 alpha = asin_clamp(h[2] / hypot_d(3, h));
1338 U = product_decompose(U).second;
1339 }
1340 T = rspintox(U);
1341 }
1342 if(back < 0) T = T * spin(M_PI), alpha = -alpha;
1343 if(prod) {
1344 if(dir == 0) ori = cspin(2, 0, alpha);
1345 if(dir == 2) ori = cspin(2, 0, alpha - M_PI/2), dir = 0;
1346 }
1347 if(dir) T = T * cspin(dir, 0, -M_PI/2);
1348 T = Position * T;
1349 return T;
1350 }
1351
spin_towards(const shiftmatrix Position,transmatrix & ori,const shiftpoint goal,int dir,int back)1352 EX shiftmatrix spin_towards(const shiftmatrix Position, transmatrix& ori, const shiftpoint goal, int dir, int back) {
1353 return shiftless(spin_towards(Position.T, ori, unshift(goal, Position.shift), dir, back), Position.shift);
1354 }
1355
ortho_error(transmatrix T)1356 EX ld ortho_error(transmatrix T) {
1357
1358 ld err = 0;
1359
1360 for(int x=0; x<3; x++) for(int y=0; y<3; y++) {
1361 ld s = 0;
1362 for(int z=0; z<3; z++) s += T[z][x] * T[z][y];
1363
1364 s -= (x==y);
1365 err += s*s;
1366 }
1367
1368 return err;
1369 }
1370
transpose(transmatrix T)1371 EX transmatrix transpose(transmatrix T) {
1372 transmatrix result;
1373 for(int i=0; i<MXDIM; i++)
1374 for(int j=0; j<MXDIM; j++)
1375 result[j][i] = T[i][j];
1376 return result;
1377 }
1378
1379 #if HDR
1380 namespace slr {
1381 hyperpoint xyz_point(ld x, ld y, ld z);
1382 hyperpoint polar(ld r, ld theta, ld phi);
1383 }
1384
cpush0(int c,ld x)1385 inline hyperpoint cpush0(int c, ld x) {
1386 hyperpoint h = Hypc;
1387 if(sl2) return slr::xyz_point(c==0?x:0, c==1?x:0, c==2?x:0);
1388 if(c == 2 && prod) {
1389 h[2] = exp(x);
1390 return h;
1391 }
1392 h[LDIM] = cos_auto(x);
1393 h[c] = sin_auto(x);
1394 return h;
1395 }
1396
xspinpush0(ld alpha,ld x)1397 inline hyperpoint xspinpush0(ld alpha, ld x) {
1398 if(sl2) return slr::polar(x, -alpha, 0);
1399 hyperpoint h = Hypc;
1400 h[LDIM] = cos_auto(x);
1401 h[0] = sin_auto(x) * cos(alpha);
1402 h[1] = sin_auto(x) * -sin(alpha);
1403 return h;
1404 }
1405
xpush0(ld x)1406 inline hyperpoint xpush0(ld x) { return cpush0(0, x); }
ypush0(ld x)1407 inline hyperpoint ypush0(ld x) { return cpush0(1, x); }
zpush0(ld x)1408 inline hyperpoint zpush0(ld x) { return cpush0(2, x); }
1409
1410 /** T * C0, optimized */
tC0(const transmatrix & T)1411 inline hyperpoint tC0(const transmatrix &T) {
1412 hyperpoint z;
1413 for(int i=0; i<MXDIM; i++) z[i] = T[i][LDIM];
1414 return z;
1415 }
1416
tC0_t(const transmatrix & T)1417 inline hyperpoint tC0_t(const transmatrix &T) { return tC0(T); }
1418
tC0(const shiftmatrix & T)1419 inline shiftpoint tC0(const shiftmatrix &T) {
1420 return shiftpoint{tC0(T.T), T.shift};
1421 }
1422 #endif
1423
1424 /** tangent vector in the given direction */
ctangent(int c,ld x)1425 EX hyperpoint ctangent(int c, ld x) { return point3(c==0?x:0, c==1?x:0, c==2?x:0); }
1426
1427 /** tangent vector in direction X */
xtangent(ld x)1428 EX hyperpoint xtangent(ld x) { return ctangent(0, x); }
1429
1430 /** tangent vector in direction Y */
ztangent(ld z)1431 EX hyperpoint ztangent(ld z) { return ctangent(2, z); }
1432
1433 /** change the length of the targent vector */
tangent_length(hyperpoint dir,ld length)1434 EX hyperpoint tangent_length(hyperpoint dir, ld length) {
1435 ld r = hypot_d(GDIM, dir);
1436 if(!r) return dir;
1437 return dir * (length / r);
1438 }
1439
1440 /** exponential function: follow the geodesic given by v */
direct_exp(hyperpoint v)1441 EX hyperpoint direct_exp(hyperpoint v) {
1442 #if CAP_SOLV
1443 if(sn::in()) return nisot::numerical_exp(v);
1444 #endif
1445 #if MAXMDIM >= 4
1446 if(nil) return nilv::formula_exp(v);
1447 if(sl2 || stretch::in()) return stretch::mstretch ? nisot::numerical_exp(v) : rots::formula_exp(v);
1448 #endif
1449 if(prod) return product::direct_exp(v);
1450 ld d = hypot_d(GDIM, v);
1451 if(d > 0) for(int i=0; i<GDIM; i++) v[i] = v[i] * sin_auto(d) / d;
1452 v[LDIM] = cos_auto(d);
1453 return v;
1454 }
1455
1456 #if HDR
1457 constexpr flagtype pfNO_INTERPOLATION = 1; /**< in tables (sol/nih geometries), do not use interpolations */
1458 constexpr flagtype pfNO_DISTANCE = 2; /**< we just need the directions -- this makes it a bit faster in sol/nih geometries */
1459 constexpr flagtype pfLOW_BS_ITER = 4; /**< low iterations in binary search (nil geometry, sl2 not affected currently) */
1460
1461 constexpr flagtype pQUICK = pfNO_INTERPOLATION | pfLOW_BS_ITER;
1462
1463 constexpr flagtype pNORMAL = 0;
1464 #endif
1465
1466 /** inverse exponential function \see hr::direct_exp */
inverse_exp(const shiftpoint h,flagtype prec IS (pNORMAL))1467 EX hyperpoint inverse_exp(const shiftpoint h, flagtype prec IS(pNORMAL)) {
1468 #if CAP_SOLV
1469 if(sn::in()) {
1470 if(nih)
1471 return sn::get_inverse_exp_nsym(h.h, prec);
1472 else
1473 return sn::get_inverse_exp_symsol(h.h, prec);
1474 }
1475 #endif
1476 if(nil) return nilv::get_inverse_exp(h.h, prec);
1477 if(sl2) return slr::get_inverse_exp(h);
1478 if(prod) return product::inverse_exp(h.h);
1479 ld d = acos_auto_clamp(h[GDIM]);
1480 hyperpoint v = Hypc;
1481 if(d && sin_auto(d)) for(int i=0; i<GDIM; i++) v[i] = h[i] * d / sin_auto(d);
1482 v[3] = 0;
1483 return v;
1484 }
1485
geo_dist(const hyperpoint h1,const hyperpoint h2,flagtype prec IS (pNORMAL))1486 EX ld geo_dist(const hyperpoint h1, const hyperpoint h2, flagtype prec IS(pNORMAL)) {
1487 if(!nonisotropic) return hdist(h1, h2);
1488 return hypot_d(3, inverse_exp(shiftless(nisot::translate(h1, -1) * h2, prec)));
1489 }
1490
geo_dist(const shiftpoint h1,const shiftpoint h2,flagtype prec IS (pNORMAL))1491 EX ld geo_dist(const shiftpoint h1, const shiftpoint h2, flagtype prec IS(pNORMAL)) {
1492 if(!nonisotropic) return hdist(h1, h2);
1493 return hypot_d(3, inverse_exp(shiftless(nisot::translate(h1.h, -1) * h2.h, h2.shift - h1.shift), prec));
1494 }
1495
geo_dist_q(const hyperpoint h1,const hyperpoint h2,flagtype prec IS (pNORMAL))1496 EX ld geo_dist_q(const hyperpoint h1, const hyperpoint h2, flagtype prec IS(pNORMAL)) {
1497 auto d = geo_dist(h1, h2, prec);
1498 if(elliptic && d > M_PI/2) return M_PI - d;
1499 return d;
1500 }
1501
lp_iapply(const hyperpoint h)1502 EX hyperpoint lp_iapply(const hyperpoint h) {
1503 return nisot::local_perspective_used() ? inverse(NLP) * h : h;
1504 }
1505
lp_apply(const hyperpoint h)1506 EX hyperpoint lp_apply(const hyperpoint h) {
1507 return nisot::local_perspective_used() ? NLP * h : h;
1508 }
1509
smalltangent()1510 EX hyperpoint smalltangent() { return xtangent(.1); }
1511
cyclefix(ld & a,ld b)1512 EX void cyclefix(ld& a, ld b) {
1513 while(a > b + M_PI) a -= 2 * M_PI;
1514 while(a < b - M_PI) a += 2 * M_PI;
1515 }
1516
raddif(ld a,ld b)1517 EX ld raddif(ld a, ld b) {
1518 ld d = a-b;
1519 if(d < 0) d = -d;
1520 if(d > 2*M_PI) d -= 2*M_PI;
1521 if(d > M_PI) d = 2 * M_PI-d;
1522 return d;
1523 }
1524
bucketer(ld x)1525 EX unsigned bucketer(ld x) {
1526 return unsigned((long long)(x * 10000 + 100000.5) - 100000);
1527 }
1528
bucketer(hyperpoint h)1529 EX unsigned bucketer(hyperpoint h) {
1530 unsigned dx = 0;
1531 if(prod) {
1532 auto d = product_decompose(h);
1533 h = d.second;
1534 dx += bucketer(d.first) * 50;
1535 }
1536 dx += bucketer(h[0]) + 1000 * bucketer(h[1]) + 1000000 * bucketer(h[2]);
1537 if(MDIM == 4) dx += bucketer(h[3]) * 1000000001;
1538 if(elliptic) dx = min(dx, -dx);
1539 return dx;
1540 }
1541
1542 #if MAXMDIM >= 4
1543 /** @brief project the origin to the triangle [h1,h2,h3] */
project_on_triangle(hyperpoint h1,hyperpoint h2,hyperpoint h3)1544 EX hyperpoint project_on_triangle(hyperpoint h1, hyperpoint h2, hyperpoint h3) {
1545 h1 /= h1[3];
1546 h2 /= h2[3];
1547 h3 /= h3[3];
1548 transmatrix T;
1549 T[0] = h1; T[1] = h2; T[2] = h3;
1550 T[3] = C0;
1551 ld det_orig = det(T);
1552 hyperpoint orthogonal = (h2 - h1) ^ (h3 - h1);
1553 T[0] = orthogonal; T[1] = h2-h1; T[2] = h3-h1;
1554 ld det_orth = det(T);
1555 hyperpoint result = orthogonal * (det_orig / det_orth);
1556 result[3] = 1;
1557 return normalize(result);
1558 }
1559 #endif
1560
lerp(hyperpoint a0,hyperpoint a1,ld x)1561 EX hyperpoint lerp(hyperpoint a0, hyperpoint a1, ld x) {
1562 return a0 + (a1-a0) * x;
1563 }
1564
linecross(hyperpoint a,hyperpoint b,hyperpoint c,hyperpoint d)1565 EX hyperpoint linecross(hyperpoint a, hyperpoint b, hyperpoint c, hyperpoint d) {
1566 a /= a[LDIM];
1567 b /= b[LDIM];
1568 c /= c[LDIM];
1569 d /= d[LDIM];
1570
1571 ld bax = b[0] - a[0];
1572 ld dcx = d[0] - c[0];
1573 ld cax = c[0] - a[0];
1574 ld bay = b[1] - a[1];
1575 ld dcy = d[1] - c[1];
1576 ld cay = c[1] - a[1];
1577
1578 hyperpoint res;
1579 res[0] = (cay * dcx * bax + a[0] * bay * dcx - c[0] * dcy * bax) / (bay * dcx - dcy * bax);
1580 res[1] = (cax * dcy * bay + a[1] * bax * dcy - c[1] * dcx * bay) / (bax * dcy - dcx * bay);
1581 res[2] = 0;
1582 res[3] = 0;
1583 res[GDIM] = 1;
1584 return normalize(res);
1585 }
1586
inner2(hyperpoint h1,hyperpoint h2)1587 EX ld inner2(hyperpoint h1, hyperpoint h2) {
1588 return
1589 hyperbolic ? h1[LDIM] * h2[LDIM] - h1[0] * h2[0] - h1[1] * h2[1] :
1590 sphere ? h1[LDIM] * h2[LDIM] + h1[0] * h2[0] + h1[1] * h2[1] :
1591 h1[0] * h2[0] + h1[1] * h2[1];
1592 }
1593
circumscribe(hyperpoint a,hyperpoint b,hyperpoint c)1594 EX hyperpoint circumscribe(hyperpoint a, hyperpoint b, hyperpoint c) {
1595 hyperpoint h = C0;
1596
1597 b = b - a;
1598 c = c - a;
1599
1600 if(euclid) {
1601 ld b2 = inner2(b, b)/2;
1602 ld c2 = inner2(c, c)/2;
1603
1604 ld det = c[1]*b[0] - b[1]*c[0];
1605
1606 h = a;
1607
1608 h[1] += (c2*b[0] - b2 * c[0]) / det;
1609 h[0] += (c2*b[1] - b2 * c[1]) / -det;
1610
1611 return h;
1612 }
1613
1614 if(inner2(b,b) < 0) {
1615 b = b / sqrt(-inner2(b, b));
1616 c = c + b * inner2(c, b);
1617 h = h + b * inner2(h, b);
1618 }
1619 else {
1620 b = b / sqrt(inner2(b, b));
1621 c = c - b * inner2(c, b);
1622 h = h - b * inner2(h, b);
1623 }
1624
1625 if(inner2(c,c) < 0) {
1626 c = c / sqrt(-inner2(c, c));
1627 h = h + c * inner2(h, c);
1628 }
1629 else {
1630 c = c / sqrt(inner2(c, c));
1631 h = h - c * inner2(h, c);
1632 }
1633
1634 if(h[LDIM] < 0) h[0] = -h[0], h[1] = -h[1], h[LDIM] = -h[LDIM];
1635
1636 ld i = inner2(h, h);
1637 if(i > 0) h /= sqrt(i);
1638 else h /= -sqrt(-i);
1639
1640 return h;
1641 }
1642
inner3(hyperpoint h1,hyperpoint h2)1643 EX ld inner3(hyperpoint h1, hyperpoint h2) {
1644 return
1645 hyperbolic ? h1[LDIM] * h2[LDIM] - h1[0] * h2[0] - h1[1] * h2[1] - h1[2]*h2[2]:
1646 sphere ? h1[LDIM] * h2[LDIM] + h1[0] * h2[0] + h1[1] * h2[1] + h1[2]*h2[2]:
1647 h1[0] * h2[0] + h1[1] * h2[1];
1648 }
1649
1650 /** circumscribe for H3 and S3 (not for E3 yet!) */
circumscribe(hyperpoint a,hyperpoint b,hyperpoint c,hyperpoint d)1651 EX hyperpoint circumscribe(hyperpoint a, hyperpoint b, hyperpoint c, hyperpoint d) {
1652
1653 array<hyperpoint, 4> ds = { b-a, c-a, d-a, C0 };
1654
1655 for(int i=0; i<3; i++) {
1656
1657 if(inner3(ds[i],ds[i]) < 0) {
1658 ds[i] = ds[i] / sqrt(-inner3(ds[i], ds[i]));
1659 for(int j=i+1; j<4; j++)
1660 ds[j] = ds[j] + ds[i] * inner3(ds[i], ds[j]);
1661 }
1662 else {
1663 ds[i] = ds[i] / sqrt(inner3(ds[i], ds[i]));
1664 for(int j=i+1; j<4; j++)
1665 ds[j] = ds[j] - ds[i] * inner3(ds[i], ds[j]);
1666 }
1667 }
1668
1669 hyperpoint& h = ds[3];
1670
1671 if(h[3] < 0) h = -h;
1672
1673 ld i = inner3(h, h);
1674 if(i > 0) h /= sqrt(i);
1675 else h /= -sqrt(-i);
1676
1677 return h;
1678 }
1679
clockwise(hyperpoint h1,hyperpoint h2)1680 EX bool clockwise(hyperpoint h1, hyperpoint h2) {
1681 return h1[0] * h2[1] > h1[1] * h2[0];
1682 }
1683
1684
1685 }
1686