diff --git a/devmods/solv-table.cpp b/devmods/solv-table.cpp
index e0f5e79c..7663bfdd 100644
--- a/devmods/solv-table.cpp
+++ b/devmods/solv-table.cpp
@@ -26,15 +26,7 @@ namespace hr {
 
 transmatrix parabolic1(ld u);
 
-namespace nisot {
-
-typedef hyperpoint pt;
-
-using sn::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])}); }
+namespace sn {
 
 template<class T> void parallelize(int threads, int Nmin, int Nmax, T action) {
   std::vector<std::thread> v;
@@ -143,49 +135,11 @@ hyperpoint iterative_solve(hyperpoint xp, hyperpoint candidate, int prec, ld min
   return at;
   }
 
-ptlow mlow(ld x, ld y, ld z) { return ptlow({float(x), float(y), float(z)}); }
+using ptlow = compressed_point;
 
-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(hyperpoint x) {
-  // azimuthal equidistant to Poincare
-  ld r = sqrt(x[0] * x[0] + x[1] * x[1] + x[2] * x[2]);
-  if(r == 0) return mlow(0,0,0);
-  ld make_r = sinh(r) / (1 + cosh(r));
-  ld d = make_r / r;
-  return mlow(x[0]*d, x[1]*d, x[2]*d);
-  }
-
-hyperpoint uncan(ptlow x) {
-  // Poincare to azimuthal equidistant
-
-  ld hr = x[0] * x[0] + x[1] * x[1] + x[2] * x[2];
-  if(hr == 0) return atxyz(0,0,0);
-  
-  ld hz = (1 + hr) / (1 - hr);
-  // println(hlog, "hz = ", hz);
-  ld d = (hz+1) * acosh(hz) / sinh(acosh(hz));
-  // println(hlog, "d = ", hz);
-  
-  // println(hlog, "returning: ", point3(x[0] * d, x[1] * d, x[2] * d));
-    
-  return point3(x[0] * d, x[1] * d, x[2] * d);
-  }
-
-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);
-  ld make_r = atanh(r);
-  println(hlog, "make_r = ", make_r);
-  ld d = make_r / r;
-  println(hlog, "d = ", d);
-  return atxyz(x[0]*d, x[1]*d, x[2]*d);
-  }
+ptlow operator +(ptlow a, ptlow b) { return make_array<float>(a[0]+b[0], a[1]+b[1], a[2]+b[2]); }
+ptlow operator -(ptlow a, ptlow b) { return make_array<float>(a[0]-b[0], a[1]-b[1], a[2]-b[2]); }
+ptlow operator *(ptlow a, ld x) { return make_array<float>(a[0]*x, a[1]*x, a[2]*x); }
 
 void fint(FILE *f, int x) { fwrite(&x, sizeof(x), 1, f); }
 void ffloat(FILE *f, float x) { fwrite(&x, sizeof(x), 1, f); }
@@ -206,28 +160,10 @@ void alloc_table(sn::tabled_inverses& tab, int X, int Y, int Z) {
   tab.tab.resize(X*Y*Z);
   }
 
-ld ix_to_x(ld ix) {
-  ld minx = 0;
-  while(x_to_ix(minx) <= ix) minx++;
-  ld maxx = minx; minx--;
-  for(int it=0; it<20; it++) {
-    ld x = (minx + maxx) / 2;
-    if(x_to_ix(x) < ix) minx = x;
-    else maxx = x;
-    }
-  return minx;
-  }
-
-ld iz_to_z(ld z) {
-  return nih ? atanh(z * 2 - 1)*4 : atanh(z); // atanh(z * 2 - 1);
-  }
-
 ld ptd(ptlow p) {
   return p[0]*p[0] + p[1]*p[1] + p[2] * p[2];
   }
 
