From 6958cbcbd91d7e0956fb43b8141e024431006fad Mon Sep 17 00:00:00 2001
From: Zeno Rogue <zeno@attnam.com>
Date: Sat, 17 Aug 2019 23:28:41 +0200
Subject: [PATCH] product:: preliminary version (no turning)

---
 archimedean.cpp   |   2 +-
 basegraph.cpp     |   2 +
 binary-tiling.cpp |   6 +-
 cell.cpp          |   9 ++-
 classes.cpp       |   2 +-
 classes.h         |   4 +-
 complex.cpp       |   2 +-
 control.cpp       |   6 +-
 drawing.cpp       |  17 +++++-
 euclid.cpp        |   4 +-
 floorshapes.cpp   |   2 +-
 geom-exp.cpp      |   4 +-
 geometry2.cpp     |  10 ++--
 graph.cpp         |  26 +++++----
 hyper.h           |   6 +-
 hyperpoint.cpp    | 105 +++++++++++++++++++++++-----------
 hypgraph.cpp      |  84 +++++++++++++++++-----------
 irregular.cpp     |  10 ++--
 nonisotropic.cpp  | 139 +++++++++++++++++++++++++++++++++++++++++++++-
 polygons.cpp      |  18 +++++-
 rug.cpp           |   4 +-
 sysconfig.h       |   3 -
 22 files changed, 349 insertions(+), 116 deletions(-)

diff --git a/archimedean.cpp b/archimedean.cpp
index 4f4d2aac..699163d4 100644
--- a/archimedean.cpp
+++ b/archimedean.cpp
@@ -555,7 +555,7 @@ struct hrmap_archimedean : hrmap {
       }
     
     if(euclid) 
-      alt = encodeId(pair_to_vec(int(T[0][GDIM]), int(T[1][GDIM])));
+      alt = encodeId(pair_to_vec(int(T[0][LDIM]), int(T[1][LDIM])));
       
     DEBB(DF_GEOM, ("look for: ", alt, " / ", T * C0));
   
