non-isotropic hyperbolic space

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;
+    GLuint invexpid = tab.get_texture_id();
-    
-    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<PRECZ*PRECY*PRECX; z++) {
-        auto& t = inverse_exp_table[z];
-        xbuffer[z] = glhr::makevertex(t[0], t[1], t[2]);
-        }
-      println(hlog, "maxd = ", maxd);
-      
-      #if !ISWEB
-      glTexImage3D(GL_TEXTURE_3D, 0, 34836 /*GL_RGBA32F*/, PRECX, PRECX, PRECZ, 0, GL_RGBA, GL_FLOAT, xbuffer);
-      #else
-      // glTexStorage3D(GL_TEXTURE_3D, 1, 34836 /*GL_RGBA32F*/, PRECX, PRECX, PRECZ);
-      // glTexSubImage3D(GL_TEXTURE_3D, 0, 0, 0, 0, PRECX, PRECY, PRECZ, GL_RGBA, GL_FLOAT, xbuffer);
-      #endif
-      delete[] xbuffer;
-      
-      }
     
     glActiveTexture(GL_TEXTURE0 + glhr::INVERSE_EXP_BINDING);
     glBindTexture(GL_TEXTURE_3D, invexpid);
 
     glActiveTexture(GL_TEXTURE0 + 0);
     
-    glhr::set_solv_prec(solv::PRECX, solv::PRECY, solv::PRECZ);    
+    glhr::set_solv_prec(tab.PRECX, tab.PRECY, tab.PRECZ);
     }
   #endif
 

binary-tiling.cpp

@@ -51,7 +51,7 @@ EX namespace binary {
     }
   #endif
   
-  void breakhere() {
+  EX void breakhere() {
     exit(1);
     }
 

classes.cpp

@@ -529,6 +529,8 @@ EX geometryinfo1 giNil     = { gcNil,        3, 3, 4, {1,1, 1,0 } };
 EX geometryinfo1 giProduct = { gcSL2,        3, 3, 4, {1,1, 1,0 } /* will be filled in product::configure() */ };
 EX geometryinfo1 giSL2     = { gcSL2,        3, 3, 4, {1,1,-1,-1} };
 
+EX geometryinfo1 giH23     = { gcNIH,        3, 3, 4, {1,1, 1,0 } };
+
 /** list of available geometries */
 vector<geometryinfo> ginf = {
   {"{7,3}", "none",     "{7,3} (standard HyperRogue map)",            "HR",       7, 3, 0,         giHyperb2,       0, {{7, 5}}, eVariation::bitruncated},
@@ -588,6 +590,7 @@ vector<geometryinfo> 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

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

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<float, 3> 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<class T> void parallelize(int threads, int Nmin, int Nmax, T action) {
   std::vector<std::thread> v;
@@ -30,35 +33,13 @@ template<class T> void parallelize(int threads, int Nmin, int Nmax, T action) {
   for(std::thread& t:v) t.join();
   }
 
-hyperpoint sol1(pt v) {
+ld solerror(hyperpoint ok, hyperpoint chk) {
-  auto [x,y,z,t] = (array<ld,4>&) v;
+  auto zok  = point3( x_to_ix(ok[0]), x_to_ix(ok[1]), z_to_iz(ok[2]) );
-  if(x == 0 && z == 0) return C0;
+  auto zchk = point3( x_to_ix(chk[0]), x_to_ix(chk[1]), z_to_iz(chk[2]) );
-  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<ld,4>&) 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]) );
   return hypot_d(3, zok - zchk);
   }
 
-ld eucerror(pt ok, pt chk) {
+hyperpoint iterative_solve(hyperpoint xp, hyperpoint candidate, int prec, ld minerr, bool debug = false) {
-  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) {
 
   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);
+  return nih ? atanh(z * 2 - 1)*4 : atanh(z); // atanh(z * 2 - 1);
-  }
-
-ld z_to_iz(ld z) {
-  return tanh(z); // (tanh(z) + 1) / 2;
   }
 
 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<pt> candidates;
