// Hyperbolic Rogue -- basic computations in non-Euclidean geometry // Copyright (C) 2011-2019 Zeno Rogue, see 'hyper.cpp' for details /** \file hyperpoint.cpp * \brief basic computations in non-Euclidean geometry * * This implements hyperpoint (a point in non-Euclidean space), transmatrix (a transformation matrix), * and various basic routines related to them: rotations, translations, inverses and determinants, etc. * For nonisotropic geometries, it rather refers to nonisotropic.cpp. */ #include "hyper.h" namespace hr { #if HDR #ifndef M_PI #define M_PI 3.14159265358979 #endif static constexpr ld A_PI = M_PI; static constexpr ld TAU = 2 * A_PI; static constexpr ld degree = A_PI / 180; static const ld golden_phi = (sqrt(5)+1)/2; static const ld log_golden_phi = log(golden_phi); constexpr ld operator"" _deg(long double deg) { return deg * A_PI / 180; } #endif eGeometry geometry; eVariation variation; #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() {} #if MAXMDIM == 4 constexpr hyperpoint(ld x, ld y, ld z, ld w) : array {{x,y,z,w}} {} #else constexpr hyperpoint(ld x, ld y, ld z, ld w) : array {{x,y,z}} {} #endif inline hyperpoint& operator *= (ld d) { for(int i=0; i1 ? 90._deg : x<-1 ? -90._deg : std::isnan(x) ? 0 : asin(x); } EX ld acos_clamp(ld x) { return x>1 ? 0 : x<-1 ? M_PI : std::isnan(x) ? 0 : acos(x); } EX ld asin_auto_clamp(ld x) { switch(cgclass) { case gcEuclid: return x; case gcHyperbolic: return asinh(x); case gcSL2: return asinh(x); case gcSphere: return asin_clamp(x); case gcProduct: return PIU(asin_auto_clamp(x)); default: return x; } } EX ld acos_auto_clamp(ld x) { switch(cgclass) { case gcHyperbolic: return x < 1 ? 0 : acosh(x); case gcSL2: return x < 1 ? 0 : acosh(x); case gcSphere: return acos_clamp(x); case gcProduct: return PIU(acos_auto_clamp(x)); default: return x; } } EX ld cos_auto(ld x) { switch(cgclass) { case gcEuclid: return 1; case gcHyperbolic: return cosh(x); case gcSL2: return cosh(x); case gcSphere: return cos(x); case gcProduct: return PIU(cos_auto(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); case gcProduct: return PIU(tan_auto(x)); case gcSL2: return tanh(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); case gcProduct: return PIU(atan_auto(x)); case gcSL2: return atanh(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 gcSL2: return atanh(y/x); case gcSphere: return atan2(y, x); case gcProduct: return PIU(atan2_auto(y, x)); default: return y/x; } } /** 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))); } EX hyperpoint hpxy(ld x, ld y) { if(embedded_plane) { geom3::light_flip(true); hyperpoint h = hpxy(x, y); geom3::light_flip(false); return cgi.emb->base_to_actual(h); } if(sl2) return hyperpoint(x, y, 0, sqrt(1+x*x+y*y)); if(rotspace) return hyperpoint(x, y, 0, sqrt(1-x*x-y*y)); return PIU(hpxyz(x,y, translatable ? 1 : sphere ? sqrt(1-x*x-y*y) : sqrt(1+x*x+y*y))); } EX hyperpoint hpxy3(ld x, ld y, ld z) { 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)); } #if HDR // a point (I hope this number needs no comments ;) ) constexpr hyperpoint Cx12 = hyperpoint(1,0,1.41421356237,0); constexpr hyperpoint Cx13 = hyperpoint(1,0,0,1.41421356237); #define Cx1 (GDIM==2?Cx12:Cx13) #endif EX bool zero_d(int d, hyperpoint h) { for(int i=0; iis_sph_in_low() && !geom3::flipped) { geom3::light_flip(true); h1 = normalize(h1); h2 = normalize(h2); transmatrix T = to_other_side(h1, h2); fix4(T); geom3::light_flip(false); return T; } if(sol && meuclid) { /* works in 4x4... */ return gpushxto0(h1) * gpushxto0(h2); } ld d = hdist(h1, h2); hyperpoint v; if(euclid) v = (h2 - h1) / d; else v = (h1 * cos_auto(d) - h2) / sin_auto(d); ld d1; if(euclid) d1 = -(v|h1) / (v|v); else d1 = atan_auto(-v[LDIM] / h1[LDIM]); hyperpoint hm = h1 * cos_auto(d1) + (sphere ? -1 : 1) * v * sin_auto(d1); return rspintox(hm) * xpush(-hdist0(hm) * 2) * spintox(hm); } /** @brief positive for a material vertex, 0 for ideal vertex, negative for ultra-ideal vertex */ EX ld material(const hyperpoint& h) { if(sphere || in_s2xe()) return intval(h, Hypc); else if(hyperbolic || in_h2xe()) return -intval(h, Hypc); #if MAXMDIM >= 4 else if(sl2) return h[2]*h[2] + h[3]*h[3] - h[0]*h[0] - h[1]*h[1]; #endif else return h[LDIM]; } EX int safe_classify_ideals(hyperpoint h) { if(hyperbolic || in_h2xe()) { h /= h[LDIM]; ld x = MDIM == 3 ? 1 - (h[0] * h[0] + h[1] * h[1]) : 1 - (h[0] * h[0] + h[1] * h[1] + h[2] * h[2]); if(x > 1e-6) return 1; if(x < -1e-6) return -1; return 0; } return 1; } EX ld ideal_limit = 10; EX ld ideal_each = degree; EX hyperpoint safe_approximation_of_ideal(hyperpoint h) { return towards_inf(C0, h, ideal_limit); } /** the point on the line ab which is closest to zero. Might not be normalized. Works even if a and b are (ultra)ideal */ EX hyperpoint closest_to_zero(hyperpoint a, hyperpoint b) { if(sqhypot_d(MDIM, a-b) < 1e-9) return a; if(isnan(a[0])) return a; a /= a[LDIM]; b /= b[LDIM]; ld mul_a = 0, mul_b = 0; for(int i=0; i eps) return false; return true; } // 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(GDIM == 2 || MDIM == 3) 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 } EX transmatrix parabolic1(ld u) { if(euclid) return ypush(u); else if(cgi.emb->is_hyp_in_solnih() && !geom3::flipped) { 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 eupush3(0, u, v); else if(cgi.emb->is_euc_in_hyp()) { ld diag = (u*u+v*v)/2; return matrix4( 1, 0, -u, u, 0, 1, -v, v, u, v, -diag+1, diag, u, v, -diag, diag+1 ); } 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 hyperpoint kleinize(hyperpoint h) { #if MAXMDIM == 3 return point3(h[0]/h[2], h[1]/h[2], 1); #else if(GDIM == 2) return point3(h[0]/h[2], h[1]/h[2], 1); else return point31(h[0]/h[3], h[1]/h[3], h[2]/h[3]); #endif } EX hyperpoint deparabolic13(hyperpoint h) { if(euclid) return h; if(cgi.emb->is_euc_in_hyp()) { h /= (1 + h[LDIM]); h[2] -= 1; h /= sqhypot_d(LDIM, h); h[2] += .5; return point3(h[0] * 2, h[1] * 2, log(2) + log(-h[2])); } h /= (1 + h[LDIM]); h[0] -= 1; h /= sqhypot_d(LDIM, h); h[0] += .5; return point3(log(2) + log(-h[0]), h[1] * 2, LDIM==3 ? h[2] * 2 : 0); } EX hyperpoint parabolic13(hyperpoint h) { if(euclid) return h; else if(cgi.emb->is_euc_in_hyp()) { return parabolic13(h[0], h[1]) * cpush0(2, h[2]); } else if(LDIM == 3) return parabolic13(h[1], h[2]) * xpush0(h[0]); else return parabolic1(h[1]) * xpush0(h[0]); } EX transmatrix parabolic13_at(hyperpoint h) { if(euclid) return rgpushxto0(h); else if(cgi.emb->is_euc_in_hyp()) { return parabolic13(h[0], h[1]) * cpush(2, h[2]); } else if(LDIM == 3) return parabolic13(h[1], h[2]) * xpush(h[0]); else return parabolic1(h[1]) * xpush(h[0]); } EX transmatrix spintoc(const hyperpoint& H, int t, int f) { transmatrix T = Id; ld R = hypot(H[f], H[t]); if(R >= 1e-15) { 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; } /** an Euclidean rotation in the axes (t,f) which rotates * the point H to the positive 't' axis */ EX transmatrix rspintoc(const hyperpoint& H, int t, int f) { transmatrix T = Id; ld R = hypot(H[f], H[t]); if(R >= 1e-15) { 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; } /** an isometry which takes the point H to the positive X axis * \see rspintox */ EX transmatrix spintox(const hyperpoint& H) { if(GDIM == 2 || gproduct) 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 || gproduct) 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 g(geometry, gSphere); fixmatrix(rot); fix4(rot); } /** determinant 2x2 */ EX ld det2(const transmatrix& T) { return T[0][0] * T[1][1] - T[0][1] * T[1][0]; } /** determinant 3x3 */ EX ld det3(const transmatrix& T) { ld det = 0; for(int i=0; i<3; i++) det += T[0][i] * T[1][(i+1)%3] * T[2][(i+2)%3]; for(int i=0; i<3; i++) det -= T[0][i] * T[1][(i+2)%3] * T[2][(i+1)%3]; return det; } /** determinant */ EX ld det(const transmatrix& T) { if(MDIM == 3) return det3(T); else { ld det = 1; transmatrix M = T; for(int a=0; a abs(M[max_at][a])) max_at = b; if(max_at != a) for(int c=a; c abs(T1[best][a])) best = b; int b = best; if(b != a) for(int c=0; c=0; a--) { for(int b=0; b product_decompose(hyperpoint h) { ld z = zlevel(h); return make_pair(z, scale_point(h, exp(-z))); } /** distance from mh and 0 */ EX ld hdist0(const hyperpoint& mh) { switch(cgclass) { case gcHyperbolic: if(mh[LDIM] < 1) return 0; return acosh(mh[LDIM]); case gcEuclid: { return hypot_d(GDIM, mh); } case gcSphere: { ld res = mh[LDIM] >= 1 ? 0 : mh[LDIM] <= -1 ? M_PI : acos(mh[LDIM]); return res; } case gcProduct: { auto d1 = product_decompose(mh); return hypot(PIU(hdist0(d1.second)), d1.first); } #if MAXMDIM >= 4 case gcSL2: { auto cosh_r = hypot(mh[2], mh[3]); auto phi = atan2(mh[2], mh[3]); return hypot(cosh_r < 1 ? 0 : acosh(cosh_r), phi); } case gcNil: { ld bz = mh[0] * mh[1] / 2; return hypot(mh[0], mh[1]) + abs(mh[2] - bz); } #endif default: return hypot_d(GDIM, mh); } } EX ld hdist0(const shiftpoint& mh) { return hdist0(unshift(mh)); } /** length of a circle of radius r */ EX ld circlelength(ld r) { switch(cgclass) { case gcEuclid: return TAU * r; case gcHyperbolic: return TAU * sinh(r); case gcSphere: return TAU * sin(r); default: return TAU * r; } } /* distance between h1 and h2 */ EX ld hdist(const hyperpoint& h1, const hyperpoint& 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); case gcProduct: { auto d1 = product_decompose(h1); auto d2 = product_decompose(h2); return hypot(PIU(hdist(d1.second, d2.second)), d1.first - d2.first); } case gcSL2: return hdist0(stretch::itranslate(h1) * h2); default: if(iv < 0) return 0; return sqrt(iv); } } EX ld hdist(const shiftpoint& h1, const shiftpoint& h2) { return hdist(h1.h, unshift(h2, h1.shift)); } /** like orthogonal_move but fol may be factor (in 2D graphics) or level (elsewhere) */ EX hyperpoint orthogonal_move_fol(const hyperpoint& h, double fol) { if(GDIM == 2) return scale_point(h, fol); else return orthogonal_move(h, fol); } /** like orthogonal_move but fol may be factor (in 2D graphics) or level (elsewhere) */ EX transmatrix orthogonal_move_fol(const transmatrix& T, double fol) { if(GDIM == 2) return scale_matrix(T, fol); else return orthogonal_move(T, fol); } /** like orthogonal_move but fol may be factor (in 2D graphics) or level (elsewhere) */ EX shiftmatrix orthogonal_move_fol(const shiftmatrix& T, double fol) { if(GDIM == 2) return scale_matrix(T, fol); else return orthogonal_move(T, fol); } /** the scaling matrix (Euclidean homogeneous scaling; also shift by log(scale) in product space */ EX transmatrix scale_matrix(const transmatrix& t, ld scale_factor) { transmatrix res; for(int i=0; i 3 for(int i=0; i<4; i++) t[3][i] = t[i][3] = i == 3; #endif return t; } EX transmatrix xyscale(const transmatrix& t, double fac) { transmatrix res; for(int i=0; i= 4 /** 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]; h2 /= h2[3]; hyperpoint w = h0; hyperpoint d1 = h1 - h0; hyperpoint d2 = h2 - h0; ld denom = (d1|d1) * (d2|d2) - (d1|d2) * (d1|d2); ld a1 = (d2|w) * (d1|d2) - (d1|w) * (d2|d2); ld a2 = (d1|w) * (d1|d2) - (d2|w) * (d1|d1); hyperpoint h = w * denom + d1 * a1 + d2 * a2; return normalize(h); } #endif EX hyperpoint hpxd(ld d, ld x, ld y, ld z) { hyperpoint H = hpxyz(d*x, d*y, z); H = mid(H, H); return H; } EX ld signum(ld x) { return x<0?-1:x>0?1:0; } EX bool asign(ld y1, ld y2) { return signum(y1) != signum(y2); } EX ld xcross(ld x1, ld y1, ld x2, ld y2) { return x1 + (x2 - x1) * y1 / (y1 - y2); } #if HDR enum eShiftMethod { smProduct, smIsotropic, smEmbedded, smLie, smGeodesic }; enum eEmbeddedShiftMethodChoice { smcNone, smcBoth, smcAuto }; enum eShiftMethodApplication { smaManualCamera, smaAutocenter, smaObject, smaWallRadar, smaAnimation }; #endif EX eEmbeddedShiftMethodChoice embedded_shift_method_choice = smcBoth; EX bool use_embedded_shift(eShiftMethodApplication sma) { switch(sma) { case smaAutocenter: case smaAnimation: return embedded_shift_method_choice != smcNone; case smaManualCamera: return embedded_shift_method_choice == smcBoth; case smaObject: return true; case smaWallRadar: return among(pmodel, mdLiePerspective, mdLieOrthogonal); default: throw hr_exception("unknown sma"); } } EX eShiftMethod shift_method(eShiftMethodApplication sma) { if(gproduct) return smProduct; if(embedded_plane && sma == smaObject) return smEmbedded; if(embedded_plane && use_embedded_shift(sma)) return smEmbedded; if(!nonisotropic && !stretch::in() && !(!nisot::geodesic_movement && hyperbolic && bt::in())) return smIsotropic; if(!nisot::geodesic_movement && !embedded_plane) return smLie; return smGeodesic; } EX transmatrix shift_object(transmatrix Position, const transmatrix& ori, const hyperpoint direction, eShiftMethod sm IS(shift_method(smaObject))) { switch(sm) { case smGeodesic: return nisot::parallel_transport(Position, direction); case smLie: return nisot::lie_transport(Position, direction); case smProduct: { hyperpoint h = product::direct_exp(ori * direction); return Position * rgpushxto0(h); } case smIsotropic: { return Position * rgpushxto0(direct_exp(direction)); } case smEmbedded: { // optimization -- should work without it if(cgi.emb->is_same_in_same()) return Position * rgpushxto0(direct_exp(direction)); if(gproduct) { return Position * cgi.emb->intermediate_to_actual_translation(ori * cgi.emb->logical_to_intermediate * direction); } auto P = inverse_shift(ggmatrix(cwt.at), shiftless(Position)); hyperpoint a = P * tile_center(); hyperpoint i0 = cgi.emb->actual_to_intermediate(a); auto l0 = cgi.emb->intermediate_to_logical * i0; auto l = l0; l[2] = 0; auto i = cgi.emb->logical_to_intermediate * l; auto rot1 = inverse(cgi.emb->intermediate_to_actual_translation(i)) * P; auto rot = cgi.emb->intermediate_to_logical_scaled * rot1 * cgi.emb->logical_scaled_to_intermediate; auto par = cgi.emb->logical_to_intermediate * rot * direction; if(cgi.emb->is_sph_in_low()) par = cspin180(0, 1) * par; // reason unknown auto tra = cgi.emb->intermediate_to_actual_translation(par); return Position * inverse(rot1) * tra * rot1; } default: throw hr_exception("unknown shift method in shift_object"); } } EX void apply_shift_object(transmatrix& Position, const transmatrix orientation, const hyperpoint direction, eShiftMethod sm IS(shift_method(smaObject))) { Position = shift_object(Position, orientation, direction, sm); } EX void rotate_object(transmatrix& Position, transmatrix& orientation, transmatrix R) { if(cgi.emb->is_euc_in_product()) orientation = orientation * R; else if(gproduct && WDIM == 3) orientation = orientation * R; else Position = Position * R; } EX transmatrix spin_towards(const transmatrix Position, transmatrix& ori, const hyperpoint goal, int dir, int back) { transmatrix T; ld alpha = 0; if(nonisotropic && nisot::geodesic_movement) T = nisot::spin_towards(Position, goal); else { hyperpoint U = inverse(Position) * goal; if(gproduct) { hyperpoint h = product::inverse_exp(U); alpha = asin_clamp(h[2] / hypot_d(3, h)); U = product_decompose(U).second; } T = rspintox(U); } if(back < 0) T = T * spin180(), alpha = -alpha; if(gproduct) { if(dir == 0) ori = cspin(2, 0, alpha); if(dir == 2) ori = cspin(2, 0, alpha - 90._deg), dir = 0; } if(dir) T = T * cspin(dir, 0, -90._deg); T = Position * T; return T; } EX shiftmatrix spin_towards(const shiftmatrix Position, transmatrix& ori, const shiftpoint goal, int dir, int back) { return shiftless(spin_towards(Position.T, ori, unshift(goal, Position.shift), dir, back), Position.shift); } 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; } EX transmatrix transpose(transmatrix T) { transmatrix result; for(int i=0; i= 4 if(nil) return nilv::formula_exp(v); if(sl2 || stretch::in()) return stretch::mstretch ? nisot::numerical_exp(v) : rots::formula_exp(v); #endif if(gproduct) return product::direct_exp(v); ld d = hypot_d(GDIM, v); if(d > 0) for(int i=0; i 90._deg) return M_PI - d; return d; } EX hyperpoint lp_iapply(const hyperpoint h) { return nisot::local_perspective_used ? inverse(NLP) * h : h; } EX hyperpoint lp_apply(const hyperpoint h) { return nisot::local_perspective_used ? NLP * h : h; } EX hyperpoint smalltangent() { return xtangent(.1); } EX void cyclefix(ld& a, ld b) { while(a > b + M_PI) a -= TAU; while(a < b - M_PI) a += TAU; } EX ld raddif(ld a, ld b) { ld d = a-b; if(d < 0) d = -d; if(d > TAU) d -= TAU; if(d > M_PI) d = TAU-d; return d; } EX int bucket_scale = 10000; EX unsigned bucketer(ld x) { return (unsigned) (long long) (floor(x * bucket_scale)); } EX unsigned bucketer(hyperpoint h) { unsigned dx = 0; if(gproduct) { auto d = product_decompose(h); h = d.second; dx += bucketer(d.first) * 50; if(cgi.emb->is_euc_in_product() && in_h2xe()) h /= h[2]; } dx += bucketer(h[0]) + 1000 * bucketer(h[1]) + 1000000 * bucketer(h[2]); if(MDIM == 4) dx += bucketer(h[3]) * 1000000001; if(elliptic) dx = min(dx, -dx); return dx; } #if MAXMDIM >= 4 /** @brief project the origin to the triangle [h1,h2,h3] */ EX hyperpoint project_on_triangle(hyperpoint h1, hyperpoint h2, hyperpoint h3) { h1 /= h1[3]; h2 /= h2[3]; h3 /= h3[3]; transmatrix T; T[0] = h1; T[1] = h2; T[2] = h3; T[3] = C0; ld det_orig = det3(T); hyperpoint orthogonal = (h2 - h1) ^ (h3 - h1); T[0] = orthogonal; T[1] = h2-h1; T[2] = h3-h1; ld det_orth = det3(T); hyperpoint result = orthogonal * (det_orig / det_orth); result[3] = 1; return normalize(result); } #endif EX hyperpoint lerp(hyperpoint a0, hyperpoint a1, ld x) { return a0 + (a1-a0) * x; } EX hyperpoint linecross(hyperpoint a, hyperpoint b, hyperpoint c, hyperpoint d) { a /= a[LDIM]; b /= b[LDIM]; c /= c[LDIM]; d /= d[LDIM]; ld bax = b[0] - a[0]; ld dcx = d[0] - c[0]; ld cax = c[0] - a[0]; ld bay = b[1] - a[1]; ld dcy = d[1] - c[1]; ld cay = c[1] - a[1]; hyperpoint res; res[0] = (cay * dcx * bax + a[0] * bay * dcx - c[0] * dcy * bax) / (bay * dcx - dcy * bax); res[1] = (cax * dcy * bay + a[1] * bax * dcy - c[1] * dcx * bay) / (bax * dcy - dcx * bay); res[2] = 0; res[3] = 0; res[GDIM] = 1; return normalize(res); } EX ld inner2(hyperpoint h1, hyperpoint h2) { return hyperbolic ? h1[LDIM] * h2[LDIM] - h1[0] * h2[0] - h1[1] * h2[1] : sphere ? h1[LDIM] * h2[LDIM] + h1[0] * h2[0] + h1[1] * h2[1] : h1[0] * h2[0] + h1[1] * h2[1]; } EX hyperpoint circumscribe(hyperpoint a, hyperpoint b, hyperpoint c) { hyperpoint h = C0; b = b - a; c = c - a; if(euclid) { ld b2 = inner2(b, b)/2; ld c2 = inner2(c, c)/2; ld det = c[1]*b[0] - b[1]*c[0]; h = a; h[1] += (c2*b[0] - b2 * c[0]) / det; h[0] += (c2*b[1] - b2 * c[1]) / -det; return h; } if(inner2(b,b) < 0) { b = b / sqrt(-inner2(b, b)); c = c + b * inner2(c, b); h = h + b * inner2(h, b); } else { b = b / sqrt(inner2(b, b)); c = c - b * inner2(c, b); h = h - b * inner2(h, b); } if(inner2(c,c) < 0) { c = c / sqrt(-inner2(c, c)); h = h + c * inner2(h, c); } else { c = c / sqrt(inner2(c, c)); h = h - c * inner2(h, c); } if(h[LDIM] < 0) h[0] = -h[0], h[1] = -h[1], h[LDIM] = -h[LDIM]; ld i = inner2(h, h); if(i > 0) h /= sqrt(i); else h /= -sqrt(-i); return h; } EX ld inner3(hyperpoint h1, hyperpoint h2) { return hyperbolic ? h1[LDIM] * h2[LDIM] - h1[0] * h2[0] - h1[1] * h2[1] - h1[2]*h2[2]: sphere ? h1[LDIM] * h2[LDIM] + h1[0] * h2[0] + h1[1] * h2[1] + h1[2]*h2[2]: h1[0] * h2[0] + h1[1] * h2[1]; } #if MAXMDIM >= 4 /** circumscribe for H3 and S3 (not for E3 yet!) */ EX hyperpoint circumscribe(hyperpoint a, hyperpoint b, hyperpoint c, hyperpoint d) { array ds = { b-a, c-a, d-a, C0 }; for(int i=0; i<3; i++) { if(inner3(ds[i],ds[i]) < 0) { ds[i] = ds[i] / sqrt(-inner3(ds[i], ds[i])); for(int j=i+1; j<4; j++) ds[j] = ds[j] + ds[i] * inner3(ds[i], ds[j]); } else { ds[i] = ds[i] / sqrt(inner3(ds[i], ds[i])); for(int j=i+1; j<4; j++) ds[j] = ds[j] - ds[i] * inner3(ds[i], ds[j]); } } hyperpoint& h = ds[3]; if(h[3] < 0) h = -h; ld i = inner3(h, h); if(i > 0) h /= sqrt(i); else h /= -sqrt(-i); return h; } #endif /** the point in distance dist from 'material' to 'dir' (usually an (ultra)ideal point) */ EX hyperpoint towards_inf(hyperpoint material, hyperpoint dir, ld dist IS(1)) { transmatrix T = gpushxto0(material); hyperpoint id = T * dir; return rgpushxto0(material) * rspintox(id) * xpush0(dist); } EX bool clockwise(hyperpoint h1, hyperpoint h2) { return h1[0] * h2[1] > h1[1] * h2[0]; } EX ld worst_precision_error; #if HDR struct hr_precision_error : hr_exception { hr_precision_error() : hr_exception("precision error") {} }; #endif /** check if a and b are the same, testing for equality. Throw an exception or warning if not sure */ EX bool same_point_may_warn(hyperpoint a, hyperpoint b) { ld d = hdist(a, b); if(d > 1e-2) return false; if(d > 1e-3) throw hr_precision_error(); if(d > 1e-6 && worst_precision_error <= 1e-6) addMessage("warning: precision errors are building up!"); if(d > worst_precision_error) worst_precision_error = d; return true; } }