PSL(2,R)

3d-models.cpp

@@ -615,8 +615,8 @@ void geometry_information::animate_bird(hpcshape& orig, hpcshape_animated& anima
   // shift_shape(orig, BIRD);
   }
 
-EX hyperpoint forward_dir(ld x) { return prod ? point3(x, 0, 0) : xpush0(x); }
+EX hyperpoint forward_dir(ld x) { return (prod || sl2) ? point3(x, 0, 0) : xpush0(x); }
-EX hyperpoint zforward_dir(ld z) { return prod ? point3(0, 0, z) : zpush0(z); }
+EX hyperpoint zforward_dir(ld z) { return (prod || sl2) ? point3(0, 0, z) : zpush0(z); }
 
 EX hyperpoint dir_to_point(hyperpoint h) { return prod ? product::direct_exp(h) : h; }
 

basegraph.cpp

@@ -303,6 +303,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 && sl2)
+      shaderside_projection = true, glhr::new_shader_projection = glhr::shader_projection::standardSL2, pers3 = true;
     if(GDIM == 3 && apply_models && pmodel == mdGeodesic && nil)
       shaderside_projection = true, glhr::new_shader_projection = glhr::shader_projection::standardNil, pers3 = true;
     if(GDIM == 3 && sphere && apply_models && pmodel == mdPerspective) {

bigstuff.cpp

@@ -246,7 +246,7 @@ EX heptagon *createAlternateMap(cell *c, int rad, hstate firststate, int special
   alt->emeraldval = 0;
   alt->zebraval = 0;
   alt->distance = 0;
-  alt->fieldval = product::current_view_level;
+  alt->fieldval = nisot::current_view_level;
   alt->c7 = NULL;
   alt->alt = alt;
   h->alt = alt;

cell.cpp

@@ -148,6 +148,8 @@ EX cell *createMov(cell *c, int d) {
   if(c->move(d)) return c->move(d);
   else if(geometry == gProduct)
     product::find_cell_connection(c, d);
+  else if(sl2)
+    slr::find_cell_connection(c, d);
   #if CAP_BT
   else if(penrose)
     kite::find_cell_connection(c, d);

classes.cpp

@@ -526,6 +526,7 @@ geometryinfo1 giSphere3 = { gcSphere,     3, 3, 4, {1,1, 1,+1} };
 geometryinfo1 giSol     = { gcSol,        3, 3, 4, {1,1, 1,0 } };
 geometryinfo1 giNil     = { gcNil,        3, 3, 4, {1,1, 1,0 } };
 geometryinfo1 giProduct = { /* will be filled in product::configure() */ };
+geometryinfo1 giSL2     = { gcSL2,        3, 3, 4, {1,1,-1,-1} };
 
 /** list of available geometries */
 vector<geometryinfo> ginf = {
@@ -584,6 +585,7 @@ vector<geometryinfo> ginf = {
   {"kd3",    "none",    "kite-and-dart on horospheres",               "kd3",     12, 3, qsBP,      giHyperb3, 0x48200, {{7, 3}}, eVariation::pure},
   {"nil",    "none",    "Nil geometry",                               "nil",      8, 3, 0,         giNil,     0x48600, {{7, 5}}, eVariation::pure},
   {"product","none",    "product space",                              "product",  7, 3, 0,         giProduct, 0x48400, {{7, 3}}, eVariation::pure},
+  {"psl2",   "psl2",    "PSL(2,R)",                                   "psl2",     7, 3, qEXPERIMENTAL, giSL2, 0x49000, {{6, 4}}, eVariation::pure},
   };
 
   // bits: 9, 10, 15, 16, (reserved for later) 17, 18

classes.h

@@ -214,10 +214,10 @@ enum eGeometry {
   gHoroTris, gHoroRec, gHoroHex,
   gField435, gField534,
   gBinary4, gSol,
-  gKiteDart2, gKiteDart3, gNil, gProduct,
+  gKiteDart2, gKiteDart3, gNil, gProduct, gSL2,
   gGUARD};
 
-enum eGeometryClass { gcHyperbolic, gcEuclid, gcSphere, gcSol, gcNil, gcProduct };
+enum eGeometryClass { gcHyperbolic, gcEuclid, gcSphere, gcSol, gcNil, gcProduct, gcSL2 };
 
 enum class eVariation { bitruncated, pure, goldberg, irregular, dual };  
 

complex.cpp

@@ -2668,7 +2668,7 @@ EX namespace sword {
     };
   
   /** dimensions available to the Sword */
-  #define SWORDDIM (prod ? 2 : WDIM)
+  #define SWORDDIM ((prod || sl2) ? 2 : WDIM)
     
   #endif
 
@@ -2695,7 +2695,7 @@ EX namespace sword {
   
   cell *pos2(cell *c, int s) {
     int t = c->type;
-    if(prod) t -= 2;
+    if(prod || sl2) t -= 2;
     s *= 2;
     s += sword_angles/t;
     s %= (2 * sword_angles);
@@ -2748,7 +2748,7 @@ EX namespace sword {
     int s2 = neighborId(c2, c1);
     if(s1 < 0 || s2 < 0) return d;
     if(SWORDDIM == 2) {
-      int sub = prod ? 2 : 0;
+      int sub = (prod || sl2) ? 2 : 0;
       int t2 = c2->type - sub;
       int t1 = c1->type - sub;
       if(c1->c.mirror(s1))

config.cpp

@@ -520,6 +520,8 @@ EX void initConfig() {
       sightranges[i] = 2 * M_PI;
     else if(ginf[i].cclass == gcEuclid)
       sightranges[i] = 10;
+    else if(ginf[i].cclass == gcSL2)
+      sightranges[i] = 1.5;
     else
       sightranges[i] = 5;
     sightranges[gArchimedean] = 10;

control.cpp

@@ -284,6 +284,13 @@ typedef SDL_Event eventtype;
 
 EX bool smooth_scrolling = false;
 
+transmatrix zforward_push(ld z) {
+  if(!sl2) return zpush(z);
+  transmatrix T = Id;
+  T[2][3] = z;
+  return T;
+  }
+
 EX void handlePanning(int sym, int uni) {
   auto& LPe = nisot::local_perspective;
   if(mousepan && dual::split([=] { handlePanning(sym, uni); })) return;
@@ -300,10 +307,10 @@ EX void handlePanning(int sym, int uni) {
   
 #if !ISPANDORA
   if(sym == SDLK_END && GDIM == 3) { 
-    View = solmul(cpush(2, -0.2*shiftmul), LPe, View), didsomething = true, playermoved = false;
+    View = solmul(zforward_push(-0.2*shiftmul), LPe, View), didsomething = true, playermoved = false;
     }
   if(sym == SDLK_HOME && GDIM == 3) { 
-    View = solmul(cpush(2, +0.2*shiftmul), LPe, View), didsomething = true, playermoved = false;
+    View = solmul(zforward_push(+0.2*shiftmul), LPe, View), didsomething = true, playermoved = false;
     }
   if(sym == SDLK_RIGHT) { 
     if(history::on)

floorshapes.cpp

@@ -942,7 +942,8 @@ void draw_shape_for_texture(floorshape* sh) {
   hyperpoint v1 = hpxyz3(0.25, 0.25, 0, 0);
   hyperpoint v2 = hpxyz3(0.25, -0.25, 0, 0);
 
-  for(int a=0; a<MAX_EDGE; a++)
+  // SL2 needs 42
+  for(int a=0; a<max(MAX_EDGE, 42); a++)
     texture_order([&] (ld x, ld y) {
       hyperpoint h = center + v1 * x + v2 * y;
       hyperpoint inmodel;

geom-exp.cpp

@@ -821,6 +821,10 @@ EX void showEuclideanMenu() {
       dialog::addSelItem(XLAT("Curvature"), XLAT("Nil"), 0);
       break;
 
+    case gcSL2:
+      dialog::addSelItem(XLAT("Curvature"), XLAT("~SL2(R)~"), 0);
+      break;
+
     case gcProduct:
       dialog::addSelItem(XLAT("Curvature"), XLAT("Product"), 0);
       break;

geometry.cpp

@@ -399,7 +399,9 @@ hpcshape
  *  based on comparing these values to the values in the other geometry.
  */
 
-ld hcrossf7 = 0.620672, hexf7 = 0.378077, tessf7 = 1.090550, hexhexdist7 = 0.566256;
+#if HDR
+static const ld hcrossf7 = 0.620672, hexf7 = 0.378077, tessf7 = 1.090550, hexhexdist7 = 0.566256;
+#endif
 
 bool scale_used() { return (shmup::on && geometry == gNormal && BITRUNCATED) ? (cheater || autocheat) : true; }
 
@@ -468,6 +470,22 @@ void geometry_information::prepare_basics() {
     goto prod_finish;
     }
 
+  if(sl2) {
+    dynamicval<eGeometry> g(geometry, gNormal);
+    check_cgi();
+    cgi.prepare_basics();
+    rhexf = cgi.rhexf/2;
+    hexf = cgi.hexf/2;
+    crossf = cgi.crossf/2;
+    hcrossf = cgi.hcrossf/2;
+    tessf = cgi.tessf/2;
+    hexvdist = cgi.hexvdist/2;
+    hexhexdist = cgi.hexhexdist/2;
+    base_distlimit = cgi.base_distlimit-1;
+    cgip = this;
+    goto prod_finish;
+    }
+  
   if((sphere || hyperbolic) && WDIM == 3 && !binarytiling) {
     rhexf = hexf = 0.378077;
     crossf = hcrossf = 0.620672;

graph.cpp

@@ -4426,7 +4426,7 @@ color_t transcolor(cell *c, cell *c2, color_t wcol) {
 
 // how much should be the d-th wall darkened in 3D
 int get_darkval(cell *c, int d) {
-  if(prod) {
+  if(prod || sl2) {
     return d >= c->type - 2 ? 4 : 0;
     }
   const int darkval_hbt[9] = {0,2,2,0,6,6,8,8,0};
@@ -4518,7 +4518,7 @@ EX int noclipped;
 void make_clipping_planes() {
 #if MAXMDIM >= 4
   clipping_planes.clear();
-  if(PIU(sphere)) return;
+  if(PIU(sphere) || experimental) return;
   auto add_clipping_plane = [] (ld x1, ld y1, ld x2, ld y2) {
     ld z1 = 1, z2 = 1;
     hyperpoint sx = point3(y1 * z2 - y2 * z1, z1 * x2 - z2 * x1, x1 * y2 - x2 * y1);
@@ -6031,7 +6031,7 @@ EX void drawcell(cell *c, transmatrix V, int spinv, bool mirrored) {
           
           for(int a=0; a<c->type; a++)
             if(c->move(a) && !isWall3(c->move(a), dummy)) {
-              if(pmodel == mdPerspective && !sphere && !quotient && !penrose && !nonisotropic && !prod) {
+              if(pmodel == mdPerspective && !sphere && !quotient && !penrose && !nonisotropic && !prod && !experimental) {
                 if(a < 4 && among(geometry, gHoroTris, gBinary3) && celldistAlt(c) >= celldistAlt(viewcenter())) continue;
                 else if(a < 2 && among(geometry, gHoroRec) && celldistAlt(c) >= celldistAlt(viewcenter())) continue;
                 else if(c->move(a)->master->distance > c->master->distance && c->master->distance > viewctr.at->distance && !quotient) continue;
@@ -6364,8 +6364,8 @@ EX void drawcell(cell *c, transmatrix V, int spinv, bool mirrored) {
             floorShadow(c, V, SHADOW_SL * sl);
             for(int s=0; s<sl; s++) 
             forCellIdEx(c2, i, c) {
-              int sl2 = snakelevel(c2);
+              int sl_2 = snakelevel(c2);
-              if(s >= sl2)
+              if(s >= sl_2)
                 if(placeSidewall(c, i, SIDE_SLEV+s, V, getSnakelevColor(c, s, sl, fd, wcol))) break;
               }
             }
@@ -8043,7 +8043,8 @@ EX void drawBug(const cellwalker& cw, color_t col) {
 
 EX cell *viewcenter() {
   if(masterless) return centerover.at;
-  else if(prod) return product::get_at(viewctr.at->c7, product::current_view_level);
+  else if(prod) return product::get_at(viewctr.at->c7, nisot::current_view_level);
+  else if(sl2) return slr::get_at(viewctr.at->c7, nisot::current_view_level);
   else return viewctr.at->c7;
   }
 

hyper.h

@@ -125,9 +125,10 @@ void addMessage(string s, char spamtype = 0);
 #define sphere (cgclass == gcSphere)
 #define sol (cgclass == gcSol)
 #define nil (cgclass == gcNil)
+#define sl2 (cgclass == gcSL2)
 #define prod (cgclass == gcProduct)
 #define hyperbolic (cgclass == gcHyperbolic)
-#define nonisotropic (sol || nil)
+#define nonisotropic (sol || nil || sl2)
 #define translatable (euclid || nonisotropic)
 #define nonorientable (ginf[geometry].flags & qNONORIENTABLE)
 #define elliptic (ginf[geometry].flags & qELLIPTIC)

hyperpoint.cpp

@@ -355,7 +355,8 @@ EX ld hypot_d(int d, const hyperpoint& h) {
   }
   
 EX ld zlevel(const hyperpoint &h) {
-  if(translatable) return h[LDIM];
+  if(sl2) return sqrt(-intval(h, Hypc));
+  else if(translatable) return h[LDIM];
   else if(sphere) return sqrt(intval(h, Hypc));
   else if(prod) return log(sqrt(abs(intval(h, Hypc)))); /* abs works with both underlying spherical and hyperbolic */
   else return (h[LDIM] < 0 ? -1 : 1) * sqrt(-intval(h, Hypc));
@@ -444,6 +445,7 @@ EX transmatrix eupush(hyperpoint h) {
   }
 
 EX transmatrix eupush3(ld x, ld y, ld z) {
+  if(sl2) return slr::translate(slr::xyz_point(x, y, z));
   return eupush(point3(x, y, z));
   }
 
@@ -827,6 +829,11 @@ EX ld hdist0(const hyperpoint& mh) {
       auto d1 = product_decompose(mh);
       return hypot(PIU(hdist0(d1.second)), d1.first);
       }
+    case gcSL2: {
+      auto cosh_r = hypot(mh[2], mh[3]);
+      auto phi = atan2(mh[2], mh[3]);
+      return hypot(cosh_r < 1 ? 0 : acosh(cosh_r), phi);
+      }
     default:
       return hypot_d(GDIM, mh);
     }
@@ -1017,8 +1024,14 @@ EX transmatrix transpose(transmatrix T) {
   }
 
 #if HDR
+namespace slr { 
+  hyperpoint xyz_point(ld x, ld y, ld z); 
+  hyperpoint polar(ld r, ld theta, ld phi); 
+  }
+
 inline hyperpoint cpush0(int c, ld x) { 
   hyperpoint h = Hypc;
+  if(sl2) return slr::xyz_point(c==0?x:0, c==1?x:0, c==2?x:0);
   if(c == 2 && prod) {
     h[2] = exp(x);
     return h;
@@ -1029,6 +1042,7 @@ inline hyperpoint cpush0(int c, ld x) {
   }
 
 inline hyperpoint xspinpush0(ld alpha, ld x) { 
+  if(sl2) return slr::polar(x, -alpha, 0);
   hyperpoint h = Hypc;
   h[LDIM] = cos_auto(x);
   h[0] = sin_auto(x) * cos(alpha);

hypgraph.cpp

@@ -1383,12 +1383,29 @@ EX void optimizeview() {
   if(subscreens::split(optimizeview)) return;
   if(dual::split(optimizeview)) return;
 
+  if(sl2) {
+    cell *c = viewcenter();
+    cell *cbest = NULL;
+    ld best = hdist0(tC0(gmatrix[c]));
+    if(isnan(best)) return;
+    forCellEx(c2, c) {
+      ld quality = hdist0(tC0(gmatrix[c2]));
+      if(quality < best) best = quality, cbest = c2;
+      }
+    if(cbest) {
+      View = View * currentmap->relative_matrix(cbest, c, C0);
+      viewctr.at = cbest->master;
+      nisot::current_view_level = slr::get_where(cbest).second;
+      }
+    return;
+    }
+  
   if(prod) {
     ld z = zlevel(tC0(View));
     View = mscale(View, -z);
     product::in_underlying_map(optimizeview);
-    if(z > cgi.plevel / 2) { product::current_view_level--; z -= cgi.plevel; }
+    if(z > cgi.plevel / 2) { nisot::current_view_level--; z -= cgi.plevel; }
-    if(z < -cgi.plevel / 2) { product::current_view_level++; z += cgi.plevel; }
+    if(z < -cgi.plevel / 2) { nisot::current_view_level++; z += cgi.plevel; }
     View = mscale(View, z);
     return;
     }
@@ -1483,7 +1500,7 @@ EX void resetview() {
   else centerover = cwt;
   cwtV = View;
   nisot::local_perspective = Id;
-  if(prod) product::current_view_level = product::get_where(cwt.at).second;
+  if(prod || sl2) nisot::current_view_level = product::get_where(cwt.at).second;
   // SDL_LockSurface(s);
   // SDL_UnlockSurface(s);
   }
@@ -1970,7 +1987,7 @@ bool limited_generation(cell *c) {
 
 EX bool do_draw(cell *c, const transmatrix& T) {
 
-  if(product::pmap) return product::in_actual([&] { return do_draw(product::get_at(c, product::current_view_level), T); });
+  if(product::pmap) return product::in_actual([&] { return do_draw(product::get_at(c, nisot::current_view_level), T); });
   if(WDIM == 3) {
     if(cells_drawn > vid.cells_drawn_limit) return false;
     if(nil && pmodel == mdGeodesic) {
@@ -1978,7 +1995,7 @@ EX bool do_draw(cell *c, const transmatrix& T) {
       if(dist > sightranges[geometry] + (vid.sloppy_3d ? 0 : 0.9)) return false;
       if(dist <= extra_generation_distance && !limited_generation(c)) return false;
       }
-    else if(pmodel == mdGeodesic && !nil) {
+    else if(pmodel == mdGeodesic && sol) {
       if(!nisot::in_table_range(tC0(T))) return false;
       if(!limited_generation(c)) return false;
       }

nonisotropic.cpp

@@ -10,6 +10,7 @@ namespace hr {
 EX namespace nisot {
   typedef array<float, 3> ptlow;
 
+  EX int current_view_level;  
   EX transmatrix local_perspective;
   #if HDR
   inline bool local_perspective_used() { return nonisotropic || prod; }
@@ -18,6 +19,8 @@ EX namespace nisot {
   EX bool geodesic_movement = true;
 
   EX transmatrix translate(hyperpoint h) {
+    if(sl2)
+      return slr::translate(h);
     transmatrix T = Id;
     for(int i=0; i<GDIM; i++) T[i][LDIM] = h[i];
     if(sol) {
@@ -564,8 +567,6 @@ EX namespace product {
     ginf[gProduct].g.graphical_dimension++;
     }
   
-  EX int current_view_level;
-  
   hrmap *pmap;
   geometry_information *pcgip;
   
@@ -634,7 +635,7 @@ EX namespace product {
     if(sphere) gmatrix[c] = V; /* some computations need gmatrix0 for underlying geometry */
     bool s = sphere;
     in_actual([&] { 
-      cell *c0 = get_at(c, current_view_level);
+      cell *c0 = get_at(c, nisot::current_view_level);
       cwall_offset = wall_offset(c0);
       if(s) cwall_mask = (1<<c->type) - 1;
       else {
@@ -652,7 +653,7 @@ EX namespace product {
       for(int z=-max_z; z<=max_z; z++) {
         if(z == 0) cwall_mask ^= (2<<c->type);
         if(z == 1) cwall_mask ^= (1<<c->type);
-        cell *c1 = get_at(c, current_view_level+z);
+        cell *c1 = get_at(c, nisot::current_view_level+z);
         setdist(c1, 7, NULL);
         drawcell(c1, V * mscale(Id, cgi.plevel * z), spinv, mirrored); 
         }
@@ -769,11 +770,360 @@ EX namespace product {
   
 EX }
 
+EX namespace slr {
+
+  /* This implementation is based on:
+  // https://pdfs.semanticscholar.org/bf46/824df892593a1b6d1c84a5f99e90eece7c54.pdf
+  // However, to make it consistent with the conventions in HyperRogue,
+  // coordinates 0<->2 and 1<->3 are swapped,
+  // then coordinates 2<->3 are swapped
+  */
+
+  EX hyperpoint from_phigans(hyperpoint h) {
+    ld r = asinh(hypot_d(2, h));
+    ld x = h[0];
+    ld y = h[1];
+    ld z = h[2];
+    return hyperpoint(x * cos(z) + y * sin(z), y * cos(z) - x * sin(z), cosh(r) * sin(z), cosh(r) * cos(z));
+    }
+  
+  EX hyperpoint to_phigans(hyperpoint h) {
+    ld z = atan2(h[2], h[3]);
+    ld x = h[0];
+    ld y = h[1];
+    return point31(x * cos(z) - y * sin(z), y * cos(z) + x * sin(z), z);
+    }
+
+  /** in the 'phigans' model */
+  hyperpoint christoffel(const hyperpoint Position, const hyperpoint Velocity, const hyperpoint Transported) {
+    ld x = Position[0];
+    ld y = Position[1];
+    ld s = x*x + y*y + 1;
+    ld x2 = x * x;
+    ld y2 = y * y;
+    ld x4 = x2 * x2;
+    ld y4 = y2 * y2;
+    return point3(
+     + Velocity[ 0 ] * Transported[ 0 ] * (x*(4*x2*y2 + 4*y4 + 9*y2 + 1))
+     + Velocity[ 0 ] * Transported[ 1 ] * (-y*(4*x4 + 4*x2*y2 + 9*x2 + 2))
+     + Velocity[ 0 ] * Transported[ 2 ] * (-x*y*(x2 + y2 + 2))
+     + Velocity[ 1 ] * Transported[ 0 ] * (-y*(4*x4 + 4*x2*y2 + 9*x2 + 2))
+     + Velocity[ 1 ] * Transported[ 1 ] * (x*(4*x4 + 4*x2*y2 + 9*x2 + 5))
+     + Velocity[ 1 ] * Transported[ 2 ] * (x4 + x2*y2 + 2*x2 + 1)
+     + Velocity[ 2 ] * Transported[ 0 ] * (-x*y*(x2 + y2 + 2))
+     + Velocity[ 2 ] * Transported[ 1 ] * (x4 + x2*y2 + 2*x2 + 1),
+  
+     + Velocity[ 0 ] * Transported[ 0 ] *  (y*(4*x2*y2 + 4*y4 + 9*y2 + 5))
+     + Velocity[ 0 ] * Transported[ 1 ] *  (-x*(4*x2*y2 + 4*y4 + 9*y2 + 2))
+     + Velocity[ 0 ] * Transported[ 2 ] *  (-x2*y2 - y4 - 2*y2 - 1)
+     + Velocity[ 1 ] * Transported[ 0 ] *  (-x*(4*x2*y2 + 4*y4 + 9*y2 + 2))
+     + Velocity[ 1 ] * Transported[ 1 ] *  (y*(4*x4 + 4*x2*y2 + 9*x2 + 1))
+     + Velocity[ 1 ] * Transported[ 2 ] *  (x*y*(x2 + y2 + 2))
+     + Velocity[ 2 ] * Transported[ 0 ] *  (-x2*y2 - y4 - 2*y2 - 1)
+     + Velocity[ 2 ] * Transported[ 1 ] *  (x*y*(x2 + y2 + 2)),
+  
+     + Velocity[ 0 ] * Transported[ 0 ] *  (-4*x*y)
+     + Velocity[ 0 ] * Transported[ 1 ] *  (2*x2 - 2*y2)
+     + Velocity[ 0 ] * Transported[ 2 ] *  x
+     + Velocity[ 1 ] * Transported[ 0 ] *  (2*x2 - 2*y2)
+     + Velocity[ 1 ] * Transported[ 1 ] *  4*x*y
+     + Velocity[ 1 ] * Transported[ 2 ] *  y
+     + Velocity[ 2 ] * Transported[ 0 ] *  x
+     + Velocity[ 2 ] * Transported[ 1 ] *  y
+       ) / s;
+    }
+
+  EX transmatrix translate(hyperpoint h) {
+    return matrix4(
+      h[3], -h[2],  h[1],  h[0],
+      h[2],  h[3], -h[0],  h[1],
+      h[1], -h[0],  h[3],  h[2],
+      h[0],  h[1], -h[2],  h[3]
+      );
+    }
+  
+  EX hyperpoint polar(ld r, ld theta, ld phi) {
+    return hyperpoint(sinh(r) * cos(theta-phi), sinh(r) * sin(theta-phi), cosh(r) * sin(phi), cosh(r) * cos(phi));
+    }
+  
+  EX hyperpoint xyz_point(ld x, ld y, ld z) {
+    ld r = hypot(x, y);
+    ld f = r ? sinh(r) / r : 1;
+    return hyperpoint(x * f * cos(z) + y * f * sin(z), y * f * cos(z) - x * f * sin(z), cosh(r) * sin(z), cosh(r) * cos(z));
+    }
+  
+  ld rootsin(ld square, ld s) {
+    if(square > 0) return sinh(sqrt(square) * s) / sqrt(square);
+    else if(square < 0) return sin(sqrt(-square) * s) / sqrt(-square);
+    else return s;
+    }
+  
+  /** it==0 is standard asin, it==1 is the next solution (PI-asin) */
+  ld asin_it(ld z, int it) {
+    auto ans = asin(z);
+    if(it & 1) ans = M_PI - ans;
+    return ans;
+    }
+  
+  ld arootsin(ld square, ld v, int it) {
+    if(square > 0) return asinh(v * sqrt(square)) / sqrt(square);
+    else if(square < 0) return asin_it(v * sqrt(-square), it) / sqrt(-square);
+    else return v;
+    }
+  
+  ld roottan(ld square, ld s) {
+    if(square > 0) return tanh(sqrt(square) * s) / sqrt(square);
+    else if(square < 0) return tan(sqrt(-square) * s) / sqrt(-square);
+    else return s;
+    }
+
+  hyperpoint geodesic_polar(ld alpha, ld beta, ld s) {
+    auto c = cos(2*alpha);
+    
+    ld t;
+    if(c > 0) 
+      t = atan(sin(alpha) * tanh(sqrt(c) * s) / sqrt(c));
+    else if(c < 0) {
+      /* the formula in the paper is roughly atan(k*tan(s))
+       * however, atan is not always to be taken in [-PI/2,PI/2]:
+       * if s is in [kPI-PI/2, kPI+PI/2], we should also increase the result by kPI
+       */
+      ld x = sqrt(-c) * s;
+      ld steps = floor(x/M_PI + 0.5);
+      t = atan(sin(alpha) * tan(sqrt(-c) * s) / sqrt(-c)) + M_PI * steps;
+      }
+    else t = atan(sin(alpha) * s);
+    
+    return polar(
+      asinh(cos(alpha) * rootsin(c, s)),
+      beta - t,
+      2*sin(alpha)*s - t
+      );
+    }
+
+  EX hyperpoint formula_exp(hyperpoint h) {
+    ld s = hypot_d(3, h);
+    ld beta = atan2(h[1], h[0]);
+    ld alpha = asin(h[2] / s);
+    return geodesic_polar(alpha, beta, s);
+    }
+
+  void find_alpha(ld phi, ld r, ld theta, ld &alpha, ld &s, ld &beta) {
+    if(phi < 0) { find_alpha(-phi, r, -theta, alpha, s, beta); alpha = -alpha; beta = -beta; return; }
+    ld mina = 0, maxa = M_PI/2;
+  
+    bool next_nan;
+    ld c;
+    
+    for(int it=0; it<40; it++) {
+      alpha = (mina + maxa) / 2;
+  
+      c = cos(2 * alpha);    
+      s = arootsin(c, sinh(r) / cos(alpha), 0);
+      if(isnan(s)) { next_nan = true, maxa = alpha; continue; }
+      ld got_phi = 2*sin(alpha)*s - atan(sin(alpha) * roottan(c, s));
+      if(got_phi > phi) next_nan = false, maxa = alpha;
+      else mina = alpha;
+      }
+    
+    if(next_nan) {
+      mina = M_PI/4;
+  
+      for(int it=0; it<40; it++) {
+        alpha = (mina + maxa) / 2;
+        c = cos(2 * alpha);    
+        s = arootsin(c, sinh(r) / cos(alpha), 1);
+        ld got_phi = 2*sin(alpha)*s - atan(sin(alpha) * roottan(c, s)) - M_PI;
+        if(got_phi < phi) maxa = alpha;
+        else mina = alpha;
+        }
+      beta = theta + atan(sin(alpha) * roottan(c, s)) + M_PI;
+      }
+    else beta = theta + atan(sin(alpha) * roottan(c, s));
+    }
+
+  EX hyperpoint get_inverse_exp(hyperpoint h) {
+    if(sqhypot_d(2, h) < 1e-12) return point3(0, 0, atan2(h[2], h[3]));
+    ld r = asinh(hypot_d(2, h));
+    ld phi = atan2(h[2], h[3]);
+    ld theta = atan2(h[1], h[0]) + phi;
+    
+    ld alpha, s, beta;
+    find_alpha(phi, r, theta, alpha, s, beta);
+    
+    return point3(s * cos(beta) * cos(alpha), s * sin(beta) * cos(alpha), s * sin(alpha));
+    }
+
+  EX string slshader = 
+    "float atan2(float y, float x) {"
+    "  if(x == 0.) return y > 0. ? PI/2. : -PI/2.;"
+    "  if(x > 0.) return atan(y / x);"
+    "  if(y >= 0.) return atan(y / x) + PI;"
+    "  if(y < 0.) return atan(y / x) - PI;"
+    "  }"
+    
+    "vec4 inverse_exp(vec4 h) {"
+      "if(h[0]*h[0] + h[1] * h[1] < 1e-6) return vec4(0, 0, atan(h[2], h[3]), 1);"
+      "float r = asinh(sqrt(h[0] * h[0] + h[1] * h[1]));"
+      "float phi = atan2(h[2], h[3]);"
+      "float theta = atan2(h[1], h[0]) + phi;"
+      "float alpha;"
+      "float s;"
+      "float beta;"
+      "float sgn = 1.;"
+      "if(phi < 0.) { phi = -phi; theta = -theta; sgn = -1.; }"
+      "float c;"
+      "s = sinh(r) / cos(PI/4.);"
+      "float gphi = 2.*sin(PI/4.)*s - atan(sin(PI/4.) * s);"
+      "if(gphi > phi) {"
+      "  float mina = 0.;"
+      "  float maxa = PI/4.;"
+      "  for(int it=0; it<40; it++) {"
+      "    alpha = (mina + maxa) / 2.;"
+      "    c = sqrt(cos(2. * alpha));"
+      "    s = asinh(sinh(r) / cos(alpha) * c) / c;"
+      "    gphi = 2.*sin(alpha)*s - atan(sin(alpha) * tanh(c * s) / c);"
+      "    if(gphi > phi) maxa = alpha;"
+      "    else mina = alpha;"
+      "    }"
+      "  beta = theta + atan(sin(alpha) * tanh(c * s) / c);"
+      "  }"
+      "else {"
+      "  int next_nan = 1;"
+      "  float mina = PI/4.;"
+      "  float maxa = PI/2.;"
+      "  for(int it=0; it<40; it++) {"
+      "    alpha = (mina + maxa) / 2.;"
+      "    c = sqrt(-cos(2. * alpha));"
+      "    if(sinh(r) * c > cos(alpha)) { next_nan = 1; maxa = alpha; continue; }"
+      "    s = asin(sinh(r) * c / cos(alpha)) / c;"
+      "    gphi = 2.*sin(alpha)*s - atan(sin(alpha) * tan(c*s) / c);"
+      "    if(gphi > phi) { next_nan = 0; maxa = alpha; }"
+      "    else mina = alpha;"
+      "    }"
+      "  if(next_nan != 0) {"
+      "    mina = PI/4.;"
+      "    for(int it=0; it<10; it++) {"
+      "      alpha = (mina + maxa) / 2.;"
+      "      c = sqrt(-cos(2. * alpha));"
+      "      s = PI - asin(sinh(r) * c / cos(alpha)) / c;"
+      "      gphi = 2.*sin(alpha)*s - atan(sin(alpha) * tan(c*s) / c) - PI;"
+      "      if(gphi < phi) maxa = alpha; else mina = alpha;"
+      "      }"
+      "    beta = theta + atan(sin(alpha) * tanh(c * s) / c) + PI;"
+      "    }"
+      "  else beta = theta + atan(sin(alpha) * tanh(c * s) / c);"
+      "  }"
+      "alpha = alpha * sgn; beta = beta * sgn;"
+      "return vec4(s * cos(beta) * cos(alpha), s * sin(beta) * cos(alpha), s * sin(alpha), 1);"
+      "}";
+
+  EX transmatrix adjmatrix(int i, int j) {
+    ld zs = 2 * M_PI / 28;
+    if(i == 7) return zpush(-zs) * spin(-2*zs);
+    if(i == 8) return zpush(+zs) * spin(+2*zs);
+    return spin(2 * M_PI * i / 7) * xpush(tessf7/2) * spin(M_PI - 2 * M_PI * j / 7);
+    }
+  
+  struct hrmap_psl2 : hrmap {
+    
+    hrmap *underlying_map;
+    
+    map<pair<cell*, int>, cell*> at;
+    map<cell*, pair<cell*, int>> where;
+    
+    heptagon *getOrigin() override { return underlying_map->getOrigin(); }
+    
+    template<class T> auto in_underlying(const T& t) -> decltype(t()) {
+      dynamicval<eGeometry> g(geometry, gNormal);
+      dynamicval<hrmap*> gu(currentmap, underlying_map);
+      return t();
+      }
+    
+    cell *getCell(cell *u, int h) {
+      h = gmod(h, 14);
+      cell*& c = at[make_pair(u, h)];
+      if(!c) { c = newCell(u->type+2, u->master); where[c] = {u, h}; }
+      return c;
+      }
+  
+    cell* gamestart() override { return getCell(underlying_map->gamestart(), 0); }
+  
+    virtual transmatrix relative_matrix(cell *c2, cell *c1, const struct hyperpoint& point_hint) override { 
+      for(int i=0; i<c1->type; i++) if(c1->move(i) == c2) return adjmatrix(i, c1->c.spin(i));
+      if(gmatrix0.count(c2) && gmatrix0.count(c1))
+        return inverse(gmatrix0[c1]) * gmatrix0[c2];
+      return Id; // not implemented yet
+      }
+
+    void draw() override {
+      set<cell*> visited;
+      
+      cell* start = viewcenter();
+      vector<pair<cell*, transmatrix>> dq;
+      
+      visited.insert(start);
+      dq.emplace_back(start, cview());
+      
+      for(int i=0; i<isize(dq); i++) {
+        cell *c = dq[i].first;
+        transmatrix V = dq[i].second;
+        
+        if(V[3][3] < 0) V = centralsym * V;
+        if(!do_draw(c, V)) continue;
+        drawcell(c, V, 0, false);
+
+        for(int i=0; i<9; i++) {
+          // note: need do cmove before c.spin
+          cell *c1 = c->cmove(i);
+          if(visited.count(c1)) continue;
+          visited.insert(c1);
+          dq.emplace_back(c1, V * adjmatrix(i, c->c.spin(i)));
+          }
+        }
+      }
+
+    hrmap_psl2() {
+      in_underlying([this] { initcells(); underlying_map = currentmap; });
+      for(hrmap*& m: allmaps) if(m == underlying_map) m = NULL;
+      }
+    
+    ~hrmap_psl2() {
+      in_underlying([this] { delete currentmap; });
+      for(auto& p: at) tailored_delete(p.second);
+      }
+  
+    };
+
+  EX cell *get_at(cell *base, int level) {
+    return ((hrmap_psl2*)currentmap)->getCell(base, level);
+    }
+  
+  EX pair<cell*, int> get_where(cell *c) { return ((hrmap_psl2*)currentmap)->where[c]; }
+  
+  void find_cell_connection(cell *c, int d) {
+    auto m = (hrmap_psl2*)currentmap;
+    if(d >= c->type - 2) {
+      cell *c1 = get_at(m->where[c].first, m->where[c].second + (d == c->type-1 ? 1 : -1));
+      c->c.connect(d, c1, c1->type - 3 + c->type - d, false);
+      }
+    else {
+      auto cu = m->where[c].first;
+      auto cu1 = m->in_underlying([&] { return cu->cmove(d); });
+      int d1 = cu->c.spin(d);
+      cell *c1 = get_at(cu1, m->where[c].second - d1*2 + d*2 + 7);
+      c->c.connect(d, c1, d1, false);
+      }
+    }
+
+EX }
+
 EX namespace nisot {
 
   EX hyperpoint christoffel(const hyperpoint at, const hyperpoint velocity, const hyperpoint transported) {
     if(sol) return solv::christoffel(at, velocity, transported);
     else if(nil) return nilv::christoffel(at, velocity, transported);
+    else if(sl2) return slr::christoffel(at, velocity, transported);
     else return point3(0, 0, 0);
     }
 
@@ -789,6 +1139,7 @@ EX namespace nisot {
   EX hyperpoint inverse_exp(const hyperpoint h, iePrecision p, bool just_direction IS(true)) {
     if(sol) return solv::get_inverse_exp(h, p == iLazy, just_direction);
     if(nil) return nilv::get_inverse_exp(h, p == iLazy ? 5 : 20);
+    if(sl2) return slr::get_inverse_exp(h);
     if(prod) return product::inverse_exp(h);
     return point3(h[0], h[1], h[2]);
     }
@@ -822,6 +1173,7 @@ EX namespace nisot {
   EX hyperpoint get_exp(hyperpoint v, int steps) {
     if(sol) return direct_exp(v, steps);
     if(nil) return nilv::formula_exp(v);
+    if(sl2) return slr::formula_exp(v);
     if(prod) return product::direct_exp(v);
     return v;
     }
@@ -831,19 +1183,35 @@ EX namespace nisot {
     hyperpoint h = tC0(T);
     h[3] = 0;
   
-    h = Pos * h;
+    auto tPos = transpose(Pos);
+    
+    const ld eps = 1e-4;
+
+    if(sl2) {
+      hyperpoint p = slr::to_phigans(tPos[3]);
+      for(int i=0; i<3; i++)
+        tPos[i] = (slr::to_phigans(tPos[3] + tPos[i] * eps) - p) / eps;
+      tPos[3] = p; 
+      h = transpose(tPos) * h;
+      }
+    else h = Pos * h;
 
     int steps = 100;
     h /= steps;
     
-    auto tPos = transpose(Pos);
-  
     for(int i=0; i<steps; i++) {
       for(int j=0; j<3; j++)
         tPos[j] += christoffel(tPos[3], h, tPos[j]);
       geodesic_step(tPos[3], h);
       }
                                                                                   
+    if(sl2) {
+      hyperpoint p = slr::from_phigans(tPos[3]);
+      for(int i=0; i<3; i++)
+        tPos[i] = (slr::from_phigans(tPos[3] + tPos[i] * eps) - p) / eps;
+      tPos[3] = p;
+      }
+
     return transpose(tPos);
     }
 
@@ -895,6 +1263,7 @@ EX namespace nisot {
   EX hrmap *new_map() { 
     if(sol) return new solv::hrmap_sol; 
     if(nil) return new nilv::hrmap_nil;
+    if(sl2) return new slr::hrmap_psl2;
     if(geometry == gProduct) return new product::hrmap_product;
     return NULL;
     }

polygons.cpp

@@ -704,7 +704,7 @@ void geometry_information::make_wall(int id, vector<hyperpoint> vertices, vector
   set_column(T, 1, vertices[1]);
   set_column(T, 2, vertices[2]);
   set_column(T, 3, C0);
-  if(det(T) < 0)
+  if(det(T) < 0 && !sl2)
     reverse(vertices.begin(), vertices.end()),
     reverse(weights.begin(), weights.end());
 
@@ -716,6 +716,8 @@ void geometry_information::make_wall(int id, vector<hyperpoint> vertices, vector
   hyperpoint center = Hypc;
   int n = isize(vertices);
 
+  bool triangles = n > 3 && hypot_d(MDIM, vertices[0] - vertices[3]) < 1e-5;
+
   vector<ld> altitudes;
   if(prod) {
     altitudes.resize(n);
@@ -735,6 +737,7 @@ void geometry_information::make_wall(int id, vector<hyperpoint> vertices, vector
   if(prod) for(int i=0; i<n; i++) center_altitude += altitudes[i] * weights[i] / w;
   
   for(int a=0; a<n; a++) {
+    if(triangles) center = normalize(vertices[a/3*3] * 5 + vertices[a/3*3+1] + vertices[a/3*3+2]);
     hyperpoint v1 = vertices[a] - center;
     hyperpoint v2 = vertices[(a+1)%n] - center;
     texture_order([&] (ld x, ld y) {
@@ -755,6 +758,7 @@ void geometry_information::make_wall(int id, vector<hyperpoint> vertices, vector
   if(true) {
     int STEP = vid.texture_step;
     for(int a=0; a<n; a++) for(int y=0; y<STEP; y++) {
+      if(triangles && (a%3 != 1)) continue;
       hyperpoint h = (vertices[a] * (STEP-y) + vertices[(a+1)%n] * y)/STEP;
       if(prod) { 
         h = zshift(normalize_flat(h), (altitudes[a] * (STEP-y) + altitudes[(a+1)%n] * y) / STEP);
@@ -799,7 +803,7 @@ void geometry_information::reserve_wall3d(int i) {
 
 void geometry_information::create_wall3d() {
   if(WDIM == 2) return;
-  reserve_wall3d(penrose ? 22 : prod ? 0 : S7);
+  reserve_wall3d(penrose ? 22 : prod ? 0 : sl2 ? 9 : S7);
   if(GDIM == 3 && binarytiling && geometry == gBinary3) {
     hyperpoint h00 = point3(-1,-1,-1);
     hyperpoint h01 = point3(-1,0,-1);
@@ -954,7 +958,7 @@ void geometry_information::create_wall3d() {
       }
     }
 
-  if(GDIM == 3 && !euclid && !binarytiling && !nil && !prod) {
+  if(GDIM == 3 && !euclid && !binarytiling && !nil && !prod && !sl2) {
     reg3::generate();
     int facesize = isize(reg3::cellshape) / S7;
     for(int w=0; w<S7; w++) {
@@ -965,6 +969,42 @@ void geometry_information::create_wall3d() {
       }
     }
 
+  if(geometry == gSL2) {
+    ld zs = 2 * M_PI / 28;
+    ld a = 2 * M_PI/7;
+    ld tf = tessf7 / 4;
+    ld halfedge = 0.283128;
+    ld he = halfedge / 2;
+    hyperpoint right_u = xpush(tf) * ypush(-he) * zpush0(zs/2);
+    hyperpoint right_d = xpush(tf) * ypush(-he) * zpush0(-zs/2);
+    hyperpoint left_u  = xpush(tf) * ypush(+he) * zpush0(zs/2);
+    hyperpoint left_d  = xpush(tf) * ypush(+he) * zpush0(-zs/2);
+    hyperpoint center_u = zpush0(zs/2);
+    hyperpoint center_d = zpush0(-zs/2);
+    for(int i=0; i<7; i++) {
+      auto s =spin(a * i);
+      make_wall(i, {s * right_u, s * right_d, s * left_d, s * left_u});
+      }
+    vector<hyperpoint> top, bot;
+    for(int i=0; i<7; i++) {
+      bot.push_back(center_d);
+      bot.push_back(spin(a*i) * left_d);
+      bot.push_back(spin(a*i) * right_d);
+      bot.push_back(center_d);
+      bot.push_back(spin(a*i) * right_d);
+      bot.push_back(spin(a*(i+1)) * left_d);
+
+      top.push_back(center_u);
+      top.push_back(spin(a*i) * left_u);
+      top.push_back(spin(a*i) * right_u);
+      top.push_back(center_u);
+      top.push_back(spin(a*i) * right_u);
+      top.push_back(spin(a*(i+1)) * left_u);
+      }
+    make_wall(7, bot);
+    make_wall(8, top);
+    }
+  
   if(geometry == gSol) {
     ld zstep = -log(2) / 2;
     ld bwh = vid.binary_width * zstep;

rug.cpp

@@ -997,7 +997,7 @@ bincode acd_bin(ld x) {
 
 bincode get_bincode(hyperpoint h) {
   switch(ginf[gwhere].cclass) {
-    case gcEuclid: case gcSol: case gcNil: case gcProduct:
+    case gcEuclid: case gcSol: case gcNil: case gcProduct: case gcSL2:
       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

@@ -53,7 +53,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,
+  standardEH2, standardSL2,
   MAX 
   };
 
@@ -635,6 +635,7 @@ void init() {
     bool sh3 = (sp == shader_projection::standardH3);
     bool ssol = (sp == shader_projection::standardSolv);
     bool snil = (sp == shader_projection::standardNil);
+    bool ssl2 = (sp == shader_projection::standardSL2);
     bool sr3 = (sp == shader_projection::standardR3);
     bool ss30 = (sp == shader_projection::standardS30);
     bool ss31 = (sp == shader_projection::standardS31);
@@ -643,7 +644,7 @@ void init() {
     bool seh2 = (sp == shader_projection::standardEH2);
     bool ss3 = ss30 || ss31 || ss32 || ss33;
     
-    bool s3 = (sh3 || sr3 || ss3 || ssol || snil || seh2);
+    bool s3 = (sh3 || sr3 || ss3 || ssol || snil || seh2 || ssl2);
     bool dim3 = s3 || among(sp, shader_projection::ball, shader_projection::halfplane3, shader_projection::band3);
     bool dim2 = !dim3;
     bool ball = (sp == shader_projection::ball);
@@ -697,6 +698,7 @@ void init() {
       
       ssol,    solv::solshader,      
       snil,    nilv::nilshader,
+      ssl2,    slr::slshader,
 
       1,       "void main() {",  
       texture,   "vTexCoord = aTexture;",
@@ -736,6 +738,7 @@ void init() {
       ssol,      "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; ",
       snil,      "t = inverse_exp(t);",
+      ssl2,      "t = inverse_exp(t);",
 
       seh2,      "float z = log(t[2] * t[2] - t[0] * t[0] - t[1] * t[1]) / 2.;",
       seh2,      "float r = sqrt(t[0] * t[0] + t[1] * t[1]);",
@@ -747,7 +750,7 @@ void init() {
       seh2,      "t[2] = z;",
 
       sh3,       "float fogs = (uFogBase - acosh(t[3]) / uFog);",
-      sr3||snil, "float fogs = (uFogBase - sqrt(t[0]*t[0] + t[1]*t[1] + t[2]*t[2]) / 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);",
 
       seh2,      "float fogs = (uFogBase - sqrt(z*z+d*d) / uFog);",

system.cpp

@@ -1187,7 +1187,7 @@ EX void set_geometry(eGeometry target) {
     if(GDIM == 3 && old_DIM == 2 && pmodel == mdDisk) pmodel = mdPerspective;
     if(nonisotropic && old_DIM == 2) pmodel = mdGeodesic;
     if(GDIM == 2 && among(pmodel, mdPerspective, mdGeodesic)) pmodel = mdDisk;
-    if(nonisotropic && old_DIM == 2 && vid.texture_step < 4) vid.texture_step = 4;
+    if(nonisotropic && old_DIM == 2 && vid.texture_step < 4 && !sl2) vid.texture_step = 4;
     if(prod) { pmodel = mdPerspective; if(vid.texture_step < 4) vid.texture_step = 4; }
     if(prod && shmup::on) shmup::on = false;
     }
@@ -1196,7 +1196,7 @@ EX void set_geometry(eGeometry target) {
 EX void set_variation(eVariation target) {
   if(variation != target) {
     stop_game();
-    if(euclid6 || binarytiling || sol || penrose) geometry = gNormal;
+    if(euclid6 || binarytiling || sol || penrose || WDIM == 3) if(!prod) geometry = gNormal;
     auto& cd = ginf[gCrystal];
     if(target == eVariation::bitruncated && cryst && cd.sides == 8 && cd.vertex == 4) {
       cd.vertex = 3;