+      vector<hyperpoint> candidates;
-      pt cand;
+      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<pt> solved_candidates;
+      vector<hyperpoint> 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<last_x; x++)
   for(int y=0; y<last_y; y++) {
     for(int z=last_z; z<_PRECZ; z++)
       solution[z][y][x] = solution[z-1][y][x] * 2 - solution[z-2][y][x];
+    if(nih)
+      solution[0][y][x] = solution[1][y][x] * 2 - solution[2][y][x];
     }
-
+  
   for(int x=0; x<last_x; x++)
   for(int y=last_y; y<_PRECY; y++)
   for(int z=0; z<_PRECZ; z++)
@@ -322,10 +302,16 @@ int main(int argc, char **argv) {
 
   println(hlog);
   
+  /*
   geometry = gSol;
 
   build_sols();
   write_table("solv-geodesics-generated.dat");
+  */
+  
+  geometry = gNIH;
+  build_sols();
+  write_table("h23-geodesics-generated.dat");
 
   exit(0);
   }

geom-exp.cpp

@@ -517,6 +517,9 @@ EX string geometry_name() {
   switch(ginf[geometry].cclass) {
     case gcHyperbolic:
       return XLAT("hyperbolic");
+
+    case gcNIH:
+      return XLAT("nonisotropic hyperbolic");
     
     case gcEuclid: 
       return XLAT("flat");

geometry.cpp

@@ -545,7 +545,7 @@ void geometry_information::prepare_basics() {
   #endif
   #if CAP_BT
   if(binarytiling) hexvdist = rhexf = 1, tessf = 1, scalefactor = 1, crossf = hcrossf7;
-  if(geometry == gHoroRec || penrose || sol || nil) hexvdist = rhexf = .5, tessf = .5, scalefactor = .5, crossf = hcrossf7/2;
+  if(geometry == gHoroRec || penrose || sol || nil || nih) hexvdist = rhexf = .5, tessf = .5, scalefactor = .5, crossf = hcrossf7/2;
   #endif
   #if CAP_BT && MAXMDIM >= 4
   if(binarytiling) binary::build_tmatrix();

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

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)

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

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;

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;

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
-  #if CAP_SOLV
+EX namespace solnihv {
-  EX int PRECX, PRECY, PRECZ;
+
+  #if HDR
+  struct tabled_inverses {
+    int PRECX, PRECY, PRECZ;
+    vector<nisot::ptlow> tab;
+    string fname;
+    bool loaded;
+    
+    void load();
+    hyperpoint get(ld ix, ld iy, ld iz, bool lazy);
   
-  EX vector<nisot::ptlow> 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";
+  void tabled_inverses::load() {
-  
+    if(loaded) return;
-  EX void load_table() {
+    FILE *f = fopen(fname.c_str(), "rb");
-    if(table_loaded) return;
-    FILE *f = fopen(solfname.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);
+    tab.resize(PRECX * PRECY * PRECZ);
-    ignore(fread(&inverse_exp_table[0], sizeof(nisot::ptlow) * PRECX * PRECY * PRECZ, 1, f));
+    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) {
+  hyperpoint tabled_inverses::get(ld ix, ld iy, ld iz, bool lazy) {
-    return hpxyz3(
+    ix *= PRECX-1;
-      -velocity[2] * transported[0] - velocity[0] * transported[2],
+    iy *= PRECY-1;
-       velocity[2] * transported[1] + velocity[1] * transported[2],
+    iz *= PRECZ-1;
-       velocity[0] * transported[0] * exp(2*at[2]) - velocity[1] * transported[1] * exp(-2*at[2]),
+    
-       0
+    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<PRECZ*PRECY*PRECX; z++) {
+      auto& t = tab[z];
+      xbuffer[z] = glhr::makevertex(t[0], t[1], t[2]);
+      }
+    
+    #if !ISWEB
+    glTexImage3D(GL_TEXTURE_3D, 0, 34836 /*GL_RGBA32F*/, PRECX, PRECX, PRECZ, 0, GL_RGBA, GL_FLOAT, xbuffer);
+    #else
+    // glTexStorage3D(GL_TEXTURE_3D, 1, 34836 /*GL_RGBA32F*/, PRECX, PRECX, PRECZ);
+    // glTexSubImage3D(GL_TEXTURE_3D, 0, 0, 0, 0, PRECX, PRECY, PRECZ, GL_RGBA, GL_FLOAT, xbuffer);
+    #endif
+    delete[] xbuffer;    
+    return texture_id;
     }
   
   EX ld x_to_ix(ld u) {
@@ -83,59 +160,9 @@ EX namespace solv {
     return 0.5 - atan((0.5-x) / y) / M_PI;
     }
 
-  EX hyperpoint get_inverse_exp(hyperpoint h, bool lazy, bool just_direction) {
+  struct hrmap_solnih : hrmap {
-    load_table();
-    
-    ld ix = h[0] >= 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 {
     hrmap *binary_map;
+    hrmap *ternary_map; /* nih only */
     unordered_map<pair<heptagon*, heptagon*>, heptagon*> at;
     unordered_map<heptagon*, pair<heptagon*, heptagon*>> coords;
     
@@ -159,9 +186,10 @@ EX namespace solv {
       return h;      
       }
 
-    hrmap_sol() {
+    hrmap_solnih() {
     
       heptagon *alt;
+      heptagon *alt3;
     
       if(true) {
         dynamicval<eGeometry> 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<eGeometry> g(geometry, gTernary); 
+        alt3 = tailored_alloc<heptagon> (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<eGeometry> g(geometry, gTernary); 
+      dynamicval<hrmap*> 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);
+      if(sol) {
-      ld bw = vid.binary_width * z;
+        ld z = log(2);
-      ld bwh = bw / 4;
+        ld bw = vid.binary_width * z;
-      switch(i) {
+        ld bwh = bw / 4;
-        case 0: return xpush(+bw);
+        switch(i) {
-        case 1: return ypush(+bw);
+          case 0: return xpush(+bw);
-        case 2: return xpush(-bwh) * zpush(+z) * ypush(j == 6 ? +bwh : -bwh);
+          case 1: return ypush(+bw);
-        case 3: return xpush(+bwh) * zpush(+z) * ypush(j == 6 ? +bwh : -bwh);
+          case 2: return xpush(-bwh) * zpush(+z) * ypush(j == 6 ? +bwh : -bwh);
-        case 4: return xpush(-bw);
+          case 3: return xpush(+bwh) * zpush(+z) * ypush(j == 6 ? +bwh : -bwh);
-        case 5: return ypush(-bw);
+          case 4: return xpush(-bw);
-        case 6: return ypush(-bwh) * zpush(-z) * xpush(j == 2 ? +bwh : -bwh);
+          case 5: return ypush(-bw);
-        case 7: return ypush(+bwh) * zpush(-z) * xpush(j == 2 ? +bwh : -bwh);
+          case 6: return ypush(-bwh) * zpush(-z) * xpush(j == 2 ? +bwh : -bwh);
-        default:return Id;
+          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<heptagon*,heptagon*> 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 string common = 
-  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<eGeometry> g(geometry, gBinary4); 
-    dynamicval<hrmap*> 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 = 
     "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<eGeometry> g(geometry, gBinary4); 
+    dynamicval<hrmap*> 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

polygons.cpp

@@ -758,7 +758,7 @@ void geometry_information::make_wall(int id, vector<hyperpoint> 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<hyperpoint> 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<S7; i++) make_wall(i, nilv::facevertices[i]);

rug.cpp

@@ -999,7 +999,7 @@ bincode acd_bin(ld x) {
 
 bincode get_bincode(hyperpoint h) {
   switch(ginf[gwhere].cclass) {
-    case gcEuclid: case gcSol: case gcNil: case gcProduct: case gcSL2:
+    case gcEuclid: case gcSol: case gcNil: case gcProduct: case gcSL2: case gcNIH:
       return acd_bin(h[0]) + acd_bin(h[1]) * sY + acd_bin(h[2]) * sZ;
     case gcHyperbolic:
       return acd_bin(hypot_d(3, h));

shaders.cpp

@@ -54,7 +54,7 @@ glvertex pointtogl(const hyperpoint& t);
 enum class shader_projection { standard, band, halfplane, standardH3, standardR3, 
   standardS30, standardS31, standardS32, standardS33, 
   ball, halfplane3, band3, flatten, standardSolv, standardNil,
-  standardEH2, standardSL2,
+  standardEH2, standardSL2, standardNIH,
   MAX 
   };
 
@@ -689,6 +689,7 @@ void init() {
     bool hp = among(sp, shader_projection::halfplane, shader_projection::halfplane3);
     bool sh3 = (sp == shader_projection::standardH3);
     bool ssol = (sp == shader_projection::standardSolv);
+    bool snih = (sp == shader_projection::standardNIH);
     bool snil = (sp == shader_projection::standardNil);
     bool ssl2 = (sp == shader_projection::standardSL2);
     bool sr3 = (sp == shader_projection::standardR3);
@@ -699,13 +700,15 @@ void init() {
     bool seh2 = (sp == shader_projection::standardEH2);
     bool ss3 = ss30 || ss31 || ss32 || ss33;
     
-    bool s3 = (sh3 || sr3 || ss3 || ssol || snil || seh2 || ssl2);
+    bool s3 = (sh3 || sr3 || ss3 || ssol || snil || seh2 || ssl2 || snih);
     bool dim3 = s3 || among(sp, shader_projection::ball, shader_projection::halfplane3, shader_projection::band3);
     bool dim2 = !dim3;
     bool ball = (sp == shader_projection::ball);
     bool flatten = (sp == shader_projection::flatten);
     
-    if(ssol && !CAP_SOLV) continue;
+    bool sson = ssol || snih;
+    
+    if(sson && !CAP_SOLV) continue;
     
     programs[i][j] = new GLprogram(stringbuilder(
       1,       "#define PI 3.14159265358979324\n",
@@ -755,6 +758,7 @@ void init() {
       
       #if CAP_SOLV
       ssol,    solv::solshader,
+      snih,    nihv::nihshader,
       #endif
       snil,    nilv::nilshader,
       ssl2,    slr::slshader,
@@ -792,10 +796,10 @@ void init() {
       hp && dim3, "t.x /= -rads; t.y /= -rads; t.z /= -rads; t[3] = 1.0;",
       
       s3,        "vec4 t = uMV * aPosition;",
-      ssol,      "t = inverse_exp(t);",
+      sson,      "t = inverse_exp(t);",
-      ssol,      "float d = sqrt(t[0] * t[0] + t[1] * t[1] + t[2] * t[2]);",
+      sson,      "float d = sqrt(t[0] * t[0] + t[1] * t[1] + t[2] * t[2]);",
-      ssol,      "float ad = (d == 0.) ? 0. : (d < 1.) ? min(atanh(d), 10.) : 10.; ",
+      sson,      "float ad = (d == 0.) ? 0. : (d < 1.) ? min(atanh(d), 10.) : 10.; ",
-      ssol,      "float m = ad / d / 11.; t[0] *= m; t[1] *= m; t[2] *= m; ",
+      sson,      "float m = ad / d / 11.; t[0] *= m; t[1] *= m; t[2] *= m; ",
       snil,      "t = inverse_exp(t);",
       ssl2,      "t = inverse_exp(t);",
 
@@ -810,7 +814,7 @@ void init() {
 
       sh3,       "float fogs = (uFogBase - acosh(t[3]) / uFog);",
       sr3||snil||ssl2, "float fogs = (uFogBase - sqrt(t[0]*t[0] + t[1]*t[1] + t[2]*t[2]) / uFog);",
-      ssol,      "float fogs = (uFogBase - ad / uFog);",
+      sson,      "float fogs = (uFogBase - ad / uFog);",
 
       seh2,      "float fogs = (uFogBase - sqrt(z*z+d*d) / uFog);",
       
@@ -821,7 +825,7 @@ void init() {
 
       s3,        "vColor.xyz = vColor.xyz * fogs + uFogColor.xyz * (1.0-fogs);",
 
-      sh3 || sr3 || ssol || ball || seh2,"t[3] = 1.0;",
+      sh3 || sr3 || sson || ball || seh2,"t[3] = 1.0;",
       
       band || hp || s3 || ball,"gl_Position = uP * t;",
       dim3 && !s3, "vColor.xyz = vColor.xyz * (0.5 - gl_Position.z / 2.0) + uFogColor.xyz * (0.5 + gl_Position.z / 2.0);",