diff --git a/basegraph.cpp b/basegraph.cpp
index 56e6f189..941bd7c1 100644
--- a/basegraph.cpp
+++ b/basegraph.cpp
@@ -297,6 +297,8 @@ void display_data::set_projection(int ed) {
       shaderside_projection = true, glhr::new_shader_projection = glhr::shader_projection::standardH3, pers3 = true;
     if(GDIM == 3 && translatable && apply_models && pmodel == mdPerspective)
       shaderside_projection = true, glhr::new_shader_projection = glhr::shader_projection::standardR3, pers3 = true;
+    if(GDIM == 3 && prod && apply_models && pmodel == mdPerspective)
+      shaderside_projection = true, glhr::new_shader_projection = glhr::shader_projection::standardR3, 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 && nil)
diff --git a/binary-tiling.cpp b/binary-tiling.cpp
index ab249c79..b79a321e 100644
--- a/binary-tiling.cpp
+++ b/binary-tiling.cpp
@@ -643,7 +643,7 @@ EX namespace binary {
 
   // on which horocycle are we
   EX ld horo_level(hyperpoint h) {
-    h /= (1 + h[GDIM]);
+    h /= (1 + h[LDIM]);
     h[0] -= 1;
     h /= sqhypot_d(GDIM, h);
     h[0] += .5;
@@ -848,7 +848,7 @@ EX void virtualRebaseSimple(heptagon*& base, transmatrix& at) {
 
   while(true) {
   
-    double currz = at[GDIM][GDIM];
+    double currz = at[LDIM][LDIM];
     
     heptagon *h = base;
     
@@ -858,7 +858,7 @@ EX void virtualRebaseSimple(heptagon*& base, transmatrix& at) {
     
     for(int d=0; d<S7; d++) {
       transmatrix V2 = itmatrix(h, d) * at;
-      double newz = V2[GDIM][GDIM];
+      double newz = V2[LDIM][LDIM];
       if(newz < currz) {
         currz = newz;
         bestV = V2;
diff --git a/cell.cpp b/cell.cpp
index b0852ccf..5c0be1a4 100644
--- a/cell.cpp
+++ b/cell.cpp
@@ -146,8 +146,8 @@ EX cell *createMov(cell *c, int d) {
     }
   
   if(c->move(d)) return c->move(d);
-  PROD( else if(geometry == gProduct)
-    product::find_cell_connection(c, d); )
+  else if(geometry == gProduct)
+    product::find_cell_connection(c, d);
   #if CAP_BT
   else if(penrose)
     kite::find_cell_connection(c, d);
@@ -263,11 +263,10 @@ EX void initcells() {
   
   hrmap* res = callhandlers((hrmap*)nullptr, hooks_newmap);
   if(res) currentmap = res;  
-  else if(nonisotropic) currentmap = nisot::new_map();
+  else if(nonisotropic || prod) currentmap = nisot::new_map();
   #if CAP_CRYSTAL
   else if(geometry == gCrystal) currentmap = crystal::new_map();
   #endif
-  PROD( else if(geometry == gProduct) currentmap = product::new_map(); )
   #if CAP_ARCM
   else if(archimedean) currentmap = arcm::new_map();
   #endif
@@ -862,7 +861,7 @@ int ld_to_int(ld x) {
 
 EX int pseudocoords(cell *c) {
   transmatrix T = arcm::archimedean_gmatrix[c->master].second;
-  return pair_to_vec(ld_to_int(T[0][GDIM]), ld_to_int((spin(60*degree) * T)[0][GDIM]));
+  return pair_to_vec(ld_to_int(T[0][LDIM]), ld_to_int((spin(60*degree) * T)[0][LDIM]));
   }
 
 EX cdata *arcmCdata(cell *c) {
diff --git a/classes.cpp b/classes.cpp
index d6ab8446..b86335cc 100644
--- a/classes.cpp
+++ b/classes.cpp
@@ -571,7 +571,7 @@ vector<geometryinfo> ginf = {
   {"kd2",    "none",    "kite-and-dart",                              "kd2",      4, 3, qPENROSE,  gcEuclid,     0x48000, {{7, 7}}, eVariation::pure},
   {"kd3",    "none",    "kite-and-dart on horospheres",               "kd3",     12, 3, qsBP,      gcHyperbolic, 0x48200, {{7, 3}}, eVariation::pure},
   {"nil",    "none",    "Nil geometry",                               "nil",      8, 3, 0,         gcNil,        0x48600, {{7, 5}}, eVariation::pure},
-//{"product","none",    "product space",                              "product",  7, 3, 0,         gcProduct,    0x48400, {{7, 3}}, eVariation::pure},
+  {"product","none",    "product space",                              "product",  7, 3, 0,         gcProduct,    0x48400, {{7, 3}}, eVariation::pure},
   };
 
   // bits: 9, 10, 15, 16, (reserved for later) 17, 18
diff --git a/classes.h b/classes.h
index f209506b..7283d6f1 100644
--- a/classes.h
+++ b/classes.h
@@ -214,10 +214,10 @@ enum eGeometry {
   gHoroTris, gHoroRec, gHoroHex,
   gField435, gField534,
   gBinary4, gSol,
-  gKiteDart2, gKiteDart3, gNil, PROD2(gProduct,)
+  gKiteDart2, gKiteDart3, gNil, gProduct,
   gGUARD};
 
-enum eGeometryClass { gcHyperbolic, gcEuclid, gcSphere, gcSol, gcNil, PROD(gcProduct) };
+enum eGeometryClass { gcHyperbolic, gcEuclid, gcSphere, gcSol, gcNil, gcProduct };
 
 enum class eVariation { bitruncated, pure, goldberg, irregular, dual };  
 
diff --git a/complex.cpp b/complex.cpp
index 274e62e4..6743e41e 100644
--- a/complex.cpp
+++ b/complex.cpp
@@ -2742,7 +2742,7 @@ EX namespace sword {
     int s1 = neighborId(c1, c2);
     int s2 = neighborId(c2, c1);
     if(s1 < 0 || s2 < 0) return d;
-    if(WDIM == 2) {
+    if(WDIM == 2 || prod) {
       if(c1->c.mirror(s1))
         d.angle = ((s2*sword_angles/c2->type - d.angle + s1*sword_angles/c1->type) + sword_angles/2) % sword_angles;
       else
diff --git a/control.cpp b/control.cpp
index 4b4df757..ef9bb882 100644
--- a/control.cpp
+++ b/control.cpp
@@ -86,7 +86,7 @@ EX movedir vectodir(const hyperpoint& P) {
   for(int i=0; i<cwt.at->type; i++) {
     transmatrix T;
     if(compute_relamatrix((cwt+i).peek(), cwt.at, i, T)) {
-      dirdist[i] = intval(U * T * C0, Centered * P);
+      dirdist[i] = quickdist(U * T * C0, Centered * P);
       }
     //xspinpush0(-i * 2 * M_PI /cwt.at->type, .5), P);
     }
@@ -115,6 +115,8 @@ EX hyperpoint move_destination_vec(int d) {
   // else if(WDIM == 2 && pmodel == mdPerspective) return cspin(0, 2, d * M_PI/4) * tC0(pushone());
   // else if(WDIM == 2) return spin(-d * M_PI/4) * tC0(pushone());
   else if(d&1) return cspin(0, 1, d > 4 ? M_PI/2 : -M_PI/2) * tC0(pushone());
+  else if(prod && d == 6) return zshift(C0, product::plevel);
+  else if(prod && d == 2) return zshift(C0, -product::plevel);
   else return cspin(0, 2, d * M_PI/4) * tC0(pushone());
   }
 
@@ -907,7 +909,7 @@ EX void handle_event(SDL_Event& ev) {
           vid.xposition += (mousex - lmousex) * 1. / current_display->scrsize,
           vid.yposition += (mousey - lmousey) * 1. / current_display->scrsize;
           }
-        else if(mouseh[GDIM] < 50 && mouseoh[GDIM] < 50) {
+        else if(mouseh[LDIM] < 50 && mouseoh[LDIM] < 50) {
           panning(mouseoh, mouseh);
           }
         }
diff --git a/drawing.cpp b/drawing.cpp
index ec6d96b4..4dde8b1d 100644
--- a/drawing.cpp
+++ b/drawing.cpp
@@ -656,7 +656,7 @@ EX ld scale_at(const transmatrix& T) {
 EX ld linewidthat(const hyperpoint& h) {
   if(!(vid.antialias & AA_LINEWIDTH)) return 1;
   else if(hyperbolic && pmodel == mdDisk && vid.alpha == 1 && !ISWEB) {
-    double dz = h[GDIM];
+    double dz = h[LDIM];
     if(dz < 1 || abs(dz-current_display->scrdist) < 1e-6) return 1;
     else {
       double dx = sqrt(dz * dz - 1);
@@ -893,6 +893,21 @@ void debug_this() { }
 void dqi_poly::draw() {
   if(flags & POLY_DEBUG) debug_this();
 
+  if(prod && vid.usingGL && pmodel == mdPerspective && (current_display->set_all(global_projection), shaderside_projection)) {
+    auto npoly = *this;
+    glcoords.clear();
+    for(int i=0; i<cnt; i++)
+      glcoords.push_back(glhr::pointtogl(product::where(V * glhr::gltopoint( (*tab)[offset+i]))));
+      
+    npoly.offset = 0;
+    npoly.tab = &glcoords;
+    npoly.V = Id;
+    set_width(1);
+    npoly.gldraw();
+    return;
+    }
+
+
   dynamicval<ld> bs(hr::band_shift, band_shift);
   if(!hyperbolic && among(pmodel, mdPolygonal, mdPolynomial)) {
     bool any = false;
diff --git a/euclid.cpp b/euclid.cpp
index 97767404..de939d2e 100644
--- a/euclid.cpp
+++ b/euclid.cpp
@@ -617,7 +617,7 @@ EX namespace euclid3 {
       tmatrix.resize(S7);
       for(int i=0; i<S7; i++) tmatrix[i] = Id;
       for(int i=0; i<S7; i++) for(int j=0; j<3; j++)
-        tmatrix[i][j][GDIM] = getcoord(shifttable[i])[j];
+        tmatrix[i][j][LDIM] = getcoord(shifttable[i])[j];
       camelot_center = NULL;
       build_torus3();
       }
@@ -1195,7 +1195,7 @@ EX namespace euclid3 {
 #endif
 
 EX ld matrixnorm(const transmatrix& Mat) {
-  return Mat[0][GDIM] * Mat[0][GDIM] + Mat[1][GDIM] * Mat[1][GDIM] + Mat[2][GDIM] * Mat[2][GDIM];
+  return Mat[0][LDIM] * Mat[0][LDIM] + Mat[1][LDIM] * Mat[1][LDIM] + Mat[2][LDIM] * Mat[2][LDIM];
   }
   
 void hrmap_euclid_any::draw() {
diff --git a/floorshapes.cpp b/floorshapes.cpp
index 7494dced..5bb5be96 100644
--- a/floorshapes.cpp
+++ b/floorshapes.cpp
@@ -224,7 +224,7 @@ void geometry_information::bshape2(hpcshape& sh, PPR prio, int shapeid, matrixli
       int mapped = 0;
       for(auto& m: matrices) {        
         hyperpoint z = m.first * h;
-        if(z[0] > -1e-5 && z[1] > -1e-5 && z[GDIM] > -1e-5) {
+        if(z[0] > -1e-5 && z[1] > -1e-5 && z[LDIM] > -1e-5) {
           nh = m.second[r] * z, mapped++;
           }
         }
diff --git a/geom-exp.cpp b/geom-exp.cpp
index b531596a..42d501ed 100644
--- a/geom-exp.cpp
+++ b/geom-exp.cpp
@@ -811,9 +811,9 @@ EX void showEuclideanMenu() {
       dialog::addSelItem(XLAT("Curvature"), XLAT("Nil"), 0);
       break;
 
-    PROD( case gcProduct:
+    case gcProduct:
       dialog::addSelItem(XLAT("Curvature"), XLAT("Product"), 0);
-      break; )
+      break;
     }    
     
   dialog::display();
diff --git a/geometry2.cpp b/geometry2.cpp
index 912eb09a..bca48217 100644
--- a/geometry2.cpp
+++ b/geometry2.cpp
@@ -10,14 +10,14 @@ namespace hr {
 transmatrix &ggmatrix(cell *c);
 
 void fixelliptic(transmatrix& at) {
-  if(elliptic && at[GDIM][GDIM] < 0) {
+  if(elliptic && at[LDIM][LDIM] < 0) {
     for(int i=0; i<MDIM; i++) for(int j=0; j<MDIM; j++)
       at[i][j] = -at[i][j];
     }
   }
 
 void fixelliptic(hyperpoint& h) {
-  if(elliptic && h[GDIM] < 0)
+  if(elliptic && h[LDIM] < 0)
     for(int i=0; i<MDIM; i++) h[i] = -h[i];
   }
 
@@ -155,7 +155,7 @@ transmatrix hrmap_standard::relative_matrix(cell *c2, cell *c1, const hyperpoint
 
 EX transmatrix &ggmatrix(cell *c) {
   transmatrix& t = gmatrix[c];
-  if(t[GDIM][GDIM] == 0) {
+  if(t[LDIM][LDIM] == 0) {
     if(euwrap && centerover.at && masterless) 
       t = calc_relative_matrix(c, centerover.at, C0);
     else if(masterless && WDIM == 2) {
@@ -352,7 +352,7 @@ EX void virtualRebaseSimple(heptagon*& base, transmatrix& at) {
 
   while(true) {
   
-    double currz = at[GDIM][GDIM];
+    double currz = at[LDIM][LDIM];
     
     heptagon *h = base;
     
@@ -364,7 +364,7 @@ EX void virtualRebaseSimple(heptagon*& base, transmatrix& at) {
       heptspin hs(h, d, false);
       heptspin hs2 = hs + wstep;
       transmatrix V2 = spin(-hs2.spin*2*M_PI/S7) * cgi.invheptmove[d] * at;
-      double newz = V2[GDIM][GDIM];
+      double newz = V2[LDIM][LDIM];
       if(newz < currz) {
         currz = newz;
         bestV = V2;
diff --git a/graph.cpp b/graph.cpp
index 4d09bccf..be409197 100644
--- a/graph.cpp
+++ b/graph.cpp
@@ -349,12 +349,14 @@ double hexshiftat(cell *c) {
   }
 
 EX transmatrix ddspin(cell *c, int d, ld bonus IS(0)) {
+  if(prod) return PIU( ddspin(c, d, bonus) );
   if(WDIM == 3 && d < c->type) return rspintox(tC0(calc_relative_matrix(c->cmove(d), c, C0))) * cspin(2, 0, bonus);
   if(WDIM == 2 && (binarytiling || penrose) && d < c->type) return spin(bonus) * rspintox(nearcorner(c, d));
   return spin(displayspin(c, d) + bonus - hexshiftat(c));
   }
 
 EX transmatrix iddspin(cell *c, int d, ld bonus IS(0)) {
+  if(prod) return PIU( iddspin(c, d, bonus) );
   if(WDIM == 3 && d < c->type) return cspin(0, 2, bonus) * spintox(tC0(calc_relative_matrix(c->cmove(d), c, C0)));
   if(WDIM == 2 && (binarytiling || penrose) && d < c->type) return spin(bonus) * spintox(nearcorner(c, d));
   return spin(hexshiftat(c) - displayspin(c, d) + bonus);
@@ -5018,7 +5020,7 @@ void drawcell_in_radar(cell *c, transmatrix V) {
 
 EX void drawcell(cell *c, transmatrix V, int spinv, bool mirrored) {
  
-  PROD( if(product::pmap) { product::drawcell_stack(c, V, spinv, mirrored); return; } )
+  if(product::pmap) { product::drawcell_stack(c, V, spinv, mirrored); return; }
 
   cells_drawn++;
 
@@ -5037,10 +5039,10 @@ EX void drawcell(cell *c, transmatrix V, int spinv, bool mirrored) {
   if(!inmirrorcount) {
     transmatrix& gm = gmatrix[c];
     orig = 
-      gm[GDIM][GDIM] == 0 ? true : 
+      gm[LDIM][LDIM] == 0 ? true : 
       euwrap ? hdist0(tC0(gm)) >= hdist0(tC0(V)) :
-      sphereflipped() ? fabs(gm[GDIM][GDIM]-1) <= fabs(V[GDIM][GDIM]-1) :
-      fabs(gm[GDIM][GDIM]-1) >= fabs(V[GDIM][GDIM]-1) - 1e-8;
+      sphereflipped() ? fabs(gm[LDIM][LDIM]-1) <= fabs(V[LDIM][LDIM]-1) :
+      fabs(gm[LDIM][LDIM]-1) >= fabs(V[LDIM][LDIM]-1) - 1e-8;
 
     if(orig) gm = V;
     }
@@ -5142,7 +5144,7 @@ EX void drawcell(cell *c, transmatrix V, int spinv, bool mirrored) {
       }
 
     if(!masterless) {
-      double dfc = euclid ? intval(tC0(V), C0) : V[GDIM][GDIM];
+      double dfc = euclid ? intval(tC0(V), C0) : V[LDIM][LDIM];
     
       if(dfc < centdist) {
         centdist = dfc;
@@ -6001,12 +6003,12 @@ EX void drawcell(cell *c, transmatrix V, int spinv, bool mirrored) {
                 ld b = vid.binary_width * log(2) / 2;
                 const ld l = log(2) / 2;
                 switch(a) {
-                  case 0: if(V[0][GDIM] >= b) continue; break;
-                  case 1: if(V[1][GDIM] >= b) continue; break;
-                  case 2: case 3: if (pmodel == mdPerspective && V[2][GDIM] >= l) continue; break;
-                  case 4: if(V[0][GDIM] <= -b) continue; break;
-                  case 5: if(V[1][GDIM] <= -b) continue; break;
-                  case 6: case 7: if (pmodel == mdPerspective && V[2][GDIM] <= -l) continue; break;
+                  case 0: if(V[0][LDIM] >= b) continue; break;
+                  case 1: if(V[1][LDIM] >= b) continue; break;
+                  case 2: case 3: if (pmodel == mdPerspective && V[2][LDIM] >= l) continue; break;
+                  case 4: if(V[0][LDIM] <= -b) continue; break;
+                  case 5: if(V[1][LDIM] <= -b) continue; break;
+                  case 6: case 7: if (pmodel == mdPerspective && V[2][LDIM] <= -l) continue; break;
                   }
                 }
               if(qfi.fshape && wmescher) {
@@ -7140,7 +7142,7 @@ EX ld wall_radar(cell *c, transmatrix T, ld max) {
 
 EX void make_actual_view() {
   sphereflip = Id;
-  if(sphereflipped()) sphereflip[GDIM][GDIM] = -1;
+  if(sphereflipped()) sphereflip[LDIM][LDIM] = -1;
   actual_view_transform = sphereflip;  
   if(vid.yshift && WDIM == 2) actual_view_transform = ypush(vid.yshift) * actual_view_transform;
   #if MAXMDIM >= 4
diff --git a/hyper.h b/hyper.h
index 274be69f..1a08dd6c 100644
--- a/hyper.h
+++ b/hyper.h
@@ -124,6 +124,7 @@ void addMessage(string s, char spamtype = 0);
 #define sphere (cgclass == gcSphere)
 #define sol (cgclass == gcSol)
 #define nil (cgclass == gcNil)
+#define prod (cgclass == gcProduct)
 #define hyperbolic (cgclass == gcHyperbolic)
 #define nonisotropic (sol || nil)
 #define translatable (euclid || nonisotropic)
@@ -325,10 +326,11 @@ extern videopar vid;
 #if MAXMDIM == 3
 #define WDIM 2
 #else
-#define WDIM ((geometry >= gBinary3 && geometry != gBinary4 && geometry != gKiteDart2 PROD(&& geometry != gProduct)) ? 3 : 2)
+#define WDIM ((geometry >= gBinary3 && geometry != gBinary4 && geometry != gKiteDart2) ? 3 : 2)
 #endif
 #define GDIM (vid.always3 ? 3 : WDIM)
-#define MDIM (GDIM+1)
+#define MDIM (prod ? 3 : GDIM+1)
+#define LDIM (MDIM-1)
 
 #define self (*this)
 
diff --git a/hyperpoint.cpp b/hyperpoint.cpp
index 72c57a44..084336a8 100644
--- a/hyperpoint.cpp
+++ b/hyperpoint.cpp
@@ -140,7 +140,7 @@ inline hyperpoint point2(ld x, ld y) { return hyperpoint(x,y,0,0); }
 
 extern const hyperpoint C02, C03;
 
-#define C0 (GDIM == 2 ? C02 : C03)
+#define C0 (MDIM == 3 ? C02 : C03)
 #endif
 
 // basic functions and types
@@ -171,13 +171,14 @@ ld inverse_tanh(ld x) { return log((1+x)/(1-x)) / 2; } */
 
 EX ld squar(ld x) { return x*x; }
 
-EX int sig(int z) { return (sphere || sol || z<GDIM)?1:-1; }
+EX int sig(int z) { return prod ? (z<2?1:-1) : (sphere || sol || z<GDIM)?1:-1; }
 
 EX int curvature() {
   switch(cgclass) {
     case gcEuclid: return 0;
     case gcHyperbolic: return -1;
     case gcSphere: return 1;
+    case gcProduct: return PIU(curvature());
     default: return 0;
     }
   }
@@ -187,6 +188,7 @@ EX ld sin_auto(ld x) {
     case gcEuclid: return x;
     case gcHyperbolic: return sinh(x);
     case gcSphere: return sin(x);
+    case gcProduct: return PIU(sin_auto(x));
     default: return x;
     }
   }
@@ -196,6 +198,7 @@ EX ld asin_auto(ld x) {
     case gcEuclid: return x;
     case gcHyperbolic: return asinh(x);
     case gcSphere: return asin(x);
+    case gcProduct: return PIU(asin_auto(x));
     default: return x;
     }
   }
@@ -204,6 +207,7 @@ EX ld acos_auto(ld x) {
   switch(cgclass) {
     case gcHyperbolic: return acosh(x);
     case gcSphere: return acos(x);
+    case gcProduct: return PIU(acos_auto(x));
     default: return x;
     }
   }
@@ -231,6 +235,7 @@ EX ld cos_auto(ld x) {
     case gcEuclid: return 1;
     case gcHyperbolic: return cosh(x);
     case gcSphere: return cos(x);
+    case gcProduct: return PIU(cos_auto(x));
     default: return 1;
     }
   }
@@ -240,6 +245,7 @@ EX ld tan_auto(ld x) {
     case gcEuclid: return x;
     case gcHyperbolic: return tanh(x);
     case gcSphere: return tan(x);
+    case gcProduct: return PIU(tan_auto(x));
     default: return 1;
     }
   }
@@ -249,6 +255,7 @@ EX ld atan_auto(ld x) {
     case gcEuclid: return x;
     case gcHyperbolic: return atanh(x);
     case gcSphere: return atan(x);
+    case gcProduct: return PIU(atan_auto(x));
     default: return x;
     }
   }
@@ -258,6 +265,7 @@ EX ld atan2_auto(ld y, ld x) {
     case gcEuclid: return y/x;
     case gcHyperbolic: return atanh(y/x);
     case gcSphere: return atan2(y, x);
+    case gcProduct: return PIU(atan2_auto(y, x));
     default: return y/x;
     }
   }
@@ -313,6 +321,11 @@ EX ld intval(const hyperpoint &h1, const hyperpoint &h2) {
   return res;
   }
 
+EX ld quickdist(const hyperpoint &h1, const hyperpoint &h2) {
+  if(prod) return hdist(h1, h2);
+  return intval(h1, h2);
+  }
+
 EX ld sqhypot_d(int d, const hyperpoint& h) {
   ld sum = 0;
   for(int i=0; i<d; i++) sum += h[i]*h[i];
@@ -324,9 +337,10 @@ EX ld hypot_d(int d, const hyperpoint& h) {
   }
   
 EX ld zlevel(const hyperpoint &h) {
-  if(translatable) return h[GDIM];
+  if(translatable) return h[LDIM];
   else if(sphere) return sqrt(intval(h, Hypc));
-  else return (h[GDIM] < 0 ? -1 : 1) * sqrt(-intval(h, Hypc));
+  else if(prod) return log(sqrt(-intval(h, Hypc)));
+  else return (h[LDIM] < 0 ? -1 : 1) * sqrt(-intval(h, Hypc));
   }
 
 EX ld hypot_auto(ld x, ld y) {
@@ -344,6 +358,7 @@ EX ld hypot_auto(ld x, ld y) {
 
 // move H back to the sphere/hyperboloid/plane
 EX hyperpoint normalize(hyperpoint H) {
+  if(prod) return H;
   ld Z = zlevel(H);
   for(int c=0; c<MDIM; c++) H[c] /= Z;
   return H;
@@ -351,11 +366,17 @@ EX hyperpoint normalize(hyperpoint H) {
 
 // get the center of the line segment from H1 to H2
 hyperpoint mid(const hyperpoint& H1, const hyperpoint& H2) {
+  if(prod) {
+    auto d1 = product_decompose(H1);
+    auto d2 = product_decompose(H2);
+    return zshift(PIU( mid(d1.second, d2.second) ), (d1.first + d2.first) / 2);
+    }
   return normalize(H1 + H2);
   }
 
 // like mid, but take 3D into account
 EX hyperpoint midz(const hyperpoint& H1, const hyperpoint& H2) {
+  if(prod) return mid(H1, H2);
   hyperpoint H3 = H1 + H2;
   
   ld Z = 2;
@@ -394,15 +415,15 @@ EX transmatrix random_spin() {
 
 EX transmatrix eupush(ld x, ld y) {
   transmatrix T = Id;
-  T[0][GDIM] = x;
-  T[1][GDIM] = y;
+  T[0][LDIM] = x;
+  T[1][LDIM] = y;
   return T;
   }
 
 EX transmatrix eupush(hyperpoint h) {
   if(nonisotropic) return nisot::translate(h);
   transmatrix T = Id;
-  for(int i=0; i<GDIM; i++) T[i][GDIM] = h[i];
+  for(int i=0; i<GDIM; i++) T[i][LDIM] = h[i];
   return T;
   }
 
@@ -430,9 +451,9 @@ EX transmatrix cpush(int cid, ld alpha) {
   transmatrix T = Id;
   if(nonisotropic) 
     return eupush3(cid == 0 ? alpha : 0, cid == 1 ? alpha : 0, cid == 2 ? alpha : 0);
-  T[GDIM][GDIM] = T[cid][cid] = cos_auto(alpha);
-  T[cid][GDIM] = sin_auto(alpha);
-  T[GDIM][cid] = -curvature() * sin_auto(alpha);
+  T[LDIM][LDIM] = T[cid][cid] = cos_auto(alpha);
+  T[cid][LDIM] = sin_auto(alpha);
+  T[LDIM][cid] = -curvature() * sin_auto(alpha);
   return T;
   }
 
@@ -550,7 +571,7 @@ EX transmatrix rspintoc(const hyperpoint& H, int t, int f) {
 
 // rotate the hyperbolic plane around C0 such that H[1] == 0 and H[0] >= 0
 EX transmatrix spintox(const hyperpoint& H) {
-  if(GDIM == 2) return spintoc(H, 0, 1); 
+  if(GDIM == 2 || prod) return spintoc(H, 0, 1); 
   transmatrix T1 = spintoc(H, 0, 1);
   return spintoc(T1*H, 0, 2) * T1;
   }
@@ -572,7 +593,7 @@ EX transmatrix build_matrix(hyperpoint h1, hyperpoint h2, hyperpoint h3, hyperpo
 
 // reverse of spintox(H)
 EX transmatrix rspintox(const hyperpoint& H) {
-  if(GDIM == 2) return rspintoc(H, 0, 1); 
+  if(GDIM == 2 || prod) return rspintoc(H, 0, 1); 
   transmatrix T1 = spintoc(H, 0, 1);
   return rspintoc(H, 0, 1) * rspintoc(T1*H, 0, 2);
   }
@@ -580,16 +601,16 @@ EX transmatrix rspintox(const hyperpoint& H) {
 // for H such that H[1] == 0, this matrix pushes H to C0
 EX transmatrix pushxto0(const hyperpoint& H) {
   transmatrix T = Id;
-  T[0][0] = +H[GDIM]; T[0][GDIM] = -H[0];
-  T[GDIM][0] = curvature() * H[0]; T[GDIM][GDIM] = +H[GDIM];
+  T[0][0] = +H[LDIM]; T[0][LDIM] = -H[0];
+  T[LDIM][0] = curvature() * H[0]; T[LDIM][LDIM] = +H[LDIM];
   return T;
   }
 
 // reverse of pushxto0(H)
 EX transmatrix rpushxto0(const hyperpoint& H) {
   transmatrix T = Id;
-  T[0][0] = +H[GDIM]; T[0][GDIM] = H[0];
-  T[GDIM][0] = -curvature() * H[0]; T[GDIM][GDIM] = +H[GDIM];
+  T[0][0] = +H[LDIM]; T[0][LDIM] = H[0];
+  T[LDIM][0] = -curvature() * H[0]; T[LDIM][LDIM] = +H[LDIM];
   return T;
   }
 
@@ -597,17 +618,21 @@ EX transmatrix ggpushxto0(const hyperpoint& H, ld co) {
   if(translatable) {
     return eupush(co * H);
     }
+  if(prod) {
+    auto d = product_decompose(H);
+    return mscale(PIU(ggpushxto0(d.second, co)), exp(d.first * co));
+    }
   transmatrix res = Id;
   if(sqhypot_d(GDIM, H) < 1e-12) return res;
-  ld fac = (H[GDIM]-1) / sqhypot_d(GDIM, H);
+  ld fac = (H[LDIM]-1) / sqhypot_d(GDIM, H);
   for(int i=0; i<GDIM; i++)
   for(int j=0; j<GDIM; j++)
     res[i][j] += H[i] * H[j] * fac;
     
   for(int d=0; d<GDIM; d++)
-    res[d][GDIM] = co * H[d],
-    res[GDIM][d] = -curvature() * co * H[d];
-  res[GDIM][GDIM] = H[GDIM];
+    res[d][LDIM] = co * H[d],
+    res[LDIM][d] = -curvature() * co * H[d];
+  res[LDIM][LDIM] = H[LDIM];
   
   return res;
   }
@@ -626,6 +651,7 @@ EX transmatrix rgpushxto0(const hyperpoint& H) {
 
 EX void fixmatrix(transmatrix& T) {
   if(nonisotropic) ; // T may be inverse... do not do that
+  else if(prod) ;
   else if(euclid) {
     for(int x=0; x<GDIM; x++) for(int y=0; y<=x; y++) {
       ld dp = 0;
@@ -635,8 +661,8 @@ EX void fixmatrix(transmatrix& T) {
       
       for(int z=0; z<GDIM; z++) T[z][x] -= dp * T[z][y];
       }
-    for(int x=0; x<GDIM; x++) T[GDIM][x] = 0;
-    T[GDIM][GDIM] = 1;
+    for(int x=0; x<GDIM; x++) T[LDIM][x] = 0;
+    T[LDIM][LDIM] = 1;
     }
   else for(int x=0; x<MDIM; x++) for(int y=0; y<=x; y++) {
     ld dp = 0;
@@ -734,20 +760,29 @@ EX transmatrix inverse(const transmatrix& T) {
     }
   }
 
+EX pair<ld, hyperpoint> product_decompose(hyperpoint h) {
+  ld z = zlevel(h);
+  return make_pair(z, mscale(h, exp(-z)));
+  }
+
 // distance between mh and 0
 EX ld hdist0(const hyperpoint& mh) {
   switch(cgclass) {
     case gcHyperbolic:
-      if(mh[GDIM] < 1) return 0;
-      return acosh(mh[GDIM]);
+      if(mh[LDIM] < 1) return 0;
+      return acosh(mh[LDIM]);
     case gcEuclid: {
       return hypot_d(GDIM, mh);
       }
     case gcSphere: {
-      ld res = mh[GDIM] >= 1 ? 0 : mh[GDIM] <= -1 ? M_PI : acos(mh[GDIM]);
+      ld res = mh[LDIM] >= 1 ? 0 : mh[LDIM] <= -1 ? M_PI : acos(mh[LDIM]);
       if(elliptic && res > M_PI/2) res = M_PI-res;
       return res;
       }
+    case gcProduct: {
+      auto d1 = product_decompose(mh);
+      return hypot(PIU(hdist0(d1.second)), d1.first);
+      }
     default:
       return hypot_d(GDIM, mh);
     }
@@ -779,6 +814,11 @@ EX ld hdist(const hyperpoint& h1, const hyperpoint& h2) {
       return 2 * asinh(sqrt(iv) / 2);
     case gcSphere:
       return 2 * asin_auto_clamp(sqrt(iv) / 2);
+    case gcProduct: {
+      auto d1 = product_decompose(h1);
+      auto d2 = product_decompose(h2);
+      return hypot(PIU(hdist(d1.second, d2.second)), d1.first - d2.first);
+      }
     default:
       if(iv < 0) return 0;
       return sqrt(iv);
@@ -786,7 +826,7 @@ EX ld hdist(const hyperpoint& h1, const hyperpoint& h2) {
   }
 
 EX hyperpoint mscale(const hyperpoint& t, double fac) {
-  if(GDIM == 3) return cpush(2, fac) * t;
+  if(GDIM == 3 && !prod) return cpush(2, fac) * t;
   hyperpoint res;
   for(int i=0; i<MDIM; i++) 
     res[i] = t[i] * fac;
@@ -794,8 +834,8 @@ EX hyperpoint mscale(const hyperpoint& t, double fac) {
   }
 
 EX transmatrix mscale(const transmatrix& t, double fac) {
-  if(GDIM == 3) {
-    // if(pmodel == mdFlatten) { transmatrix u = t; u[2][GDIM] -= fac; return u; }
+  if(GDIM == 3 && !prod) {
+    // if(pmodel == mdFlatten) { transmatrix u = t; u[2][LDIM] -= fac; return u; }
     return t * cpush(2, fac);
     }
   transmatrix res;
@@ -816,7 +856,7 @@ EX transmatrix xyzscale(const transmatrix& t, double fac, double facz) {
   for(int i=0; i<MDIM; i++) for(int j=0; j<GDIM; j++)
     res[i][j] = t[i][j] * fac;
   for(int i=0; i<MDIM; i++) 
-    res[i][GDIM] = t[i][GDIM] * facz;
+    res[i][LDIM] = t[i][LDIM] * facz;
   return res;
   }
 
@@ -870,6 +910,7 @@ EX hyperpoint orthogonal_of_C0(hyperpoint h0, hyperpoint h1, hyperpoint h2) {
 
 EX hyperpoint zshift(hyperpoint x, ld z) {
   if(GDIM == 3 && WDIM == 2) return rgpushxto0(x) * cpush0(2, z);
+  else if(prod) return mscale(x, exp(z));
   else return mscale(x, z);
   }
 
@@ -931,14 +972,14 @@ EX transmatrix transpose(transmatrix T) {
 #if HDR
 inline hyperpoint cpush0(int c, ld x) { 
   hyperpoint h = Hypc;
-  h[GDIM] = cos_auto(x);
+  h[LDIM] = cos_auto(x);
   h[c] = sin_auto(x);
   return h;
   }
 
 inline hyperpoint xspinpush0(ld alpha, ld x) { 
   hyperpoint h = Hypc;
-  h[GDIM] = cos_auto(x);
+  h[LDIM] = cos_auto(x);
   h[0] = sin_auto(x) * cos(alpha);
   h[1] = sin_auto(x) * -sin(alpha);
   return h;
@@ -950,7 +991,7 @@ inline hyperpoint ypush0(ld x) { return cpush0(1, x); }
 // T * C0, optimized
 inline hyperpoint tC0(const transmatrix &T) {
   hyperpoint z;
-  for(int i=0; i<MDIM; i++) z[i] = T[i][GDIM];
+  for(int i=0; i<MDIM; i++) z[i] = T[i][LDIM];
   return z;
   }
 #endif
diff --git a/hypgraph.cpp b/hypgraph.cpp
index df23cc34..1c966cb9 100644
--- a/hypgraph.cpp
+++ b/hypgraph.cpp
@@ -66,7 +66,7 @@ hyperpoint perspective_to_space(hyperpoint h, ld alpha, eGeometryClass gc) {
   hyperpoint H;
   H[0] = hx * (hz+vid.alpha);
   H[1] = hy * (hz+vid.alpha);
-  H[GDIM] = hz;
+  H[LDIM] = hz;
   
   return H;  
   }
@@ -74,7 +74,7 @@ hyperpoint perspective_to_space(hyperpoint h, ld alpha, eGeometryClass gc) {
 hyperpoint space_to_perspective(hyperpoint z, ld alpha = vid.alpha);
 
 hyperpoint space_to_perspective(hyperpoint z, ld alpha) {
-  ld s = 1 / (alpha + z[GDIM]);
+  ld s = 1 / (alpha + z[LDIM]);
   z[0] *= s;
   z[1] *= s;
   if(GDIM == 3) {
@@ -103,7 +103,7 @@ EX hyperpoint gethyper(ld x, ld y) {
 
 void ballmodel(hyperpoint& ret, double alpha, double d, double zl) {
   hyperpoint H = ypush(vid.camera) * xpush(d) * ypush(zl) * C0;
-  ld tzh = vid.ballproj + H[GDIM];
+  ld tzh = vid.ballproj + H[LDIM];
   ld ax = H[0] / tzh;
   ld ay = H[1] / tzh;
   
@@ -166,7 +166,7 @@ ld find_zlev(hyperpoint& H) {
   }
 
 ld get_tz(hyperpoint H) {
-  ld tz = euclid ? (1+vid.alpha) : vid.alpha+H[GDIM];
+  ld tz = euclid ? (1+vid.alpha) : vid.alpha+H[LDIM];
   if(tz < BEHIND_LIMIT && tz > -BEHIND_LIMIT) tz = BEHIND_LIMIT;
   return tz;
   }
@@ -201,8 +201,8 @@ template<class T> void makeband(hyperpoint H, hyperpoint& ret, const T& f) {
   y = asin_auto(H[1]);
   x = asin_auto_clamp(H[0] / cos_auto(y));
   if(sphere) {
-    if(H[GDIM] < 0 && x > 0) x = M_PI - x;
-    else if(H[GDIM] < 0 && x <= 0) x = -M_PI - x;
+    if(H[LDIM] < 0 && x > 0) x = M_PI - x;
+    else if(H[LDIM] < 0 && x <= 0) x = -M_PI - x;
     }
   x += band_shift;
   hypot_zlev(zlev, y, yf, zf);
@@ -495,9 +495,9 @@ EX void applymodel(hyperpoint H, hyperpoint& ret) {
     case mdFisheye: {
       ld zlev = find_zlev(H);
       H = space_to_perspective(H);
-      H[GDIM] = zlev;
+      H[LDIM] = zlev;
       ret = H / sqrt(1 + sqhypot_d(GDIM+1, H));
-      if(GDIM == 3) ret[GDIM] = zlev;
+      if(GDIM == 3) ret[LDIM] = zlev;
       break;
       }
     
@@ -866,11 +866,11 @@ EX bool behindsphere(const hyperpoint& h) {
   if(mdBandAny()) return false;
 
   if(vid.alpha > 1) {
-     if(h[GDIM] > -1/vid.alpha) return true;
+     if(h[LDIM] > -1/vid.alpha) return true;
      }  
   
   if(vid.alpha <= 1) {
-    if(h[GDIM] < .2-vid.alpha) return true;
+    if(h[LDIM] < .2-vid.alpha) return true;
     }
   
   return false;
@@ -948,12 +948,16 @@ EX bool invalid_matrix(const transmatrix T) {
   for(int i=0; i<GDIM; i++) for(int j=0; j<GDIM; j++)
     if(std::isnan(T[i][j]) || T[i][j] > 1e8 || T[i][j] < -1e8 || std::isinf(T[i][j]))
       return true;
-  for(int i=0; i<GDIM; i++) for(int j=0; j<GDIM; j++) if(T[i][j] > .5 || T[i][j] < -.5) return false;
+  if(prod) {
+    for(int i=0; i<GDIM; i++) for(int j=0; j<GDIM; j++) if(abs(T[i][j]) > 1e-6) return false;
+    }
+  else 
+    for(int i=0; i<GDIM; i++) for(int j=0; j<GDIM; j++) if(T[i][j] > .5 || T[i][j] < -.5) return false;
   return true;
   }
   
 EX bool invalid_point(const hyperpoint h) {
-  return std::isnan(h[GDIM]) || h[GDIM] > 1e8 || std::isinf(h[GDIM]);
+  return std::isnan(h[LDIM]) || h[LDIM] > 1e8 || std::isinf(h[LDIM]);
   }
 
 EX bool in_smart_range(const transmatrix& T) {
@@ -1175,25 +1179,25 @@ int mindx=-7, mindy=-7, maxdx=7, maxdy=7;
   
 EX transmatrix eumove(ld x, ld y) {
   transmatrix Mat = Id;
-  Mat[GDIM][GDIM] = 1;
+  Mat[LDIM][LDIM] = 1;
   
   if(a4) {
-    Mat[0][GDIM] += x * cgi.crossf;
-    Mat[1][GDIM] += y * cgi.crossf;
+    Mat[0][LDIM] += x * cgi.crossf;
+    Mat[1][LDIM] += y * cgi.crossf;
     }
   else {
-    Mat[0][GDIM] += (x + y * .5) * cgi.crossf;
-    // Mat[GDIM][0] += (x + y * .5) * cgi.crossf;
-    Mat[1][GDIM] += y * q3 /2 * cgi.crossf;
-    // Mat[GDIM][1] += y * q3 /2 * cgi.crossf;
+    Mat[0][LDIM] += (x + y * .5) * cgi.crossf;
+    // Mat[LDIM][0] += (x + y * .5) * cgi.crossf;
+    Mat[1][LDIM] += y * q3 /2 * cgi.crossf;
+    // Mat[LDIM][1] += y * q3 /2 * cgi.crossf;
     }
   
   ld v = a4 ? 1 : q3;
 
-  while(Mat[0][GDIM] <= -16384 * cgi.crossf) Mat[0][GDIM] += 32768 * cgi.crossf;
-  while(Mat[0][GDIM] >= 16384 * cgi.crossf) Mat[0][GDIM] -= 32768 * cgi.crossf;
-  while(Mat[1][GDIM] <= -16384 * v * cgi.crossf) Mat[1][GDIM] += 32768 * v * cgi.crossf;
-  while(Mat[1][GDIM] >= 16384 * v * cgi.crossf) Mat[1][GDIM] -= 32768 * v * cgi.crossf;
+  while(Mat[0][LDIM] <= -16384 * cgi.crossf) Mat[0][LDIM] += 32768 * cgi.crossf;
+  while(Mat[0][LDIM] >= 16384 * cgi.crossf) Mat[0][LDIM] -= 32768 * cgi.crossf;
+  while(Mat[1][LDIM] <= -16384 * v * cgi.crossf) Mat[1][LDIM] += 32768 * v * cgi.crossf;
+  while(Mat[1][LDIM] >= 16384 * v * cgi.crossf) Mat[1][LDIM] -= 32768 * v * cgi.crossf;
   return Mat;
   }
 
@@ -1320,7 +1324,7 @@ EX void centerpc(ld aspd) {
     }
   #endif
   hyperpoint H = inverse(actual_view_transform) * tC0(T);
-  ld R = zero_d(GDIM, H) ? 0 : hdist0(H);
+  ld R = (zero_d(GDIM, H) && !prod) ? 0 : hdist0(H);
   if(R < 1e-9) {
     // either already centered or direction unknown
     /* if(playerfoundL && playerfoundR) {
@@ -1338,12 +1342,20 @@ EX void centerpc(ld aspd) {
     if(aspd > R) aspd = R;
     
     for(int i=0; i<GDIM; i++)
-      View[i][GDIM] -= H[i] * aspd / R;
+      View[i][LDIM] -= H[i] * aspd / R;
         
     }
   
   else {
     aspd *= (1+R+(shmup::on?1:0));
+    
+    if(prod) {
+      auto d = product_decompose(H);
+      ld dist = -d.first / R * aspd;
+      if(abs(dist) > abs(d.first)) dist = -d.first;
+      View = mscale(View, exp(dist));
+      aspd *= sqrt(R*R - d.first * d.first) / R;
+      }
 
     if(R < aspd) {
       View = solmul(gpushxto0(H), View);
@@ -1363,6 +1375,14 @@ EX void optimizeview() {
   if(subscreens::split(optimizeview)) return;
   if(dual::split(optimizeview)) return;
   
+  if(prod) {
+    ld z = zlevel(tC0(View));
+    View = mscale(View, exp(-z));
+    product::in_underlying_geometry(optimizeview);
+    View = mscale(View, exp(z));
+    return;
+    }
+  
   #if CAP_ANIMATIONS
   if(centerover.at && inmirror(centerover.at)) {
     anims::reflect_view();
@@ -1375,7 +1395,7 @@ EX void optimizeview() {
   
   transmatrix TB = Id;
   
-  if(0) ;
+  if(false) ;
 
   #if CAP_BT || CAP_ARCM || MAXMDIM == 4
   else if(binarytiling || archimedean || penrose || WDIM == 3) {
@@ -1408,7 +1428,7 @@ EX void optimizeview() {
       ld trot = -i * M_PI * 2 / (S7+.0);
       transmatrix T = i < 0 ? Id : spin(trot) * xpush(cgi.tessf) * pispin;
       hyperpoint H = View * tC0(T);
-      if(H[GDIM] < best) best = H[GDIM], turn = i, TB = T;
+      if(H[LDIM] < best) best = H[LDIM], turn = i, TB = T;
       }
   
     if(turn >= 0) {
@@ -1665,7 +1685,7 @@ void queuestraight(hyperpoint X, int style, color_t lc, color_t fc, PPR p) {
 EX void draw_boundary(int w) {
 
   if(w == 1) return;
-  if(nonisotropic || euclid) return;
+  if(nonisotropic || euclid || prod) return;
 
   dynamicval<ld> lw(vid.linewidth, vid.linewidth * vid.multiplier_ring);
 
@@ -1864,7 +1884,7 @@ EX void draw_boundary(int w) {
 
 EX ld band_shift = 0;
 EX void fix_the_band(transmatrix& T) {
-  if(((mdinf[pmodel].flags & mf::uses_bandshift) && T[GDIM][GDIM] > 1e6) || (sphere && pmodel == mdSpiral)) {
+  if(((mdinf[pmodel].flags & mf::uses_bandshift) && T[LDIM][LDIM] > 1e6) || (sphere && pmodel == mdSpiral)) {
     hyperpoint H = tC0(T);
     find_zlev(H);
     models::apply_orientation(H[0], H[1]);
@@ -1872,8 +1892,8 @@ EX void fix_the_band(transmatrix& T) {
     ld y = asin_auto(H[1]);
     ld x = asin_auto_clamp(H[0] / cos_auto(y));
     if(sphere) {
-      if(H[GDIM] < 0 && x > 0) x = M_PI - x;
-      else if(H[GDIM] < 0 && x <= 0) x = -M_PI - x;
+      if(H[LDIM] < 0 && x > 0) x = M_PI - x;
+      else if(H[LDIM] < 0 && x <= 0) x = -M_PI - x;
       }
     band_shift += x;
     T = xpush(-x) * T;
@@ -1936,7 +1956,7 @@ bool limited_generation(cell *c) {
 
 EX bool do_draw(cell *c, const transmatrix& T) {
 
-  PROD( if(product::pmap) return product::in_actual([&] { return do_draw(product::get_at(c, product::plevel), T); }); )
+  if(product::pmap) return product::in_actual([&] { return do_draw(product::get_at(c, product::plevel), T); });
   if(WDIM == 3) {
     if(cells_drawn > vid.cells_drawn_limit) return false;
     if(nil && pmodel == mdGeodesic) {
diff --git a/irregular.cpp b/irregular.cpp
index 1670d7c0..e109273c 100644
--- a/irregular.cpp
+++ b/irregular.cpp
@@ -41,8 +41,8 @@ vector<cellinfo> cells;
 
 ld inner(hyperpoint h1, hyperpoint h2) {
   return 
-    hyperbolic ? h1[GDIM] * h2[GDIM] - h1[0] * h2[0] - h1[1] * h2[1] :
-    h1[GDIM] * h2[GDIM] + h1[0] * h2[0] + h1[1] * h2[1];
+    hyperbolic ? h1[LDIM] * h2[LDIM] - h1[0] * h2[0] - h1[1] * h2[1] :
+    h1[LDIM] * h2[LDIM] + h1[0] * h2[0] + h1[1] * h2[1];
   }
 
 hyperpoint circumscribe(hyperpoint a, hyperpoint b, hyperpoint c) {
@@ -71,7 +71,7 @@ hyperpoint circumscribe(hyperpoint a, hyperpoint b, hyperpoint c) {
     h = h - c * inner(h, c);
     }
   
-  if(h[GDIM] < 0) h[0] = -h[0], h[1] = -h[1], h[GDIM] = -h[GDIM];
+  if(h[LDIM] < 0) h[0] = -h[0], h[1] = -h[1], h[LDIM] = -h[LDIM];
 
   ld i = inner(h, h);
   if(i > 0) h /= sqrt(i);
@@ -328,7 +328,7 @@ bool step(int delta) {
           hyperpoint best_h;
           for(int k=0; k<isize(cells); k++) if(k != i && k != j && k != oldj) {
             hyperpoint h = circumscribe(C0, p1.jpoints[j], p1.jpoints[k]);
-            if(h[GDIM] < 0) continue;
+            if(h[LDIM] < 0) continue;
             if(!clockwise(t, h)) continue;
             if(best_k == -1)
               best_k = k, best_h = h;
@@ -839,7 +839,7 @@ bool save_map(const string& fname) {
     println(f, origcells);
     for(auto i: cells_of_heptagon[h->master]) if(cells[i].generation == 0) {
       auto &ci = cells[i];
-      println(f, spaced(ci.p[0], ci.p[1], ci.p[GDIM]));
+      println(f, spaced(ci.p[0], ci.p[1], ci.p[LDIM]));
       }
     }
   return true;
diff --git a/nonisotropic.cpp b/nonisotropic.cpp
index 0add6b5a..3f939eb5 100644
--- a/nonisotropic.cpp
+++ b/nonisotropic.cpp
@@ -19,7 +19,7 @@ EX namespace nisot {
 
   EX transmatrix translate(hyperpoint h) {
     transmatrix T = Id;
-    for(int i=0; i<GDIM; i++) T[i][GDIM] = h[i];
+    for(int i=0; i<GDIM; i++) T[i][LDIM] = h[i];
     if(sol) {
       T[0][0] = exp(-h[2]);
       T[1][1] = exp(+h[2]);
@@ -541,6 +541,136 @@ EX namespace nilv {
     }
 EX }
 
+EX namespace product {
+
+  EX eGeometry underlying;
+  EX geometry_information *underlying_cgip;
+  
+  void configure() {
+    check_cgi();
+    cgi.prepare_basics();
+    underlying = geometry;
+    underlying_cgip = cgip;
+    geometry = gProduct;
+    ginf[gProduct] = ginf[underlying];
+    ginf[gProduct].sides += 2;
+    ginf[gProduct].cclass = gcProduct;
+    ginf[gProduct].flags |= qEXPERIMENTAL;
+    pmodel = mdPerspective;
+    }
+  
+  EX ld plevel = 1;
+  
+  hrmap *pmap;
+  geometry_information *pcgip;
+  
+  template<class T> auto in_actual(const T& t) {
+    dynamicval<eGeometry> g(geometry, gProduct);
+    dynamicval<geometry_information*> gc(cgip, pcgip);
+    dynamicval<hrmap*> gu(currentmap, pmap);
+    dynamicval<hrmap*> gup(pmap, NULL);
+    return t();  
+    }
+  
+  struct hrmap_product : 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) {
+      pcgip = cgip; pmap = this;
+      dynamicval<eGeometry> g(geometry, underlying);
+      dynamicval<geometry_information*> gc(cgip, underlying_cgip);
+      dynamicval<hrmap*> gu(currentmap, underlying_map);
+      return t();
+      }
+    
+    cell *getCell(cell *u, int h) {
+      cell*& c = at[{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); }
+  
+    transmatrix relative_matrix(cell *c2, cell *c1, const hyperpoint& point_hint) override {
+      return in_underlying([&] { return calc_relative_matrix(where[c2].first, where[c1].first, point_hint); }) * mscale(Id, exp(plevel * (where[c2].second - where[c1].second)));
+      }
+  
+    hrmap_product() {
+      in_underlying([this] { initcells(); underlying_map = currentmap; });
+      }
+  
+    void draw() override {
+      in_underlying([this] { currentmap->draw(); });
+      }
+    };
+  
+  cell *get_at(cell *base, int level) {
+    return ((hrmap_product*)currentmap)->getCell(base, level);
+    }
+  
+  void drawcell_stack(cell *c, transmatrix V, int spinv, bool mirrored) {
+    in_actual([&] { for(int z=-5; z<=5; z++) drawcell(get_at(c, z), V * mscale(Id, exp(plevel * z)), spinv, mirrored); });
+    }
+  
+  void find_cell_connection(cell *c, int d) {
+    auto m = (hrmap_product*)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); });
+      cell *c1 = get_at(cu1, m->where[c].second);
+      c->c.connect(d, c1, cu->c.spin(d), cu->c.mirror(d));
+      }
+    }
+  
+  EX hyperpoint get_corner(cell *c, int i, ld z) {
+    dynamicval<eGeometry> g(geometry, underlying);
+    dynamicval<geometry_information*> gc(cgip, underlying_cgip);
+    return mscale(get_corner_position(c, i), exp(plevel * z/2));
+    }
+  
+  EX hyperpoint where(hyperpoint h) {
+    hyperpoint res;
+    res[2] = zlevel(h);
+    h = mscale(h, exp(-res[2]));
+    ld r = hypot_d(2, h);
+    if(r < 1e-6 || h[2] < 1) {
+      res[0] = h[0];
+      res[1] = h[1];
+      }
+    else {
+      r = acosh(h[2]) / r;
+      res[0] = h[0] * r;
+      res[1] = h[1] * r;
+      }
+    return res;
+    }
+  
+  EX void in_underlying_map(const reaction_t& f) {
+    ((hrmap_product*)currentmap)->in_underlying(f);
+    }
+
+  #if HDR
+  template<class T> auto in_underlying_geometry(const T& f) {
+    dynamicval<eGeometry> g(geometry, underlying);
+    dynamicval<geometry_information*> gc(cgip, underlying_cgip);
+    return f();
+    }
+  
+  #define PIU(x) hr::product::in_underlying_geometry([&] { return (x); })
+  #endif
+  
+EX }
+
 EX namespace nisot {
 
   EX hyperpoint christoffel(const hyperpoint at, const hyperpoint velocity, const hyperpoint transported) {
@@ -657,6 +787,7 @@ EX namespace nisot {
   EX hrmap *new_map() { 
     if(sol) return new solv::hrmap_sol; 
     if(nil) return new nilv::hrmap_nil;
+    if(geometry == gProduct) return new product::hrmap_product;
     return NULL;
     }
   
@@ -686,6 +817,12 @@ EX namespace nisot {
       pmodel = mdPerspective;
       return 0;
       }
+    else if(argis("-product")) {
+      PHASEFROM(2);
+      stop_game();
+      product::configure();
+      return 0;
+      }
     return 1;
     });
     
diff --git a/polygons.cpp b/polygons.cpp
index ce3a92f4..700f09c9 100644
--- a/polygons.cpp
+++ b/polygons.cpp
@@ -718,6 +718,7 @@ void geometry_information::make_wall(int id, vector<hyperpoint> vertices, vector
       hyperpoint h = center + v1 * x + v2 * y;
       if(nil && (x || y))
         h = nilv::on_geodesic(center, nilv::on_geodesic(v1+center, v2+center, y / (x+y)), x + y);
+      if(prod) { hpcpush(h); return; }
       if(sol || !binarytiling) { hpcpush(normalize(h)); return; }
       hyperpoint res = binary::parabolic3(h[0], h[1]) * xpush0(yy*h[2]);
       hpcpush(res);
@@ -729,6 +730,7 @@ void geometry_information::make_wall(int id, vector<hyperpoint> vertices, vector
     int STEP = vid.texture_step;
     for(int a=0; a<n; a++) for(int y=0; y<STEP; y++) {
       hyperpoint h = (vertices[a] * (STEP-y) + vertices[(a+1)%n] * y)/STEP;
+      if(prod) { hpcpush(h); continue; }
       if(nil)
         h = nilv::on_geodesic(vertices[a], vertices[(a+1)%n], y * 1. / STEP);
       if(sol || !binarytiling) { hpcpush(normalize(h)); continue; }
@@ -852,6 +854,20 @@ void geometry_information::create_wall3d() {
     make_wall(12, {point3(3*h,r3,z), point3(0,2*r3,z), point3(-3*h,r3,z)});
     make_wall(13, {point3(3*h,r3,z), point3(3*h,-r3,z), point3(0,-2*r3,z), point3(-3*h,-r3,z), point3(-3*h,r3,z)});
     }
+  
+  if(prod) {
+    cell model;
+    model.type = S7-2;
+    for(int i=0; i<S7-2; i++)
+      make_wall(i, {product::get_corner(&model, i, -1), product::get_corner(&model, i, +1), product::get_corner(&model, i+1, +1), product::get_corner(&model, i+1, -1)});
+    for(int a: {0,1}) {
+      vector<hyperpoint> l;
+      int z = a ? 1 : -1;
+      for(int i=0; i<S7-2; i++)
+        l.push_back(product::get_corner(&model, i, z));
+      make_wall(S7-2+a, l);
+      }    
+    }
 
   if(GDIM == 3 && euclid && S7 == 6) {
     for(int w=0; w<6; w++) {
@@ -907,7 +923,7 @@ void geometry_information::create_wall3d() {
       }
     }
 
-  if(GDIM == 3 && !euclid && !binarytiling && !nil) {
+  if(GDIM == 3 && !euclid && !binarytiling && !nil && !prod) {
     reg3::generate();
     int facesize = isize(reg3::cellshape) / S7;
     for(int w=0; w<S7; w++) {
diff --git a/rug.cpp b/rug.cpp
index 333ac423..57c3f214 100644
--- a/rug.cpp
+++ b/rug.cpp
@@ -326,7 +326,7 @@ EX rugpoint *addRugpoint(hyperpoint h, double dist) {
   
   else {
     m->flat = h;
-    ld hd = h[GDIM];
+    ld hd = h[LDIM];
     for(int d=GDIM; d<rugdim; d++)
       m->flat[d] = (hd - .99) * (rand() % 1000 - rand() % 1000) / 1000;
     }
@@ -997,7 +997,7 @@ bincode acd_bin(ld x) {
 
 bincode get_bincode(hyperpoint h) {
   switch(ginf[gwhere].cclass) {
-    case gcEuclid: case gcSol: case gcNil: PROD( case gcProduct: )
+    case gcEuclid: case gcSol: case gcNil: case gcProduct:
       return acd_bin(h[0]) + acd_bin(h[1]) * sY + acd_bin(h[2]) * sZ;
     case gcHyperbolic:
       return acd_bin(hypot_d(3, h));
diff --git a/sysconfig.h b/sysconfig.h
index 315d3bac..16315b79 100644
--- a/sysconfig.h
+++ b/sysconfig.h
@@ -512,7 +512,4 @@ union SDL_Event;
 #define CAP_MEMORY_RESERVE (!ISMOBILE && !ISWEB)
 #endif
 
-#define PROD(x) /* unimplemented */
-#define PROD2(x,y) /* unimplemented */
-
 #undef TRANSPARENT