nisot:: nicer precision setting

This commit is contained in:
Zeno Rogue 2020-04-11 15:08:24 +02:00
parent c5d7479d00
commit e4a8b73819
11 changed files with 81 additions and 67 deletions

View File

@ -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) {

View File

@ -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; i<cwt.at->type; 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);
}

View File

@ -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;
}

View File

@ -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;

View File

@ -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) {

View File

@ -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));
}

View File

@ -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<parallels; ai++) {
ld a = ai * 360 / parallels;
@ -1524,7 +1524,7 @@ namespace mapeditor {
displayfr(vid.xres-8, vid.yres-8-fs*5, 2, vid.fsize, XLAT("z: %1", fts(mh[2],4)), 0xC0C0C0, 16);
if(MDIM == 4)
displayfr(vid.xres-8, vid.yres-8-fs*4, 2, vid.fsize, XLAT("w: %1", fts(mh[3],4)), 0xC0C0C0, 16);
mh = inverse_exp(mh, iTable, false);
mh = inverse_exp(mh);
displayfr(vid.xres-8, vid.yres-8-fs*3, 2, vid.fsize, XLAT("r: %1", fts(hypot_d(3, mh),4)), 0xC0C0C0, 16);
if(GDIM == 3) {
displayfr(vid.xres-8, vid.yres-8-fs, 2, vid.fsize, XLAT("ϕ: %1°", fts(-atan2(mh[2], hypot_d(2, mh)) / degree,4)), 0xC0C0C0, 16);

View File

@ -552,7 +552,7 @@ EX namespace sn {
}
}
EX hyperpoint get_inverse_exp_symsol(hyperpoint h, bool lazy, bool just_direction) {
EX hyperpoint get_inverse_exp_symsol(hyperpoint h, flagtype flags) {
auto& s = get_tabled();
s.load();
@ -562,18 +562,17 @@ EX namespace sn {
if(h[2] < 0.) { iz = -iz; swap(ix, iy); }
hyperpoint res = s.get(ix, iy, iz, lazy);
hyperpoint res = s.get(ix, iy, iz, flags & pfNO_INTERPOLATION);
if(h[2] < 0.) { swap(res[0], res[1]); res[2] = -res[2]; }
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 hyperpoint get_inverse_exp_nsym(hyperpoint h, bool lazy, bool just_direction) {
EX hyperpoint get_inverse_exp_nsym(hyperpoint h, flagtype flags) {
auto& s = get_tabled();
s.load();
@ -581,14 +580,13 @@ EX namespace sn {
ld iy = h[1] >= 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<steps; i++) geodesic_step(at, v);
for(int i=0; i<rk_steps; i++) geodesic_step(at, v);
return at;
}
@ -2045,12 +2048,12 @@ EX namespace nisot {
return parallel_transport_bare(P, direction);
}
EX transmatrix spin_towards(const transmatrix Position, const hyperpoint goal) {
EX transmatrix spin_towards(const transmatrix Position, const hyperpoint goal, flagtype prec IS(pNORMAL)) {
hyperpoint at = tC0(Position);
transmatrix push_back = inverse(translate(at));
hyperpoint back_goal = push_back * goal;
back_goal = inverse_exp(back_goal, iTable);
back_goal = inverse_exp(back_goal, prec);
transmatrix back_Position = push_back * Position;
@ -2133,6 +2136,11 @@ EX namespace nisot {
shift_arg_formula(nilv::nilwidth);
return 0;
}
else if(argis("-rk-steps")) {
PHASEFROM(2);
shift(); rk_steps = argi();
return 0;
}
else if(argis("-nilv")) {
PHASEFROM(2);
if(nil) stop_game();

View File

@ -26,7 +26,7 @@ pair<bool, hyperpoint> 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};

View File

@ -122,8 +122,8 @@ EX namespace reg3 {
vertex_distance = binsearch(0, M_PI, [&] (ld d) {
// sometimes breaks in elliptic
dynamicval<eGeometry> 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; i<face; i++) mid += cspin(1, 2, 2*i*M_PI/face) * v2;

View File

@ -1104,7 +1104,7 @@ void movePlayer(monster *m, int delta) {
int i0 = i;
for(int a=0; a<3; a++) v[a] = (i0 % 3) - 1, i0 /= 3;
v = v * .1 / hypot_d(3, v);
transmatrix T1 = (i == 13) ? nat : parallel_transport(nat, m->ori, 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)