From e4a8b738198bd4f2c944518eec960c060c32829b Mon Sep 17 00:00:00 2001 From: Zeno Rogue Date: Sat, 11 Apr 2020 15:08:24 +0200 Subject: [PATCH] nisot:: nicer precision setting --- 3d-models.cpp | 8 ++++---- control.cpp | 6 +++--- drawing.cpp | 2 +- graph.cpp | 10 +++++----- hyperpoint.cpp | 28 +++++++++++++++++----------- hypgraph.cpp | 20 ++++++++++---------- mapeditor.cpp | 8 ++++---- nonisotropic.cpp | 46 +++++++++++++++++++++++++++------------------- radar.cpp | 2 +- reg3.cpp | 6 +++--- shmup.cpp | 12 ++++++------ 11 files changed, 81 insertions(+), 67 deletions(-) diff --git a/3d-models.cpp b/3d-models.cpp index 4f911df5..d57c5f3c 100644 --- a/3d-models.cpp +++ b/3d-models.cpp @@ -623,16 +623,16 @@ void geometry_information::slimetriangle(hyperpoint a, hyperpoint b, hyperpoint texture_order([&] (ld x, ld y) { ld z = 1-x-y; ld r = scalefactor * hcrossf7 * (0 + pow(max(x,max(y,z)), .3) * 0.8) * (hybri ? .5 : 1); - hyperpoint h = direct_exp(tangent_length(a*x+b*y+c*z, r), 10); + hyperpoint h = direct_exp(tangent_length(a*x+b*y+c*z, r)); hpcpush(h); }); } void geometry_information::balltriangle(hyperpoint a, hyperpoint b, hyperpoint c, ld rad, int lev) { if(lev == 0) { - hpcpush(direct_exp(a, 10)); - hpcpush(direct_exp(b, 10)); - hpcpush(direct_exp(c, 10)); + hpcpush(direct_exp(a)); + hpcpush(direct_exp(b)); + hpcpush(direct_exp(c)); } else { auto midpoint = [&] (hyperpoint h1, hyperpoint h2) { diff --git a/control.cpp b/control.cpp index add772e7..932cf127 100644 --- a/control.cpp +++ b/control.cpp @@ -76,7 +76,7 @@ EX movedir vectodir(hyperpoint P) { transmatrix U = ggmatrix(cwt.at); if(GDIM == 3 && WDIM == 2) U = radar_transform * U; - P = direct_exp(lp_iapply(P), 100); + P = direct_exp(lp_iapply(P)); hyperpoint H = sphereflip * tC0(U); transmatrix Centered = sphereflip * rgpushxto0(H); @@ -87,8 +87,8 @@ EX movedir vectodir(hyperpoint P) { for(int i=0; itype; i++) { transmatrix T = currentmap->adj(cwt.at, (cwt + i).spin); - ld d1 = geo_dist(U * T * C0, Centered * P, iTable); - ld d2 = geo_dist(U * T * C0, Centered * C0, iTable); + ld d1 = geo_dist(U * T * C0, Centered * P); + ld d2 = geo_dist(U * T * C0, Centered * C0); dirdist[i] = d1 - d2; //xspinpush0(-i * 2 * M_PI /cwt.at->type, .5), P); } diff --git a/drawing.cpp b/drawing.cpp index 441a5625..6ae2969f 100644 --- a/drawing.cpp +++ b/drawing.cpp @@ -409,7 +409,7 @@ void coords_to_poly() { bool behind3(hyperpoint h) { if(pmodel == mdGeodesic) - h = lp_apply(inverse_exp(h, iTable)); + h = lp_apply(inverse_exp(h)); return h[2] < 0; } diff --git a/graph.cpp b/graph.cpp index 257ecf60..1547a860 100644 --- a/graph.cpp +++ b/graph.cpp @@ -2884,7 +2884,7 @@ void apply_joukowsky_aura(hyperpoint& h) { h = ret; } if(nonisotropic) { - h = lp_apply(inverse_exp(h, iTable, true)); + h = lp_apply(inverse_exp(h, pfNO_DISTANCE)); } } @@ -3667,7 +3667,7 @@ bool celldrawer::cell_clipped() { hyperpoint H = tC0(V); if(abs(H[0]) <= 3 && abs(H[1]) <= 3 && abs(H[2]) <= 3 ) ; else { - hyperpoint H2 = inverse_exp(H, iLazy); + hyperpoint H2 = inverse_exp(H, pQUICK); for(hyperpoint& cpoint: clipping_planes) if((H2|cpoint) < -.4) return true; } noclipped++; @@ -3676,7 +3676,7 @@ bool celldrawer::cell_clipped() { hyperpoint H = tC0(V); if(abs(H[0]) <= 3 && abs(H[1]) <= 3 && abs(H[2]) <= 3 ) ; else { - hyperpoint H2 = inverse_exp(H, iLazy); + hyperpoint H2 = inverse_exp(H, pQUICK); for(hyperpoint& cpoint: clipping_planes) if((H2|cpoint) < -2) return true; } noclipped++; @@ -4470,12 +4470,12 @@ EX void precise_mouseover() { if(WDIM == 3) { mouseover2 = mouseover = centerover; ld best = HUGE_VAL; - hyperpoint h = direct_exp(lp_iapply(ztangent(0.01)), 100); + hyperpoint h = direct_exp(lp_iapply(ztangent(0.01))); transmatrix cov = ggmatrix(mouseover2); forCellIdEx(c1, i, mouseover2) { hyperpoint h1 = tC0(cov * currentmap->adj(mouseover2, i)); - ld dist = geo_dist(h, h1, iTable) - geo_dist(C0, h1, iTable); + ld dist = geo_dist(h, h1) - geo_dist(C0, h1); if(dist < best) mouseover = c1, best = dist; } return; diff --git a/hyperpoint.cpp b/hyperpoint.cpp index b4c98c41..a348cf04 100644 --- a/hyperpoint.cpp +++ b/hyperpoint.cpp @@ -1030,13 +1030,13 @@ 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); } -EX transmatrix parallel_transport(const transmatrix Position, const transmatrix& ori, const hyperpoint direction, int precision IS(100)) { +EX transmatrix parallel_transport(const transmatrix Position, const transmatrix& ori, const hyperpoint direction) { if(nonisotropic) return nisot::parallel_transport(Position, direction); else if(prod) { hyperpoint h = product::direct_exp(ori * direction); return Position * rgpushxto0(h); } - else return Position * rgpushxto0(direct_exp(direction, precision)); + else return Position * rgpushxto0(direct_exp(direction)); } EX void apply_parallel_transport(transmatrix& Position, const transmatrix orientation, const hyperpoint direction) { @@ -1151,8 +1151,8 @@ EX hyperpoint tangent_length(hyperpoint dir, ld length) { } /** exponential function: follow the geodesic given by v */ -EX hyperpoint direct_exp(hyperpoint v, int steps) { - if(sn::in()) return nisot::numerical_exp(v, steps); +EX hyperpoint direct_exp(hyperpoint v) { + if(sn::in()) return nisot::numerical_exp(v); if(nil) return nilv::formula_exp(v); if(sl2) return slr::formula_exp(v); if(prod) return product::direct_exp(v); @@ -1163,20 +1163,26 @@ EX hyperpoint direct_exp(hyperpoint v, int steps) { } #if HDR -enum iePrecision { iLazy, iTable }; +constexpr flagtype pfNO_INTERPOLATION = 1; /**< in tables (sol/nih geometries), do not use interpolations */ +constexpr flagtype pfNO_DISTANCE = 2; /**< we just need the directions -- this makes it a bit faster in sol/nih geometries */ +constexpr flagtype pfLOW_BS_ITER = 4; /**< low iterations in binary search (nil geometry, sl2 not affected currently) */ + +constexpr flagtype pQUICK = pfNO_INTERPOLATION | pfLOW_BS_ITER; + +constexpr flagtype pNORMAL = 0; #endif /** inverse exponential function \see hr::direct_exp */ -EX hyperpoint inverse_exp(const hyperpoint h, iePrecision p, bool just_direction IS(true)) { +EX hyperpoint inverse_exp(const hyperpoint h, flagtype prec IS(pNORMAL)) { #if CAP_SOLV if(sn::in()) { if(nih) - return sn::get_inverse_exp_nsym(h, p == iLazy, just_direction); + return sn::get_inverse_exp_nsym(h, prec); else - return sn::get_inverse_exp_symsol(h, p == iLazy, just_direction); + return sn::get_inverse_exp_symsol(h, prec); } #endif - if(nil) return nilv::get_inverse_exp(h, p == iLazy ? 5 : 20); + if(nil) return nilv::get_inverse_exp(h, prec); if(sl2) return slr::get_inverse_exp(h); if(prod) return product::inverse_exp(h); ld d = acos_auto_clamp(h[GDIM]); @@ -1186,9 +1192,9 @@ EX hyperpoint inverse_exp(const hyperpoint h, iePrecision p, bool just_direction return v; } -EX ld geo_dist(const hyperpoint h1, const hyperpoint h2, iePrecision p) { +EX ld geo_dist(const hyperpoint h1, const hyperpoint h2, flagtype prec IS(pNORMAL)) { if(!nonisotropic) return hdist(h1, h2); - return hypot_d(3, inverse_exp(inverse(nisot::translate(h1)) * h2, p, false)); + return hypot_d(3, inverse_exp(inverse(nisot::translate(h1)) * h2, prec)); } EX hyperpoint lp_iapply(const hyperpoint h) { diff --git a/hypgraph.cpp b/hypgraph.cpp index f060efea..90f41a92 100644 --- a/hypgraph.cpp +++ b/hypgraph.cpp @@ -337,7 +337,7 @@ EX void applymodel(hyperpoint H, hyperpoint& ret) { } case mdGeodesic: { - auto S = lp_apply(inverse_exp(H, iTable)); + auto S = lp_apply(inverse_exp(H, pNORMAL | pfNO_DISTANCE)); ld ratio = vid.xres / current_display->tanfov / current_display->radius / 2; ret[0] = S[0]/S[2] * ratio; ret[1] = S[1]/S[2] * ratio; @@ -360,7 +360,7 @@ EX void applymodel(hyperpoint H, hyperpoint& ret) { case mdDisk: { if(nonisotropic) { - ret = lp_apply(inverse_exp(H, iTable, true)); + ret = lp_apply(inverse_exp(H, pNORMAL | pfNO_DISTANCE)); ld w; if(sn::in()) { // w = 1 / sqrt(1 - sqhypot_d(3, ret)); @@ -546,7 +546,7 @@ EX void applymodel(hyperpoint H, hyperpoint& ret) { case mdFisheye: { ld zlev; if(nonisotropic) { - H = lp_apply(inverse_exp(H, iTable, false)); + H = lp_apply(inverse_exp(H)); zlev = 1; } else { @@ -768,7 +768,7 @@ EX void applymodel(hyperpoint H, hyperpoint& ret) { case mdEquidistant: case mdEquiarea: case mdEquivolume: { if(nonisotropic || prod) { - ret = lp_apply(inverse_exp(H, iTable, false)); + ret = lp_apply(inverse_exp(H)); ret[3] = 1; break; } @@ -1002,7 +1002,7 @@ EX transmatrix actualV(const heptspin& hs, const transmatrix& V) { EX bool point_behind(hyperpoint h) { if(sphere) return false; if(!in_perspective()) return false; - if(pmodel == mdGeodesic) h = inverse_exp(h, iLazy); + if(pmodel == mdGeodesic) h = inverse_exp(h, pQUICK); if(pmodel == mdPerspective && prod) h = product::inverse_exp(h); h = lp_apply(h); return h[2] < 1e-8; @@ -1995,7 +1995,7 @@ EX bool do_draw(cell *c, const transmatrix& T) { if(cells_drawn > vid.cells_drawn_limit) return false; if(cells_drawn < 50) { limited_generation(c); return true; } if(nil && pmodel == mdGeodesic) { - ld dist = hypot_d(3, inverse_exp(tC0(T), iLazy)); + ld dist = hypot_d(3, inverse_exp(tC0(T), pQUICK)); if(dist > sightranges[geometry] + (vid.sloppy_3d ? 0 : 0.9)) return false; if(dist <= extra_generation_distance && !limited_generation(c)) return false; } @@ -2004,7 +2004,7 @@ EX bool do_draw(cell *c, const transmatrix& T) { if(!limited_generation(c)) return false; } else if(pmodel == mdGeodesic && nih) { - hyperpoint h = inverse_exp(tC0(T), iLazy, false); + hyperpoint h = inverse_exp(tC0(T), pQUICK); ld dist = hypot_d(3, h); if(dist > sightranges[geometry] + (vid.sloppy_3d ? 0 : cgi.corner_bonus)) return false; if(dist <= extra_generation_distance && !limited_generation(c)) return false; @@ -2096,7 +2096,7 @@ EX void rotate_view(transmatrix T) { /** shift the view according to the given tangent vector */ EX transmatrix get_shift_view_of(const hyperpoint H, const transmatrix V) { if(!nonisotropic) { - return rgpushxto0(direct_exp(lp_iapply(H), 100)) * V; + return rgpushxto0(direct_exp(lp_iapply(H))) * V; } else if(!nisot::geodesic_movement) { transmatrix IV = inverse(V); @@ -2132,7 +2132,7 @@ void multiply_view(transmatrix T) { EX void shift_view_to(hyperpoint H) { if(!nonisotropic) multiply_view(gpushxto0(H)); - else shift_view(-inverse_exp(H, iTable, false)); + else shift_view(-inverse_exp(H)); } EX void shift_view_towards(hyperpoint H, ld l) { @@ -2141,7 +2141,7 @@ EX void shift_view_towards(hyperpoint H, ld l) { else if(nonisotropic && !nisot::geodesic_movement) shift_view(tangent_length(H-C0, -l)); else { - hyperpoint ie = inverse_exp(H, iTable, true); + hyperpoint ie = inverse_exp(H, pNORMAL | pfNO_DISTANCE); if(prod) ie = lp_apply(ie); shift_view(tangent_length(ie, -l)); } diff --git a/mapeditor.cpp b/mapeditor.cpp index 2258e102..d547daaf 100644 --- a/mapeditor.cpp +++ b/mapeditor.cpp @@ -1188,7 +1188,7 @@ namespace mapeditor { unsigned gridcolor = 0xC0C0C040; hyperpoint in_front_dist(ld d) { - return direct_exp(lp_iapply(ztangent(d)), 100); + return direct_exp(lp_iapply(ztangent(d))); } hyperpoint find_mouseh3() { @@ -1203,7 +1203,7 @@ namespace mapeditor { ld d1 = front_edit; hyperpoint h1 = in_front_dist(d); if(front_config == eFront::sphere_center) - d1 = geo_dist(drawtrans * C0, h1, iTable); + d1 = geo_dist(drawtrans * C0, h1); if(front_config == eFront::equidistants) { hyperpoint h = idt * in_front_dist(d); d1 = asin_auto(h[2]); @@ -1246,7 +1246,7 @@ namespace mapeditor { } if(front_config == eFront::sphere_center) for(int i=0; i<4; i+=2) { auto pt = [&] (ld a, ld b) { - return d2 * direct_exp(spin(a*degree) * cspin(0, 2, b*degree) * xtangent(front_edit), 100); + return d2 * direct_exp(spin(a*degree) * cspin(0, 2, b*degree) * xtangent(front_edit)); }; for(int ai=0; ai= 0. ? sn::x_to_ix(h[1]) : sn::x_to_ix(-h[1]); ld iz = sn::z_to_iz(h[2]); - hyperpoint res = s.get(ix, iy, iz, lazy); + hyperpoint res = s.get(ix, iy, iz, flags & pfNO_INTERPOLATION); if(h[0] < 0.) res[0] = -res[0]; if(h[1] < 0.) res[1] = -res[1]; - if(!just_direction) return table_to_azeq(res); - - return res; + if(flags & pfNO_DISTANCE) return res; + return table_to_azeq(res); } EX string shader_symsol = sn::common + @@ -745,7 +743,7 @@ EX namespace nilv { ); } - EX hyperpoint get_inverse_exp(hyperpoint h, int iterations) { + EX hyperpoint get_inverse_exp(hyperpoint h, flagtype prec IS(pNORMAL)) { ld wmin, wmax; ld side = h[2] - h[0] * h[1] / 2; @@ -769,11 +767,13 @@ EX namespace nilv { ld s = sin(2 * alpha_total); + int max_iter = (prec & pfLOW_BS_ITER) ? 5 : 20; + for(int it=0;; it++) { ld w = (wmin + wmax) / 2; ld z = b * b * (s + (sin(w) - w)/(cos(w) - 1)) + w; - if(it == iterations) { + if(it == max_iter) { ld alpha = alpha_total - w/2; ld c = b / sin(w/2); return point3(c * w * cos(alpha), c * w * sin(alpha), w); @@ -967,7 +967,7 @@ EX namespace nilv { EX hyperpoint on_geodesic(hyperpoint s0, hyperpoint s1, ld x) { hyperpoint local = inverse(nisot::translate(s0)) * s1; - hyperpoint h = get_inverse_exp(local, 100); + hyperpoint h = get_inverse_exp(local); return nisot::translate(s0) * formula_exp(h * x); } @@ -1822,7 +1822,7 @@ EX namespace rots { hybrid::in_underlying_geometry([&] { hyperpoint h = tC0(T); Spin = inverse(gpushxto0(h) * T); - d = hr::inverse_exp(h, iTable); + d = hr::inverse_exp(h); alpha = atan2(Spin[0][1], Spin[0][0]); distance = hdist0(h); beta = atan2(h[1], h[0]); @@ -1986,11 +1986,14 @@ EX namespace nisot { vel += (acc1+2*acc2+2*acc3+acc4)/6; } - EX hyperpoint numerical_exp(hyperpoint v, int steps) { + EX int rk_steps = 20; + EX int pt_steps = 100; + + EX hyperpoint numerical_exp(hyperpoint v) { hyperpoint at = point31(0, 0, 0); - v /= steps; + v /= rk_steps; v[3] = 0; - for(int i=0; i makeradar(hyperpoint h) { ld d = hdist0(h); if(sol && nisot::geodesic_movement) { - h = inverse_exp(h, iLazy); + h = inverse_exp(h, pQUICK); ld r = hypot_d(3, h); if(r < 1) h = h * (atanh(r) / r); else return {false, h}; diff --git a/reg3.cpp b/reg3.cpp index 2c689701..da4360fc 100644 --- a/reg3.cpp +++ b/reg3.cpp @@ -122,8 +122,8 @@ EX namespace reg3 { vertex_distance = binsearch(0, M_PI, [&] (ld d) { // sometimes breaks in elliptic dynamicval g(geometry, elliptic ? gCell120 : geometry); - hyperpoint v2 = direct_exp(dir_v2 * d, iTable); - hyperpoint v3 = direct_exp(dir_v3 * d, iTable); + hyperpoint v2 = direct_exp(dir_v2 * d); + hyperpoint v3 = direct_exp(dir_v3 * d); return hdist(v2, v3) >= edge_length; }); } @@ -131,7 +131,7 @@ EX namespace reg3 { DEBB(DF_GEOM, ("vertex_distance = ", vertex_distance)); /* actual vertex */ - hyperpoint v2 = direct_exp(dir_v2 * vertex_distance, iTable); + hyperpoint v2 = direct_exp(dir_v2 * vertex_distance); hyperpoint mid = Hypc; for(int i=0; iori, v, 2); + transmatrix T1 = (i == 13) ? nat : parallel_transport(nat, m->ori, v); cell *c3 = c2; while(true) { cell *c4 = findbaseAround(tC0(T1), c3, 1); @@ -1619,10 +1619,10 @@ void moveBullet(monster *m, int delta) { m->dead = true; if(inertia_based) { - nat = parallel_transport(nat, m->ori, m->inertia * delta, 10); + nat = parallel_transport(nat, m->ori, m->inertia * delta); } else - nat = parallel_transport(nat, m->ori, fronttangent(delta * SCALE * m->vel / speedfactor()), 10); + nat = parallel_transport(nat, m->ori, fronttangent(delta * SCALE * m->vel / speedfactor())); cell *c2 = m->findbase(nat, 1); if(m->parent && isPlayer(m->parent) && markOrb(itOrbLava) && c2 != m->base && !isPlayerOn(m->base)) @@ -2106,14 +2106,14 @@ void moveMonster(monster *m, int delta) { if(inertia_based) { if(igo) return; - nat = parallel_transport(nat, m->ori, m->inertia * delta, 10); + nat = parallel_transport(nat, m->ori, m->inertia * delta); } else if(WDIM == 3 && igo) { ld fspin = rand() % 1000; - nat = parallel_transport(nat0, m->ori, cspin(1,2,fspin) * spin(igospan[igo]) * xtangent(step), 10); + nat = parallel_transport(nat0, m->ori, cspin(1,2,fspin) * spin(igospan[igo]) * xtangent(step)); } else { - nat = parallel_transport(nat0, m->ori, spin(igospan[igo]) * xtangent(step), 10); + nat = parallel_transport(nat0, m->ori, spin(igospan[igo]) * xtangent(step)); } if(m->type != moRagingBull && !peace::on)