From 429909bd561db5c77e921561d0da8e7c2cbdbf4a Mon Sep 17 00:00:00 2001 From: Zeno Rogue Date: Thu, 22 Aug 2019 13:45:51 +0200 Subject: [PATCH] improved documentation in hyperpoint --- hyperpoint.cpp | 179 ++++++++++++++++++++++++++++--------------------- 1 file changed, 101 insertions(+), 78 deletions(-) diff --git a/hyperpoint.cpp b/hyperpoint.cpp index b6faeb7c..f6c8dbd0 100644 --- a/hyperpoint.cpp +++ b/hyperpoint.cpp @@ -19,9 +19,30 @@ eGeometry geometry; eVariation variation; -// hyperbolic points and matrices - #if HDR +/** \brief A point in our continuous space + * Originally used for representing points in the hyperbolic plane. + * Currently used for all kinds of supported spaces, as well as + * for all vector spaces (up to 4 dimensions). We are using + * the normalized homogeneous coordinates, which allows us to work with most + * geometries in HyperRogue in a uniform way. + + * In the hyperbolic plane, this is the Minkowski hyperboloid model: + * (x,y,z) such that x*x+y*y-z*z == -1 and z > 0. + * + * In spherical geometry, we have x*x+y*y+z*z == 1. + * + * In Euclidean geometry, we have z = 1. + * + * In isotropic 3D geometries an extra coordinate is added. + * + * In nonisotropic coordinates h[3] == 1. + * + * In product geometries the 'z' coordinate is modelled by multiplying all + * three coordinates with exp(z). + * + */ + struct hyperpoint : array { hyperpoint() {} @@ -74,6 +95,12 @@ struct hyperpoint : array { } }; +/** \brief A matrix acting on hr::hyperpoint + * Since we are using homogeneous coordinates for hr::hyperpoint, + * rotations and translations can be represented + * as matrix multiplications. Other applications of matrices in HyperRogue + * (in dimension up to 4) are also implemented using transmatrix. + */ struct transmatrix { ld tab[MAXMDIM][MAXMDIM]; hyperpoint& operator [] (int i) { return (hyperpoint&)tab[i][0]; } @@ -99,7 +126,7 @@ struct transmatrix { } }; - +/** returns a diagonal matrix */ constexpr transmatrix diag(ld a, ld b, ld c, ld d) { #if MAXMDIM==3 return transmatrix{{{a,0,0}, {0,b,0}, {0,0,c}}}; @@ -110,26 +137,28 @@ constexpr transmatrix diag(ld a, ld b, ld c, ld d) { const static hyperpoint Hypc = hyperpoint(0, 0, 0, 0); -// identity matrix +/** identity matrix */ const static transmatrix Id = diag(1,1,1,1); -// zero matrix +/** zero matrix */ const static transmatrix Zero = diag(0,0,0,0); -// mirror image +/** mirror image */ const static transmatrix Mirror = diag(1,-1,1,1); + +/** mirror image: flip in the Y coordinate */ const static transmatrix MirrorY = diag(1,-1,1,1); -// mirror image +/** mirror image: flip in the X coordinate */ const static transmatrix MirrorX = diag(-1,1,1,1); -// mirror image +/** mirror image: flip in the Z coordinate */ const static transmatrix MirrorZ = diag(1,1,-1,1); -// rotate by PI +/** rotate by PI in the XY plane */ const static transmatrix pispin = diag(-1,-1,1,1); -// central symmetry +/** central symmetry matrix */ const static transmatrix centralsym = diag(-1,-1,-1,-1); inline hyperpoint hpxyz(ld x, ld y, ld z) { return MDIM == 3 ? hyperpoint(x,y,z,0) : hyperpoint(x,y,0,z); } @@ -140,31 +169,13 @@ inline hyperpoint point2(ld x, ld y) { return hyperpoint(x,y,0,0); } extern const hyperpoint C02, C03; +/** C0 is the origin in our space */ #define C0 (MDIM == 3 ? C02 : C03) #endif // basic functions and types //=========================== -#ifdef SINHCOSH -// ld sinh(ld alpha) { return (exp(alpha) - exp(-alpha)) / 2; } -// ld cosh(ld alpha) { return (exp(alpha) + exp(-alpha)) / 2; } - -/* ld inverse_sinh(ld z) { - return log(z+sqrt(1+z*z)); - } - -double inverse_cos(double c) { - double s = sqrt(1-c*c); - double r = atan(s/c); - if(r < 0) r = -r; - return r; - } - -// ld tanh(ld x) { return sinh(x) / cosh(x); } -ld inverse_tanh(ld x) { return log((1+x)/(1-x)) / 2; } */ - -#endif #ifndef M_PI #define M_PI 3.14159265358979 #endif @@ -279,18 +290,13 @@ EX ld atan2_auto(ld y, ld x) { } } -// cosine rule -- edge opposite alpha +/** This function returns the length of the edge opposite the angle alpha in + * a triangle with angles alpha, beta, gamma. This is called the cosine rule, + * and of course works only in non-Euclidean geometry. */ EX ld edge_of_triangle_with_angles(ld alpha, ld beta, ld gamma) { return acos_auto((cos(alpha) + cos(beta) * cos(gamma)) / (sin(beta) * sin(gamma))); } -// hyperbolic point: -//=================== - -// we represent the points on the hyperbolic plane -// by points in 3D space (Minkowski space) such that x^2+y^2-z^2 == -1, z > 0 -// (this is analogous to representing a sphere with points such that x^2+y^2+z^2 == 1) - hyperpoint hpxy(ld x, ld y) { return hpxyz(x,y, translatable ? 1 : sphere ? sqrt(1-x*x-y*y) : sqrt(1+x*x+y*y)); } @@ -314,10 +320,11 @@ EX bool zero_d(int d, hyperpoint h) { return true; } -// this function returns approximate square of distance between two points -// (in the spherical analogy, this would be the distance in the 3D space, -// through the interior, not on the surface) -// also used to verify whether a point h1 is on the hyperbolic plane by using Hypc for h2 +/** this function returns approximate square of distance between two points + * (in the spherical analogy, this would be the distance in the 3D space, + * through the interior, not on the surface) + * also used to verify whether a point h1 is on the hyperbolic plane by using Hypc for h2 + */ EX ld intval(const hyperpoint &h1, const hyperpoint &h2) { ld res = 0; @@ -335,12 +342,14 @@ EX ld quickdist(const hyperpoint &h1, const hyperpoint &h2) { return intval(h1, h2); } +/** square Euclidean hypotenuse in the first d dimensions */ EX ld sqhypot_d(int d, const hyperpoint& h) { ld sum = 0; for(int i=0; i= 0 +/** an isometry which takes the point H to the positive X axis + * \see rspintox + */ EX transmatrix spintox(const hyperpoint& H) { if(GDIM == 2 || prod) return spintoc(H, 0, 1); transmatrix T1 = spintoc(H, 0, 1); return spintoc(T1*H, 0, 2) * T1; } +/** inverse of hr::spintox + */ +EX transmatrix rspintox(const hyperpoint& H) { + if(GDIM == 2 || prod) return rspintoc(H, 0, 1); + transmatrix T1 = spintoc(H, 0, 1); + return rspintoc(H, 0, 1) * rspintoc(T1*H, 0, 2); + } + +/** for H on the X axis, this matrix pushes H to C0 + * \see gpushxto0 + */ + +EX transmatrix pushxto0(const hyperpoint& H) { + transmatrix T = Id; + T[0][0] = +H[LDIM]; T[0][LDIM] = -H[0]; + T[LDIM][0] = curvature() * H[0]; T[LDIM][LDIM] = +H[LDIM]; + return T; + } + +/** set the i-th column of T to H */ EX void set_column(transmatrix& T, int i, const hyperpoint& H) { for(int j=0; j product_decompose(hyperpoint h) { return make_pair(z, mscale(h, -z)); } -// distance between mh and 0 +/** distance from mh and 0 */ EX ld hdist0(const hyperpoint& mh) { switch(cgclass) { case gcHyperbolic: @@ -805,6 +832,7 @@ EX ld hdist0(const hyperpoint& mh) { } } +/** length of a circle of radius r */ EX ld circlelength(ld r) { switch(cgclass) { case gcEuclid: @@ -818,9 +846,8 @@ EX ld circlelength(ld r) { } } -// distance between two points +/* distance between h1 and h2 */ EX ld hdist(const hyperpoint& h1, const hyperpoint& h2) { - // return hdist0(gpushxto0(h1) * h2); ld iv = intval(h1, h2); switch(cgclass) { case gcEuclid: @@ -879,8 +906,6 @@ EX transmatrix xyzscale(const transmatrix& t, double fac, double facz) { return res; } -// double downspin_zivory; - EX transmatrix mzscale(const transmatrix& t, double fac) { if(GDIM == 3) return t * cpush(2, fac); // take only the spin @@ -897,8 +922,6 @@ EX transmatrix mzscale(const transmatrix& t, double fac) { EX transmatrix pushone() { return xpush(sphere?.5 : 1); } -// rotation matrix in R^3 - EX hyperpoint mid3(hyperpoint h1, hyperpoint h2, hyperpoint h3) { return mid(h1+h2+h3, h1+h2+h3); } @@ -912,7 +935,7 @@ EX hyperpoint mid_at_actual(hyperpoint h, ld v) { return rspintox(h) * xpush0(hdist0(h) * v); } -// in 3D, an orthogonal projection of C0 on the given triangle +/** in 3D, an orthogonal projection of C0 on the given triangle */ EX hyperpoint orthogonal_of_C0(hyperpoint h0, hyperpoint h1, hyperpoint h2) { h0 /= h0[3]; h1 /= h1[3]; @@ -1017,7 +1040,7 @@ inline hyperpoint xpush0(ld x) { return cpush0(0, x); } inline hyperpoint ypush0(ld x) { return cpush0(1, x); } inline hyperpoint zpush0(ld x) { return cpush0(2, x); } -// T * C0, optimized +/** T * C0, optimized */ inline hyperpoint tC0(const transmatrix &T) { hyperpoint z; for(int i=0; i