// Hyperbolic Rogue // This file contains hyperbolic points and matrices. // Copyright (C) 2011-2018 Zeno Rogue, see 'hyper.cpp' for details namespace hr { #if HDR static const ld degree = M_PI / 180; #endif eGeometry geometry; eVariation variation; // hyperbolic points and matrices // 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 EX ld squar(ld x) { return x*x; } EX int sig(int z) { return (sphere || sol || z1 ? M_PI/2 : x<-1 ? -M_PI/2 : std::isnan(x) ? 0 : asin(x); default: return x; } } EX ld cos_auto(ld x) { switch(cgclass) { case gcEuclid: return 1; case gcHyperbolic: return cosh(x); case gcSphere: return cos(x); default: return 1; } } EX ld tan_auto(ld x) { switch(cgclass) { case gcEuclid: return x; case gcHyperbolic: return tanh(x); case gcSphere: return tan(x); default: return 1; } } EX ld atan_auto(ld x) { switch(cgclass) { case gcEuclid: return x; case gcHyperbolic: return atanh(x); case gcSphere: return atan(x); default: return x; } } EX ld atan2_auto(ld y, ld x) { switch(cgclass) { case gcEuclid: return y/x; case gcHyperbolic: return atanh(y/x); case gcSphere: return atan2(y, x); default: return y/x; } } // cosine rule -- edge opposite alpha 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)); } hyperpoint hpxy3(ld x, ld y, ld z) { return hpxyz3(x,y,z, translatable ? 1 : sphere ? sqrt(1-x*x-y*y-z*z) : sqrt(1+x*x+y*y+z*z)); } // origin of the hyperbolic plane const hyperpoint C02 = hyperpoint(0,0,1,0); const hyperpoint C03 = hyperpoint(0,0,0,1); // a point (I hope this number needs no comments ;) ) const hyperpoint Cx12 = hyperpoint(1,0,1.41421356237,0); const hyperpoint Cx13 = hyperpoint(1,0,0,1.41421356237); #define Cx1 (GDIM==2?Cx12:Cx13) EX bool zero_d(int d, hyperpoint h) { for(int i=0; i eps) return false; return true; } #if MAXMDIM >= 4 // in the 3D space, move the point h orthogonally to the (x,y) plane by z units EX hyperpoint orthogonal_move(const hyperpoint& h, ld z) { if(!hyperbolic) return rgpushxto0(h) * cpush(2, z) * C0; if(nil) return nisot::translate(h) * cpush0(2, z); if(translatable) return hpxy3(h[0], h[1], h[2] + z); ld u = 1; if(h[2]) z += asin_auto(h[2]), u /= acos_auto(z); u *= cos_auto(z); return hpxy3(h[0] * u, h[1] * u, sinh(z)); } #endif // push alpha units vertically EX transmatrix ypush(ld alpha) { return cpush(1, alpha); } EX transmatrix zpush(ld z) { return cpush(2, z); } EX transmatrix matrix3(ld a, ld b, ld c, ld d, ld e, ld f, ld g, ld h, ld i) { #if MAXMDIM==3 return transmatrix {{{a,b,c},{d,e,f},{g,h,i}}}; #else if(DIM == 2) return transmatrix {{{a,b,c,0},{d,e,f,0},{g,h,i,0},{0,0,0,1}}}; else return transmatrix {{{a,b,0,c},{d,e,0,f},{0,0,1,0},{g,h,0,i}}}; #endif } 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) { #if MAXMDIM==3 return transmatrix {{{a,b,d},{e,f,h},{m,n,p}}}; #else return transmatrix {{{a,b,c,d},{e,f,g,h},{i,j,k,l},{m,n,o,p}}}; #endif } #if MAXMDIM >= 4 EX void swapmatrix(transmatrix& T) { for(int i=0; i<4; i++) swap(T[i][2], T[i][3]); for(int i=0; i<4; i++) swap(T[2][i], T[3][i]); if(DIM == 3) { for(int i=0; i<4; i++) T[i][2] = T[2][i] = 0; T[2][2] = 1; } fixmatrix(T); for(int i=0; i<4; i++) for(int j=0; j<4; j++) if(isnan(T[i][j])) T = Id; } EX void swapmatrix(hyperpoint& h) { swap(h[2], h[3]); } #endif EX transmatrix parabolic1(ld u) { if(euclid) return ypush(u); else { ld diag = u*u/2; return matrix3( -diag+1, u, diag, -u, 1, u, -diag, u, diag+1 ); } } EX transmatrix parabolic13(ld u, ld v) { if(euclid) return ypush(u); else { ld diag = (u*u+v*v)/2; return matrix4( -diag+1, u, v, diag, -u, 1, 0, u, -v, 0, 1, v, -diag, u, v, diag+1 ); } } EX transmatrix spintoc(const hyperpoint& H, int t, int f) { transmatrix T = Id; ld R = hypot(H[f], H[t]); if(R >= 1e-12) { T[t][t] = +H[t]/R; T[t][f] = +H[f]/R; T[f][t] = -H[f]/R; T[f][f] = +H[t]/R; } return T; } EX transmatrix rspintoc(const hyperpoint& H, int t, int f) { transmatrix T = Id; ld R = hypot(H[f], H[t]); if(R >= 1e-12) { T[t][t] = +H[t]/R; T[t][f] = -H[f]/R; T[f][t] = +H[f]/R; T[f][f] = +H[t]/R; } return T; } // rotate the hyperbolic plane around C0 such that H[1] == 0 and H[0] >= 0 EX transmatrix spintox(const hyperpoint& H) { if(GDIM == 2) return spintoc(H, 0, 1); transmatrix T1 = spintoc(H, 0, 1); return spintoc(T1*H, 0, 2) * T1; } EX void set_column(transmatrix& T, int i, const hyperpoint& H) { for(int j=0; j abs(T1[best][a])) best = b; int b = best; if(b != a) for(int c=0; c=0; a--) { for(int b=0; b= 1 ? 0 : mh[GDIM] <= -1 ? M_PI : acos(mh[GDIM]); if(elliptic && res > M_PI/2) res = M_PI-res; return res; } default: return hypot_d(GDIM, mh); } } EX ld circlelength(ld r) { switch(cgclass) { case gcEuclid: return 2 * M_PI * r; case gcHyperbolic: return 2 * M_PI * sinh(r); case gcSphere: return 2 * M_PI * sin(r); default: return 2 * M_PI * r; } } // distance between two points EX ld hdist(const hyperpoint& h1, const hyperpoint& h2) { // return hdist0(gpushxto0(h1) * h2); ld iv = intval(h1, h2); switch(cgclass) { case gcEuclid: if(iv < 0) return 0; return sqrt(iv); case gcHyperbolic: if(iv < 0) return 0; return 2 * asinh(sqrt(iv) / 2); case gcSphere: return 2 * asin_auto_clamp(sqrt(iv) / 2); default: if(iv < 0) return 0; return sqrt(iv); } } EX hyperpoint mscale(const hyperpoint& t, double fac) { if(GDIM == 3) return cpush(2, fac) * t; hyperpoint res; for(int i=0; i0?1:0; } bool asign(ld y1, ld y2) { return signum(y1) != signum(y2); } ld xcross(ld x1, ld y1, ld x2, ld y2) { return x1 + (x2 - x1) * y1 / (y1 - y2); } EX transmatrix solmul(const transmatrix T, const transmatrix V) { if(nonisotropic) return nisot::transport_view(T, V); else return T * V; } EX transmatrix solmul_pt(const transmatrix Position, const transmatrix T) { if(nonisotropic) return nisot::parallel_transport(Position, T); else return Position * T; } EX transmatrix spin_towards(const transmatrix Position, const hyperpoint goal, int dir, int back) { transmatrix T = nonisotropic ? nisot::spin_towards(Position, goal) : rspintox(inverse(Position) * goal); if(back < 0) T = T * spin(M_PI); if(dir) T = T * cspin(dir, 0, -M_PI/2); T = Position * T; return T; } EX ld ortho_error(transmatrix T) { ld err = 0; for(int x=0; x<3; x++) for(int y=0; y<3; y++) { ld s = 0; for(int z=0; z<3; z++) s += T[z][x] * T[z][y]; s -= (x==y); err += s*s; } return err; } }