diff --git a/basegraph.cpp b/basegraph.cpp index 8b0d117f..220ce91f 100644 --- a/basegraph.cpp +++ b/basegraph.cpp @@ -314,6 +314,8 @@ void display_data::set_projection(int ed) { shaderside_projection = true, glhr::new_shader_projection = glhr::shader_projection::standardEH2, pers3 = true; if(GDIM == 3 && apply_models && pmodel == mdGeodesic && sol) shaderside_projection = true, glhr::new_shader_projection = glhr::shader_projection::standardSolv, pers3 = true; + if(GDIM == 3 && apply_models && pmodel == mdGeodesic && nih) + shaderside_projection = true, glhr::new_shader_projection = glhr::shader_projection::standardNIH, pers3 = true; if(GDIM == 3 && apply_models && pmodel == mdGeodesic && sl2) shaderside_projection = true, glhr::new_shader_projection = glhr::shader_projection::standardSL2, pers3 = true; if(GDIM == 3 && apply_models && pmodel == mdGeodesic && nil) @@ -333,56 +335,17 @@ void display_data::set_projection(int ed) { if(pmodel == mdRug) return; #if CAP_SOLV - if(glhr::new_shader_projection == glhr::shader_projection::standardSolv) { + if(among(glhr::new_shader_projection, glhr::shader_projection::standardSolv, glhr::shader_projection::standardNIH)) { + auto &tab = ((glhr::new_shader_projection == glhr::shader_projection::standardSolv) ? solv::solt : nihv::niht); - static bool toload = true; - - static GLuint invexpid = 0; - - if(toload) { - solv::load_table(); - if(!solv::table_loaded) { pmodel = mdPerspective; set_projection(ed); return; } - - println(hlog, "installing table"); - using namespace solv; - toload = false; - - if(invexpid == 0) glGenTextures(1, &invexpid); - - glBindTexture( GL_TEXTURE_3D, invexpid); - - glTexParameteri(GL_TEXTURE_3D, GL_TEXTURE_MIN_FILTER, GL_LINEAR); - glTexParameteri(GL_TEXTURE_3D, GL_TEXTURE_MAG_FILTER, GL_LINEAR); - - glTexParameteri(GL_TEXTURE_3D, GL_TEXTURE_WRAP_S, GL_CLAMP_TO_EDGE); - glTexParameteri(GL_TEXTURE_3D, GL_TEXTURE_WRAP_T, GL_CLAMP_TO_EDGE); - glTexParameteri(GL_TEXTURE_3D, GL_TEXTURE_WRAP_R, GL_CLAMP_TO_EDGE); - - auto xbuffer = new glvertex[PRECZ*PRECY*PRECX]; - - ld maxd = 0; - for(int z=0; z ginf = { {"product","none", "product space", "product", 7, 3, qHYBRID, giProduct, 0x00000, {{7, 3}}, eVariation::pure}, {"twisted","none", "rotation space", "twisted", 7, 3, qHYBRID, giSL2, 0x00000, {{6, 4}}, eVariation::pure}, {"ternary","none", "standard ternary tiling", "ternary", 6, 3, qBINARY, giHyperb2, 0x48400, {{6, 4}}, eVariation::pure}, + {"NIH", "none", "non-isotropic hyperbolic", "NIH", 11, 3, qBINARY | qEXPERIMENTAL, giH23, 0x49000, {{6, 4}}, eVariation::pure}, }; // bits: 9, 10, 15, 16, (reserved for later) 17, 18 diff --git a/classes.h b/classes.h index 51bc2a20..065705fd 100644 --- a/classes.h +++ b/classes.h @@ -215,10 +215,10 @@ enum eGeometry { gField435, gField534, gBinary4, gSol, gKiteDart2, gKiteDart3, gNil, gProduct, gRotSpace, - gTernary, + gTernary, gNIH, gGUARD}; -enum eGeometryClass { gcHyperbolic, gcEuclid, gcSphere, gcSol, gcNil, gcProduct, gcSL2 }; +enum eGeometryClass { gcHyperbolic, gcEuclid, gcSphere, gcSol, gcNil, gcProduct, gcSL2, gcNIH }; enum class eVariation { bitruncated, pure, goldberg, irregular, dual }; diff --git a/devmods/solv-table.cpp b/devmods/solv-table.cpp index 6c238740..09f6b8d0 100644 --- a/devmods/solv-table.cpp +++ b/devmods/solv-table.cpp @@ -14,12 +14,15 @@ const int _PRECZ = 64; transmatrix parabolic1(ld u); -namespace solv { +namespace nisot { typedef hyperpoint pt; -typedef array ptlow; -ptlow be_low(pt x) { return ptlow({float(x[0]), float(x[1]), float(x[2])}); } +using solnihv::x_to_ix; + +ld z_to_iz(ld z) { if(sol) return tanh(z); else return tanh(z/4)/2 + .5; } + +ptlow be_low(hyperpoint x) { return ptlow({float(x[0]), float(x[1]), float(x[2])}); } template void parallelize(int threads, int Nmin, int Nmax, T action) { std::vector v; @@ -30,35 +33,13 @@ template void parallelize(int threads, int Nmin, int Nmax, T action) { for(std::thread& t:v) t.join(); } -hyperpoint sol1(pt v) { - auto [x,y,z,t] = (array&) v; - if(x == 0 && z == 0) return C0; - hyperpoint h = parabolic1(x) * xpush(-z) * C0; - ld d = acosh(h[2]) / sqrt(h[0] * h[0] + h[1] * h[1]); - return hyperpoint({h[1]*d, 0, -h[0]*d,1}); - } - -hyperpoint sol2(pt v) { - auto [x,y,z,t] = (array&) v; - if(y == 0 && z == 0) return C0; - hyperpoint h = parabolic1(y) * xpush(z) * C0; - ld d = acosh(h[2]) / sqrt(h[0] * h[0] + h[1] * h[1]); - return hyperpoint({0, h[1]*d, +h[0]*d,1}); - } - -ld x_to_ix(ld u); - -ld solerror(pt ok, pt chk) { - auto zok = point3( x_to_ix(ok[0]), x_to_ix(ok[1]), tanh(ok[2]) ); - auto zchk = point3( x_to_ix(chk[0]), x_to_ix(chk[1]), tanh(chk[2]) ); +ld solerror(hyperpoint ok, hyperpoint chk) { + auto zok = point3( x_to_ix(ok[0]), x_to_ix(ok[1]), z_to_iz(ok[2]) ); + auto zchk = point3( x_to_ix(chk[0]), x_to_ix(chk[1]), z_to_iz(chk[2]) ); return hypot_d(3, zok - zchk); } -ld eucerror(pt ok, pt chk) { - return pow(ok[0]-chk[0], 2) + pow(ok[1]-chk[1], 2) + pow(ok[2]-chk[2], 2); - } - -pt iterative_solve(pt xp, pt candidate, int prec, ld minerr, bool debug = false) { +hyperpoint iterative_solve(hyperpoint xp, hyperpoint candidate, int prec, ld minerr, bool debug = false) { transmatrix T = Id; T[0][1] = 8; T[2][2] = 5; @@ -70,7 +51,7 @@ pt iterative_solve(pt xp, pt candidate, int prec, ld minerr, bool debug = false) ld eps = 1e-6; - pt c[3]; + hyperpoint c[3]; for(int a=0; a<3; a++) c[a] = point3(a==0, a==1, a==2); while(err > minerr) { @@ -145,13 +126,13 @@ ptlow solution[_PRECZ][_PRECY][_PRECX]; ptlow mlow(ld x, ld y, ld z) { return ptlow({float(x), float(y), float(z)}); } -pt atxyz(ld x, ld y, ld z) { return hyperpoint({x, y, z, 1}); } +hyperpoint atxyz(ld x, ld y, ld z) { return hyperpoint({x, y, z, 1}); } ptlow operator +(ptlow a, ptlow b) { return mlow(a[0]+b[0], a[1]+b[1], a[2]+b[2]); } ptlow operator -(ptlow a, ptlow b) { return mlow(a[0]-b[0], a[1]-b[1], a[2]-b[2]); } ptlow operator *(ptlow a, ld x) { return mlow(a[0]*x, a[1]*x, a[2]*x); } -ptlow can(pt x) { +ptlow can(hyperpoint x) { // azimuthal equidistant to Klein ld r = sqrt(x[0] * x[0] + x[1] * x[1] + x[2] * x[2]); if(r == 0) return mlow(0,0,0); @@ -160,7 +141,7 @@ ptlow can(pt x) { return mlow(x[0]*d, x[1]*d, x[2]*d); } -pt uncan(ptlow x) { +hyperpoint uncan(ptlow x) { ld r = sqrt(x[0] * x[0] + x[1] * x[1] + x[2] * x[2]); if(r == 0) return atxyz(0,0,0); ld make_r = atanh(r); @@ -169,7 +150,7 @@ pt uncan(ptlow x) { return atxyz(x[0]*d, x[1]*d, x[2]*d); } -pt uncan_info(ptlow x) { +hyperpoint uncan_info(ptlow x) { ld r = sqrt(x[0] * x[0] + x[1] * x[1] + x[2] * x[2]); println(hlog, "r = ", r); if(r == 0) return atxyz(0,0,0); @@ -213,11 +194,7 @@ ld ix_to_x(ld ix) { } ld iz_to_z(ld z) { - return atanh(z); // atanh(z * 2 - 1); - } - -ld z_to_iz(ld z) { - return tanh(z); // (tanh(z) + 1) / 2; + return nih ? atanh(z * 2 - 1)*4 : atanh(z); // atanh(z * 2 - 1); } int last_x = _PRECX-1, last_y = _PRECY-1, last_z = _PRECZ-1; @@ -232,6 +209,7 @@ void build_sols() { std::mutex file_mutex; ld max_err = 0; auto act = [&] (int tid, int iz) { + if((nih && iz == 0) || iz == _PRECZ-1) return; auto solve_at = [&] (int ix, int iy) { ld x = ix_to_x(ix / (_PRECX-1.)); @@ -240,18 +218,18 @@ void build_sols() { auto v = hyperpoint ({x,y,z,1}); - vector candidates; - pt cand; + vector candidates; + hyperpoint cand; candidates.push_back(atxyz(0,0,0)); static constexpr int prec = 100; - // sort(candidates.begin(), candidates.end(), [&] (pt a, pt b) { return solerror(v, direct_exp(a, prec)) > solerror(v, direct_exp(b, prec)); }); + // sort(candidates.begin(), candidates.end(), [&] (hyperpoint a, hyperpoint b) { return solerror(v, direct_exp(a, prec)) > solerror(v, direct_exp(b, prec)); }); // cand_best = candidates.back(); - vector solved_candidates; + vector solved_candidates; for(auto c: candidates) { auto solt = iterative_solve(v, c, prec, 1e-6); @@ -259,7 +237,7 @@ void build_sols() { if(solerror(v, nisot::numerical_exp(solt, prec)) < 1e-9) break; } - sort(solved_candidates.begin(), solved_candidates.end(), [&] (pt a, pt b) { return solerror(v, nisot::numerical_exp(a, prec)) > solerror(v, nisot::numerical_exp(b, prec)); }); + sort(solved_candidates.begin(), solved_candidates.end(), [&] (hyperpoint a, hyperpoint b) { return solerror(v, nisot::numerical_exp(a, prec)) > solerror(v, nisot::numerical_exp(b, prec)); }); cand = solved_candidates.back(); @@ -299,14 +277,16 @@ void build_sols() { } }; - parallelize(last_z, 0, last_z, act); + parallelize(_PRECZ, 0, _PRECZ, act); for(int x=0; x= 4 if(binarytiling) binary::build_tmatrix(); diff --git a/graph.cpp b/graph.cpp index 17e8605b..97e04d52 100644 --- a/graph.cpp +++ b/graph.cpp @@ -4462,6 +4462,7 @@ int get_darkval(cell *c, int d) { const int darkval_sol[8] = {0,2,4,4,0,2,4,4}; const int darkval_penrose[12] = {0, 2, 0, 2, 4, 4, 6, 6, 6, 6, 6, 6}; const int darkval_nil[8] = {6,6,0,3,6,6,0,3}; + const int darkval_nih[11] = {0,2,0,2,4,6,6,6,6,6,6}; if(sphere) return darkval_s12[d]; if(euclid && S7 == 6) return darkval_e6[d]; if(euclid && S7 == 12) return darkval_e12[d]; @@ -4470,6 +4471,7 @@ int get_darkval(cell *c, int d) { if(geometry == gHoroRec) return darkval_hrec[d]; if(penrose) return darkval_penrose[d]; if(sol) return darkval_sol[d]; + if(nih) return darkval_nih[d]; if(binarytiling) return darkval_hbt[d]; if(hyperbolic && S7 == 6) return darkval_e6[d]; if(hyperbolic && S7 == 12) return darkval_s12[d]; diff --git a/hyper.h b/hyper.h index e217a96f..0062b766 100644 --- a/hyper.h +++ b/hyper.h @@ -130,12 +130,14 @@ void addMessage(string s, char spamtype = 0); #define euclid (cgclass == gcEuclid) #define sphere (cgclass == gcSphere) #define sol (cgclass == gcSol) +#define nih (cgclass == gcNIH) +#define solnih (sol || nih) #define nil (cgclass == gcNil) #define sl2 (cgclass == gcSL2) #define prod (cgclass == gcProduct) #define hybri (ginf[geometry].flags & qHYBRID) #define hyperbolic (cgclass == gcHyperbolic) -#define nonisotropic (sol || nil || sl2) +#define nonisotropic (sol || nil || sl2 || nih) #define translatable (euclid || nonisotropic) #define nonorientable (ginf[geometry].flags & qNONORIENTABLE) #define elliptic (ginf[geometry].flags & qELLIPTIC) diff --git a/hyperpoint.cpp b/hyperpoint.cpp index 9532045a..61907a47 100644 --- a/hyperpoint.cpp +++ b/hyperpoint.cpp @@ -1114,7 +1114,7 @@ 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(sol) return nisot::numerical_exp(v, steps); + if(solnih) return nisot::numerical_exp(v, steps); if(nil) return nilv::formula_exp(v); if(sl2) return slr::formula_exp(v); if(prod) return product::direct_exp(v); @@ -1132,6 +1132,7 @@ enum iePrecision { iLazy, iTable }; EX hyperpoint inverse_exp(const hyperpoint h, iePrecision p, bool just_direction IS(true)) { #if CAP_SOLV if(sol) return solv::get_inverse_exp(h, p == iLazy, just_direction); + if(nih) return nihv::get_inverse_exp(h, p == iLazy, just_direction); #endif if(nil) return nilv::get_inverse_exp(h, p == iLazy ? 5 : 20); if(sl2) return slr::get_inverse_exp(h); diff --git a/hypgraph.cpp b/hypgraph.cpp index 70bf0810..7dee9bb3 100644 --- a/hypgraph.cpp +++ b/hypgraph.cpp @@ -967,7 +967,7 @@ EX bool invalid_point(const hyperpoint h) { EX bool in_smart_range(const transmatrix& T) { hyperpoint h = tC0(T); if(invalid_point(h)) return false; - if(nil) return true; + if(nil || nih) return true; #if CAP_SOLV if(pmodel == mdGeodesic) return solv::in_table_range(h); #endif @@ -2011,6 +2011,13 @@ EX bool do_draw(cell *c, const transmatrix& T) { if(!nisot::in_table_range(tC0(T))) return false; if(!limited_generation(c)) return false; } + else if(pmodel == mdGeodesic && nih) { + hyperpoint h = inverse_exp(tC0(T), iLazy, false); + ld dist = hypot_d(3, h); + if(isnan(dist)) binary::breakhere(); + if(dist > sightranges[geometry] + (vid.sloppy_3d ? 0 : cgi.corner_bonus)) return false; + if(dist <= extra_generation_distance && !limited_generation(c)) return false; + } else if(pmodel == mdGeodesic && sl2) { if(hypot(tC0(T)[2], tC0(T)[3]) > cosh(slr::range_xy)) return false; if(!limited_generation(c)) return false; diff --git a/landgen.cpp b/landgen.cpp index 0980726c..29e8dcd6 100644 --- a/landgen.cpp +++ b/landgen.cpp @@ -2595,7 +2595,7 @@ EX void setdist(cell *c, int d, cell *from) { if(d >= BARLEV) { - if(binarytiling && WDIM == 3 && !c->land && !sol) { + if(binarytiling && WDIM == 3 && !c->land && !solnih) { ld z = vid.binary_width; cell *cseek = c; int step = 0; diff --git a/nonisotropic.cpp b/nonisotropic.cpp index 5239a69d..52a63edf 100644 --- a/nonisotropic.cpp +++ b/nonisotropic.cpp @@ -29,6 +29,10 @@ EX namespace nisot { T[0][0] = exp(-h[2]); T[1][1] = exp(+h[2]); } + if(nih) { + T[0][0] = pow(2, h[2]); + T[1][1] = pow(3, h[2]); + } if(nil) T[2][1] = h[0]; return T; @@ -36,37 +40,110 @@ EX namespace nisot { EX } -EX namespace solv { - #if CAP_SOLV - EX int PRECX, PRECY, PRECZ; +#if CAP_SOLV +EX namespace solnihv { + + #if HDR + struct tabled_inverses { + int PRECX, PRECY, PRECZ; + vector tab; + string fname; + bool loaded; + + void load(); + hyperpoint get(ld ix, ld iy, ld iz, bool lazy); - EX vector inverse_exp_table; + GLuint texture_id; + bool toload; + + GLuint get_texture_id(); - EX bool table_loaded; + tabled_inverses(string s) : fname(s), texture_id(0), toload(true) {} + }; + #endif - EX string solfname = "solv-geodesics.dat"; - - EX void load_table() { - if(table_loaded) return; - FILE *f = fopen(solfname.c_str(), "rb"); + void tabled_inverses::load() { + if(loaded) return; + FILE *f = fopen(fname.c_str(), "rb"); // if(!f) f = fopen("/usr/lib/soltable.dat", "rb"); if(!f) { addMessage(XLAT("geodesic table missing")); pmodel = mdPerspective; return; } ignore(fread(&PRECX, 4, 1, f)); ignore(fread(&PRECY, 4, 1, f)); ignore(fread(&PRECZ, 4, 1, f)); - inverse_exp_table.resize(PRECX * PRECY * PRECZ); - ignore(fread(&inverse_exp_table[0], sizeof(nisot::ptlow) * PRECX * PRECY * PRECZ, 1, f)); + tab.resize(PRECX * PRECY * PRECZ); + ignore(fread(&tab[0], sizeof(nisot::ptlow) * PRECX * PRECY * PRECZ, 1, f)); fclose(f); - table_loaded = true; + loaded = true; } - hyperpoint christoffel(const hyperpoint at, const hyperpoint velocity, const hyperpoint transported) { - return hpxyz3( - -velocity[2] * transported[0] - velocity[0] * transported[2], - velocity[2] * transported[1] + velocity[1] * transported[2], - velocity[0] * transported[0] * exp(2*at[2]) - velocity[1] * transported[1] * exp(-2*at[2]), - 0 - ); + hyperpoint tabled_inverses::get(ld ix, ld iy, ld iz, bool lazy) { + ix *= PRECX-1; + iy *= PRECY-1; + iz *= PRECZ-1; + + hyperpoint res; + + if(lazy) { + auto r = tab[(int(iz)*PRECY+int(iy))*PRECX+int(ix)]; + for(int i=0; i<3; i++) res[i] = r[i]; + } + + else { + + if(ix >= PRECX-1) ix = PRECX-2; + if(iy >= PRECX-1) iy = PRECX-2; + if(iz >= PRECZ-1) iz = PRECZ-2; + + int ax = ix, bx = ax+1; + int ay = iy, by = ay+1; + int az = iz, bz = az+1; + + #define S0(x,y,z) tab[(z*PRECY+y)*PRECX+x][t] + #define S1(x,y) (S0(x,y,az) * (bz-iz) + S0(x,y,bz) * (iz-az)) + #define S2(x) (S1(x,ay) * (by-iy) + S1(x,by) * (iy-ay)) + + for(int t=0; t<3; t++) + res[t] = S2(ax) * (bx-ix) + S2(bx) * (ix-ax); + } + + return res; + } + + GLuint tabled_inverses::get_texture_id() { + if(!toload) return texture_id; + + load(); + if(!loaded) return 0; + + println(hlog, "installing table"); + toload = false; + + if(texture_id == 0) glGenTextures(1, &texture_id); + + glBindTexture( GL_TEXTURE_3D, texture_id); + + glTexParameteri(GL_TEXTURE_3D, GL_TEXTURE_MIN_FILTER, GL_LINEAR); + glTexParameteri(GL_TEXTURE_3D, GL_TEXTURE_MAG_FILTER, GL_LINEAR); + + glTexParameteri(GL_TEXTURE_3D, GL_TEXTURE_WRAP_S, GL_CLAMP_TO_EDGE); + glTexParameteri(GL_TEXTURE_3D, GL_TEXTURE_WRAP_T, GL_CLAMP_TO_EDGE); + glTexParameteri(GL_TEXTURE_3D, GL_TEXTURE_WRAP_R, GL_CLAMP_TO_EDGE); + + auto xbuffer = new glvertex[PRECZ*PRECY*PRECX]; + + for(int z=0; z= 0. ? x_to_ix(h[0]) : x_to_ix(-h[0]); - ld iy = h[1] >= 0. ? x_to_ix(h[1]) : x_to_ix(-h[1]); - ld iz = tanh(h[2]); - - if(h[2] < 0.) { iz = -iz; swap(ix, iy); } - - ix *= PRECX-1; - iy *= PRECY-1; - iz *= PRECZ-1; - - hyperpoint res = C0; - - if(lazy) { - auto r = inverse_exp_table[(int(iz)*PRECY+int(iy))*PRECX+int(ix)]; - for(int i=0; i<3; i++) res[i] = r[i]; - } - - else { - - if(ix >= PRECX-1) ix = PRECX-2; - if(iy >= PRECX-1) iy = PRECX-2; - if(iz >= PRECZ-1) iz = PRECZ-2; - - int ax = ix, bx = ax+1; - int ay = iy, by = ay+1; - int az = iz, bz = az+1; - - #define S0(x,y,z) inverse_exp_table[(z*PRECY+y)*PRECX+x][t] - #define S1(x,y) (S0(x,y,az) * (bz-iz) + S0(x,y,bz) * (iz-az)) - #define S2(x) (S1(x,ay) * (by-iy) + S1(x,by) * (iy-ay)) - - for(int t=0; t<3; t++) - res[t] = S2(ax) * (bx-ix) + S2(bx) * (ix-ax); - } - - 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) { - ld r = hypot_d(3, res); - if(r == 0.) return res; - return res * atanh(r) / r; - } - - return res; - } - - struct hrmap_sol : hrmap { + struct hrmap_solnih : hrmap { hrmap *binary_map; + hrmap *ternary_map; /* nih only */ unordered_map, heptagon*> at; unordered_map> coords; @@ -159,9 +186,10 @@ EX namespace solv { return h; } - hrmap_sol() { + hrmap_solnih() { heptagon *alt; + heptagon *alt3; if(true) { dynamicval g(geometry, gBinary4); @@ -176,7 +204,24 @@ EX namespace solv { binary_map = binary::new_alt_map(alt); } - origin = get_at(alt, alt); + if(nih) { + dynamicval g(geometry, gTernary); + alt3 = tailored_alloc (S7); + alt3->s = hsOrigin; + alt3->alt = alt3; + alt3->cdata = NULL; + alt3->c7 = NULL; + alt3->zebraval = 0; + alt3->distance = 0; + alt3->emeraldval = 0; + ternary_map = binary::new_alt_map(alt3); + } + else { + alt3 = alt; + ternary_map = nullptr; + } + + origin = get_at(alt, alt3); } heptagon *altstep(heptagon *h, int d) { @@ -185,6 +230,12 @@ EX namespace solv { return h->cmove(d); } + heptagon *altstep3(heptagon *h, int d) { + dynamicval g(geometry, gTernary); + dynamicval cm(currentmap, ternary_map); + return h->cmove(d); + } + heptagon *create_step(heptagon *parent, int d) override { auto p = coords[parent]; auto pf = p.first, ps = p.second; @@ -194,7 +245,7 @@ EX namespace solv { return g; }; - switch(d) { + if(sol) switch(d) { case 0: // right return rule(altstep(pf, 2), ps, 4); case 1: // up @@ -214,27 +265,57 @@ EX namespace solv { default: return NULL; } + + else switch(d) { + case 0: // right + return rule(altstep(pf, 2), ps, 2); + case 1: // up + return rule(pf, altstep3(ps, 3), 3); + case 2: // left + return rule(altstep(pf, 4), ps, 0); + case 3: // down + return rule(pf, altstep3(ps, 5), 1); + case 4: // back + return rule(altstep(pf, 3), altstep3(ps, 4), 5 + pf->zebraval + 2 * ps->zebraval); + default: + return rule(altstep(pf, (d-5) % 2), altstep3(ps, (d-5)/2), 4); + } } - ~hrmap_sol() { + ~hrmap_solnih() { delete binary_map; + if(ternary_map) delete ternary_map; for(auto& p: at) clear_heptagon(p.second); } - + transmatrix adjmatrix(int i, int j) { - ld z = log(2); - ld bw = vid.binary_width * z; - ld bwh = bw / 4; - switch(i) { - case 0: return xpush(+bw); - case 1: return ypush(+bw); - case 2: return xpush(-bwh) * zpush(+z) * ypush(j == 6 ? +bwh : -bwh); - case 3: return xpush(+bwh) * zpush(+z) * ypush(j == 6 ? +bwh : -bwh); - case 4: return xpush(-bw); - case 5: return ypush(-bw); - case 6: return ypush(-bwh) * zpush(-z) * xpush(j == 2 ? +bwh : -bwh); - case 7: return ypush(+bwh) * zpush(-z) * xpush(j == 2 ? +bwh : -bwh); - default:return Id; + if(sol) { + ld z = log(2); + ld bw = vid.binary_width * z; + ld bwh = bw / 4; + switch(i) { + case 0: return xpush(+bw); + case 1: return ypush(+bw); + case 2: return xpush(-bwh) * zpush(+z) * ypush(j == 6 ? +bwh : -bwh); + case 3: return xpush(+bwh) * zpush(+z) * ypush(j == 6 ? +bwh : -bwh); + case 4: return xpush(-bw); + case 5: return ypush(-bw); + case 6: return ypush(-bwh) * zpush(-z) * xpush(j == 2 ? +bwh : -bwh); + case 7: return ypush(+bwh) * zpush(-z) * xpush(j == 2 ? +bwh : -bwh); + default:return Id; + } + } + else { + ld bw = vid.binary_width; + switch(i) { + case 0: return xpush(+bw); + case 1: return ypush(+bw); + case 2: return xpush(-bw); + case 3: return ypush(-bw); + case 4: return xpush(-((j-5)%2-.5)*bw) * ypush(-((j-5)/2-1)*bw) * zpush(1); + default: + return zpush(-1) * xpush(((i-5)%2-.5)*bw) * ypush(((i-5)/2-1)*bw); + } } } @@ -271,32 +352,16 @@ EX namespace solv { }; EX pair getcoord(heptagon *h) { - return ((hrmap_sol*)currentmap)->coords[h]; + return ((hrmap_solnih*)currentmap)->coords[h]; } EX heptagon *get_at(heptagon *h1, heptagon *h2, bool gen) { - auto m = ((hrmap_sol*)currentmap); + auto m = ((hrmap_solnih*)currentmap); if(!gen && !m->at.count(make_pair(h1, h2))) return nullptr; return m->get_at(h1, h2); } - EX ld solrange_xy = 15; - EX ld solrange_z = 4; - - EX bool in_table_range(hyperpoint h) { - return abs(h[0]) < solrange_xy && abs(h[1]) < solrange_xy && abs(h[2]) < solrange_z; - } - - EX int approx_distance(heptagon *h1, heptagon *h2) { - auto m = (hrmap_sol*) currentmap; - dynamicval g(geometry, gBinary4); - dynamicval cm(currentmap, m->binary_map); - int d1 = binary::celldistance3_approx(m->coords[h1].first, m->coords[h2].first); - int d2 = binary::celldistance3_approx(m->coords[h1].second, m->coords[h2].second); - return d1 + d2 - abs(h1->distance - h2->distance); - } - - EX string solshader = + EX string common = "uniform mediump sampler3D tInvExpTable;" "uniform mediump float PRECX, PRECY, PRECZ;" @@ -312,8 +377,55 @@ EX namespace solv { " y /= (1.+z);" " return 0.5 - atan((0.5-x) / y) / 3.1415926535897932384626433832795;" - " }" + " }"; + +EX } +EX namespace solv { + + EX ld solrange_xy = 15; + EX ld solrange_z = 4; + + EX bool in_table_range(hyperpoint h) { + return abs(h[0]) < solrange_xy && abs(h[1]) < solrange_xy && abs(h[2]) < solrange_z; + } + + EX solnihv::tabled_inverses solt = solnihv::tabled_inverses("solv-geodesics.dat"); + + hyperpoint christoffel(const hyperpoint at, const hyperpoint velocity, const hyperpoint transported) { + return hpxyz3( + -velocity[2] * transported[0] - velocity[0] * transported[2], + velocity[2] * transported[1] + velocity[1] * transported[2], + velocity[0] * transported[0] * exp(2*at[2]) - velocity[1] * transported[1] * exp(-2*at[2]), + 0 + ); + } + + EX hyperpoint get_inverse_exp(hyperpoint h, bool lazy, bool just_direction) { + solt.load(); + + ld ix = h[0] >= 0. ? solnihv::x_to_ix(h[0]) : solnihv::x_to_ix(-h[0]); + ld iy = h[1] >= 0. ? solnihv::x_to_ix(h[1]) : solnihv::x_to_ix(-h[1]); + ld iz = tanh(h[2]); + + if(h[2] < 0.) { iz = -iz; swap(ix, iy); } + + hyperpoint res = solt.get(ix, iy, iz, lazy); + + 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) { + ld r = hypot_d(3, res); + if(r == 0.) return res; + return res * atanh(r) / r; + } + + return res; + } + + EX string solshader = solnihv::common + "vec4 inverse_exp(vec4 h) {" @@ -343,9 +455,81 @@ EX namespace solv { "return res;" "}"; - #endif + + EX int approx_distance(heptagon *h1, heptagon *h2) { + auto m = (solnihv::hrmap_solnih*) currentmap; + dynamicval g(geometry, gBinary4); + dynamicval cm(currentmap, m->binary_map); + int d1 = binary::celldistance3_approx(m->coords[h1].first, m->coords[h2].first); + int d2 = binary::celldistance3_approx(m->coords[h1].second, m->coords[h2].second); + return d1 + d2 - abs(h1->distance - h2->distance); + } EX } +EX namespace nihv { + + EX solnihv::tabled_inverses niht = solnihv::tabled_inverses("h23-geodesics.dat"); + + hyperpoint christoffel(const hyperpoint at, const hyperpoint velocity, const hyperpoint transported) { + static const ld l2 = log(2); + static const ld l3 = log(3); + return hpxyz3( + (velocity[2] * transported[0] + velocity[0] * transported[2]) * l2, + (velocity[2] * transported[1] + velocity[1] * transported[2]) * l3, + -(velocity[0] * transported[0] * exp(-2*l2*at[2]) * l2 + velocity[1] * transported[1] * exp(-2*l3*at[2]) * l3), + +// (-velocity[2] * transported[0] - velocity[0] * transported[2]) * l2, +// (-velocity[2] * transported[1] - velocity[1] * transported[2]) * l3, +// velocity[0] * transported[0] * exp(2*l2*at[2]) * l2 + velocity[1] * transported[1] * exp(2*l3*at[2]) * l3, + 0 + ); + } + + EX hyperpoint get_inverse_exp(hyperpoint h, bool lazy, bool just_direction) { + niht.load(); + + ld ix = h[0] >= 0. ? solnihv::x_to_ix(h[0]) : solnihv::x_to_ix(-h[0]); + ld iy = h[1] >= 0. ? solnihv::x_to_ix(h[1]) : solnihv::x_to_ix(-h[1]); + ld iz = (tanh(h[2]/4)+1)/2; + + hyperpoint res = niht.get(ix, iy, iz, lazy); + + if(h[0] < 0.) res[0] = -res[0]; + if(h[1] < 0.) res[1] = -res[1]; + + if(!just_direction) { + ld r = hypot_d(3, res); + if(r == 0.) return res; + return res * atanh(r) / r; + } + + return res; + } + + EX string nihshader = solnihv::common + + + "vec4 inverse_exp(vec4 h) {" + + "float ix = h[0] >= 0. ? x_to_ix(h[0]) : x_to_ix(-h[0]);" + "float iy = h[1] >= 0. ? x_to_ix(h[1]) : x_to_ix(-h[1]);" + "float iz = (tanh(h[2]/4.)+1.) / 2.;" + + "vec4 res;" + + "float cx = ix*(1.-1./PRECX) + .5/PRECX;" + "float cy = iy*(1.-1./PRECY) + .5/PRECY;" + "float cz = iz*(1.-1./PRECZ) + .5/PRECZ;" + + "res = texture3D(tInvExpTable, vec3(cx, cy, cz));" + + "if(h[0] < 0.) res[0] = -res[0];" + "if(h[1] < 0.) res[1] = -res[1];" + + "return res;" + "}"; +EX } +#endif + EX namespace nilv { hyperpoint christoffel(const hyperpoint Position, const hyperpoint Velocity, const hyperpoint Transported) { @@ -1330,6 +1514,7 @@ EX namespace nisot { if(nil) return nilv::christoffel(at, velocity, transported); #if CAP_SOLV else if(sol) return solv::christoffel(at, velocity, transported); + else if(nih) return nihv::christoffel(at, velocity, transported); #endif else if(sl2) return slr::christoffel(at, velocity, transported); else return point3(0, 0, 0); @@ -1337,7 +1522,7 @@ EX namespace nisot { EX bool in_table_range(hyperpoint h) { #if CAP_SOLV - if(sol) return solv::in_table_range(h); + if(solnih) return solv::in_table_range(h); #endif return true; } @@ -1428,7 +1613,7 @@ EX namespace nisot { EX hrmap *new_map() { #if CAP_SOLV - if(sol) return new solv::hrmap_sol; + if(solnih) return new solnihv::hrmap_solnih; #endif if(nil) return new nilv::hrmap_nil; if(prod) return new product::hrmap_product; @@ -1451,7 +1636,11 @@ EX namespace nisot { } #if CAP_SOLV else if(argis("-fsol")) { - shift(); solv::solfname = args(); + shift(); solv::solt.fname = args(); + return 0; + } + else if(argis("-nihsol")) { + shift(); nihv::niht.fname = args(); return 0; } #endif diff --git a/polygons.cpp b/polygons.cpp index 46fd2e37..f182d800 100644 --- a/polygons.cpp +++ b/polygons.cpp @@ -758,7 +758,7 @@ void geometry_information::make_wall(int id, vector vertices, vector h = zshift(normalize_flat(h), center_altitude * (1-x-y) + altitudes[a] * x + altitudes[(a+1)%n] * y); hpcpush(h); return; } - if(sol || !binarytiling) { hpcpush(normalize(h)); return; } + if(solnih || !binarytiling) { hpcpush(normalize(h)); return; } hyperpoint res = binary::parabolic3(h[0], h[1]) * xpush0(yy*h[2]); hpcpush(res); }); @@ -776,7 +776,7 @@ void geometry_information::make_wall(int id, vector vertices, vector } if(nil) h = nilv::on_geodesic(vertices[a], vertices[(a+1)%n], y * 1. / STEP); - if(sol || !binarytiling) { hpcpush(normalize(h)); continue; } + if(solnih || !binarytiling) { hpcpush(normalize(h)); continue; } hyperpoint res = binary::parabolic3(h[0], h[1]) * xpush0(yy*h[2]); hpcpush(res); } @@ -993,6 +993,26 @@ void geometry_information::create_wall3d() { make_wall(6, {pt(-1,+1,+1), pt(+1,+1,+1), pt(+1,00,+1), pt(-1,00,+1)}); make_wall(7, {pt(-1,00,+1), pt(+1,00,+1), pt(+1,-1,+1), pt(-1,-1,+1)}); } + + if(geometry == gNIH) { + ld zstep = .5; + ld bwh = vid.binary_width / 6; + auto pt = [&] (int x, int y, int z) { return xpush(bwh*x) * ypush(bwh*y) * zpush(zstep*z) * C0; }; + println(hlog, xpush(1) * ypush(1) * zpush(1) * xpush(-2) * ypush(-3) * zpush(-1) * C0); + println(hlog, xpush(1) * ypush(1) * zpush(-1) * xpush(-2) * ypush(-3) * zpush(1) * C0); + make_wall(0, {pt(+3,-3,-1), pt(+3,-3,+1), pt(+3,+3,+1), pt(+3,+3,-1) }); + make_wall(1, {pt(-3,+3,-1), pt(-3,+3,+1), pt(+3,+3,+1), pt(+3,+3,-1)}); + make_wall(2, {pt(-3,-3,-1), pt(-3,-3,+1), pt(-3,+3,+1), pt(-3,+3,-1) }); + make_wall(3, {pt(-3,-3,-1), pt(-3,-3,+1), pt(+3,-3,+1), pt(+3,-3,-1)}); + + make_wall(4, {pt(-3,-3,+1), pt(-3,+3,+1), pt(+3,+3,+1), pt(+3,-3,+1)}); + + for(int i=0; i<6; i++) { + int x = -3 + (i%2) * 3; + int y = -3 + (i/2) * 2; + make_wall(5+i, {pt(x,y,-1), pt(x+3,y,-1), pt(x+3,y+2,-1), pt(x,y+2,-1)}); + } + } if(geometry == gNil) { for(int i=0; i