-ptlow zflip(ptlow x) { return mlow(x[1], x[0], -x[2]); }
-
 void fix_boundaries(sn::tabled_inverses& tab, int last_x, int last_y, int last_z) {
   int PRECX = tab.PRECX;
   int PRECY = tab.PRECY;
@@ -270,7 +206,7 @@ void build_sols(int PRECX, int PRECY, int PRECZ) {
       vector<hyperpoint> candidates;
       hyperpoint cand;
       
-      candidates.push_back(atxyz(0,0,0)); 
+      candidates.push_back(point3(0,0,0)); 
       
       static constexpr int prec = 100;
       
@@ -301,23 +237,15 @@ void build_sols(int PRECX, int PRECY, int PRECZ) {
         println(hlog, "f(?) = ", v);
         println(hlog, "f(", cand, ") = ", nisot::numerical_exp(cand, prec));
         println(hlog, "error = ", xerr);
-        println(hlog, "canned = ", can(cand));
+        println(hlog, "canned = ", compress(azeq_to_table(cand)));
         max_err = xerr;
-        /*
-        hyperpoint h1 = uncan(solution[iz][iy-1][ix]);
-        hyperpoint h2 = uncan(solution[iz][iy][ix-1]);
-        hyperpoint h3 = uncan(solution[iz][iy-1][ix-1]);
-        hyperpoint h4 = h1 + h2 - h3;
-        solution[iz][iy][ix] = can(h4);
-        */
         return;
         }
 
       auto& so = tab.get_int(ix, iy, iz);
+      
+      so = compress(azeq_to_table(cand));
 
-      so = can(cand);
-      
-      
       for(int z=0; z<3; z++) if(isnan(so[z]) || isinf(so[z])) {
         println(hlog, cand, "canned to ", so);
         exit(4);
@@ -486,7 +414,7 @@ hyperpoint find_optimal_geodesic(hyperpoint res) {
       max_iter = 1000;
       auto h1 = iterative_solve(res, p.first * quality(p) / hypot_d(3, p.first), 100, 1e-6);
 
-      if(deb) println(hlog, "h1 returns ", h1, " of length ", hypot_d(3, h1), " and error ", hypot_d(3, numerical_exp(h1, 100) - res));
+      if(deb) println(hlog, "h1 returns ", h1, " of length ", hypot_d(3, h1), " and error ", hypot_d(3, nisot::numerical_exp(h1, 100) - res));
 
       if(h1 == fail) return;
       
@@ -509,7 +437,7 @@ hyperpoint find_optimal_geodesic(hyperpoint res) {
   }
 
 void fix_bugs(sn::tabled_inverses& tab) {
-  auto bug = can(fail);
+  auto bug = compress(azeq_to_table(fail));
   for(int iz=0; iz<tab.PRECZ; iz++) 
     for(int iy=0; iy<tab.PRECY; iy++)
       for(int ix=0; ix<tab.PRECX; ix++)  {
@@ -565,9 +493,7 @@ void improve(sn::tabled_inverses& tab) {
       
       if(h2 != fail) {
         auto& so = tab.get_int(ix, iy, iz);
-        so = can(h2);        
-        // println(hlog, "can: ", can(h2), " @ ", so[0]*so[0]+so[1]*so[1]+so[2]*so[2]);
-        // println(hlog, h2, " vs ", uncan(can(h2)));
+        so = compress(azeq_to_table(h2)); 
         }
       }
     };
diff --git a/nonisotropic.cpp b/nonisotropic.cpp
index a45fb1cf..6a6b1a16 100644
--- a/nonisotropic.cpp
+++ b/nonisotropic.cpp
@@ -9,9 +9,6 @@
 namespace hr {
 
 EX namespace nisot {
-  #if HDR
-  typedef array<float, 3> ptlow;
-  #endif
 
   #if HDR
   inline bool local_perspective_used() { return nonisotropic || prod; }
@@ -60,16 +57,21 @@ EX namespace sn {
     }
 
   #if HDR
+  typedef array<float, 3> compressed_point;
+  
+  inline hyperpoint decompress(compressed_point p) { return point3(p[0], p[1], p[2]); }
+  inline compressed_point compress(hyperpoint h) { return make_array<float>(h[0], h[1], h[2]); }
+
   struct tabled_inverses {
     int PRECX, PRECY, PRECZ;
-    vector<nisot::ptlow> tab;
+    vector<compressed_point> tab;
     string fname;
     bool loaded;
     
     void load();
     hyperpoint get(ld ix, ld iy, ld iz, bool lazy);
     
-    nisot::ptlow& get_int(int ix, int iy, int iz) { return tab[(iz*PRECY+iy)*PRECX+ix]; }
+    compressed_point& get_int(int ix, int iy, int iz) { return tab[(iz*PRECY+iy)*PRECX+ix]; }
   
     GLuint texture_id;
     bool toload;
@@ -89,7 +91,7 @@ EX namespace sn {
     ignore(fread(&PRECY, 4, 1, f));
     ignore(fread(&PRECZ, 4, 1, f));
     tab.resize(PRECX * PRECY * PRECZ);
-    ignore(fread(&tab[0], sizeof(nisot::ptlow) * PRECX * PRECY * PRECZ, 1, f));
+    ignore(fread(&tab[0], sizeof(compressed_point) * PRECX * PRECY * PRECZ, 1, f));
     fclose(f);
     loaded = true;    
     }
@@ -102,8 +104,7 @@ EX namespace sn {
     hyperpoint res;
     
     if(lazy) {
-      auto r = get_int(int(ix), int(iy), int(iz));
-      for(int i=0; i<3; i++) res[i] = r[i];
+      return decompress(get_int(int(ix), int(iy), int(iz)));
       }
     
     else {
@@ -180,6 +181,56 @@ EX namespace sn {
     return 0.5 - atan((0.5-x) / y) / M_PI;
     }
 
+  EX ld ix_to_x(ld ix) {
+    ld minx = 0;
+    while(x_to_ix(minx) <= ix) minx++;
+    ld maxx = minx; minx--;
+    for(int it=0; it<20; it++) {
+      ld x = (minx + maxx) / 2;
+      if(x_to_ix(x) < ix) minx = x;
+      else maxx = x;
+      }
+    return minx;
+    }
+
+  EX ld z_to_iz(ld z) { 
+    z = sinh(z) / (1 + cosh(z));
+    if(nih) z = (z+1) / 2;
+    return z;
+    }
+
+  EX ld iz_to_z(ld iz) {
+    ld minz = 0;
+    while(z_to_iz(minz) <= iz) minz++;
+    while(z_to_iz(minz) > iz) minz--;
+    ld maxz = minz + 1;
+    for(int it=0; it<20; it++) {
+      ld z = (minz + maxz) / 2;
+      if(z_to_iz(z) < iz) minz = z;
+      else maxz = z;
+      }
+    return (minz+maxz) / 2;
+    }
+
+  EX hyperpoint azeq_to_table(hyperpoint x) {
+    // azimuthal equidistant to Poincare
+    ld r = hypot_d(3, x);
+    if(r == 0) return point3(0,0,0);
+    ld make_r = sinh(r) / (1 + cosh(r));
+    ld d = make_r / r;
+    return x * d;
+    }
+
+  EX hyperpoint table_to_azeq(hyperpoint x) {
+    // Poincare to azimuthal equidistant  
+    ld hr = sqhypot_d(3, x);
+    if(hr == 0) return point3(0,0,0);    
+    ld hz = (1 + hr) / (1 - hr);
+    ld d = (hz+1) * acosh(hz) / sinh(acosh(hz));    
+    return x * d;
+    }
+
+
   struct hrmap_solnih : hrmap {
     hrmap *binary_map;
     hrmap *ternary_map; /* nih only */
@@ -459,7 +510,16 @@ EX namespace sn {
     "  y /= (1.+z);"
     
     "  return 0.5 - atan((0.5-x) / y) / 3.1415926535897932384626433832795;"
-    "  }";
+    "  }"
+    
+    "float z_to_iz_s(float z) {"
+      "return sinh(z) / (1 + cosh(z));"
+      "}"
+    
+    "float z_to_iz_ns(float z) {"
+      "z = sinh(z) / (1 + cosh(z));"
+      "return (z+1)/2;"
+      "}";    
 
   hyperpoint christoffel(const hyperpoint at, const hyperpoint velocity, const hyperpoint transported) {
     const ld l2 = log(2);
@@ -507,14 +567,7 @@ EX namespace sn {
     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 point3(0,0,0);
-      ld make_r = atanh(r);
-      if(r == 1) make_r = 30;
-      ld d = make_r / r;
-      return point3(res[0]*d, res[1]*d, res[2]*d);
-      }
+    if(!just_direction) return table_to_azeq(res);
 
     return res;
     }
@@ -532,11 +585,7 @@ EX namespace sn {
     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;
-      }
+    if(!just_direction) return table_to_azeq(res);
 
     return res;
     }
@@ -547,7 +596,7 @@ EX namespace sn {
     
     "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]);"
+    "float iz = z_to_iz_s(h[2]);"
     
     "if(h[2] < 1e-6) { iz = -iz; float s = ix; ix = iy; iy = s; }"
     "if(iz < 0.) iz = 0.;"
@@ -578,7 +627,7 @@ EX namespace sn {
       
     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.;
+    float iz = z_to_iz_ns(h[2]);
       
     vec4 res;