diff --git a/binary-tiling.cpp b/binary-tiling.cpp
index 810e6364..24453202 100644
--- a/binary-tiling.cpp
+++ b/binary-tiling.cpp
@@ -4,7 +4,6 @@ namespace hr {
 namespace binary {
 #if CAP_BT
 
-#if DIM == 2
   enum bindir {
     bd_right = 0,
     bd_up_right = 1,
@@ -15,20 +14,16 @@ namespace binary {
     bd_down_left = 5, /* for cells of degree 7 */
     bd_down_right = 6 /* for cells of degree 7 */
     };
-#endif
   
   int type_of(heptagon *h) {
     return h->c7->type;
     }
 
-#if DIM == 2  
   // 0 - central, -1 - left, +1 - right
   int mapside(heptagon *h) {
     return h->zebraval;
     }
-#endif
   
-#if DIM == 2
   #if DEBUG_BINARY_TILING
   map<heptagon*, long long> xcode;
   map<long long, heptagon*> rxcode;
@@ -46,7 +41,6 @@ namespace binary {
     breakhere();
     }
   #endif
-#endif
   
   void breakhere() {
     exit(1);
@@ -83,16 +77,7 @@ namespace binary {
     return h1;
     }
   
-#if DIM == 2
   heptagon *build(heptagon *parent, int d, int d1, int t, int side, int delta) {
-#else
-  heptagon *build(heptagon *parent, int d, int d1, int delta) {
-    int t = 9; 
-    int side = 0;
-    if(d < 4) side = (parent->zebraval * 2 + d) % 5;
-    if(d == 8) side = ((parent->zebraval-d1) * 3) % 5;
-    
-#endif
     auto h = buildHeptagon1(tailored_alloc<heptagon> (t), parent, d, hsOrigin, d1);
     h->distance = parent->distance + delta;
     h->c7 = newCell(t, h);
@@ -108,8 +93,16 @@ namespace binary {
     #endif
     return h;
     }
+
+  #if MAXDIM==4
+  heptagon *build3(heptagon *parent, int d, int d1, int delta) {
+    int side = 0;
+    if(d < 4) side = (parent->zebraval * 2 + d) % 5;
+    if(d == 8) side = ((parent->zebraval-d1) * 3) % 5;
+    return build(parent, d, d1, 9, side, delta);
+    }
+  #endif
   
-  #if DIM == 2
   heptagon *createStep(heptagon *parent, int d) {
     auto h = parent;
     switch(d) {
@@ -163,15 +156,16 @@ namespace binary {
     breakhere();
     return NULL;
     }
-  #else
-  heptagon *createStep(heptagon *parent, int d) {
+
+  #if MAXDIM==4
+  heptagon *createStep3(heptagon *parent, int d) {
     auto h = parent;
     switch(d) {
       case 0: case 1:
       case 2: case 3:
-        return build(parent, d, 8, 1);
+        return build3(parent, d, 8, 1);
       case 8:
-        return build(parent, 8, hrand(4), -1);
+        return build3(parent, 8, hrand(4), -1);
       case 4:
         parent->cmove(8);
         if(parent->c.spin(8) & 1)
@@ -227,7 +221,7 @@ namespace binary {
       return direct_tmatrix[dir];
     }
   
-  #if DIM == 3
+  #if MAXDIM == 4
 
   void queuecube(const transmatrix& V, ld size, color_t linecolor, color_t facecolor) {
     ld yy = log(2) / 2;
@@ -287,29 +281,28 @@ namespace binary {
       
       
       cell *c = h->c7;
-      #if DIM==2
-      if(!do_draw(c, V)) continue;
-      #endif
-      #if DIM==3
-      if(V[DIM][DIM] > btrange_cosh) continue;
-      setdist(c, 7, c);
-      #endif
+      if(DIM == 2 && !do_draw(c, V)) continue;
+      if(DIM == 3) {
+        if(V[DIM][DIM] > btrange_cosh) continue;
+        setdist(c, 7, c);
+        }
       drawcell(c, V, 0, false);
 
-      #if DIM==2
-      dq::enqueue(h->move(bd_up), V * xpush(-log(2)));
-      dq::enqueue(h->move(bd_right), V * parabolic(1));
-      dq::enqueue(h->move(bd_left), V * parabolic(-1));
-      if(c->type == 6)
-        dq::enqueue(h->move(bd_down), V * xpush(log(2)));
-      if(c->type == 7) {
-        dq::enqueue(h->move(bd_down_left), V * parabolic(-1) * xpush(log(2)));
-        dq::enqueue(h->move(bd_down_right), V * parabolic(1) * xpush(log(2)));
+      if(DIM == 2) {
+        dq::enqueue(h->move(bd_up), V * xpush(-log(2)));
+        dq::enqueue(h->move(bd_right), V * parabolic(1));
+        dq::enqueue(h->move(bd_left), V * parabolic(-1));
+        if(c->type == 6)
+          dq::enqueue(h->move(bd_down), V * xpush(log(2)));
+        if(c->type == 7) {
+          dq::enqueue(h->move(bd_down_left), V * parabolic(-1) * xpush(log(2)));
+          dq::enqueue(h->move(bd_down_right), V * parabolic(1) * xpush(log(2)));
+          }
+        }
+      else {
+        for(int i=0; i<9; i++)
+          dq::enqueue(h->move(i), V * tmatrix(h, i));
         }
-      #else
-      for(int i=0; i<9; i++)
-        dq::enqueue(h->move(i), V * tmatrix(h, i));
-      #endif
       }
     }
   
@@ -319,28 +312,28 @@ namespace binary {
     transmatrix gm = Id, where = Id;
     while(h1 != h2) {
       if(h1->distance <= h2->distance) {
-        #if DIM==2
-        if(type_of(h2) == 6)
-          h2 = hr::createStep(h2, bd_down), where = xpush(-log(2)) * where;
-        else if(mapside(h2) == 1)
-          h2 = hr::createStep(h2, bd_left), where = parabolic(+1) * where;
-        else if(mapside(h2) == -1)
-          h2 = hr::createStep(h2, bd_right), where = parabolic(-1) * where;
-        #else
-        h2 = hr::createStep(h2, 8), where = inverse(tmatrix(h2, 8)) * where;
-        #endif
+        if(DIM == 3)
+          h2 = hr::createStep(h2, 8), where = inverse(tmatrix(h2, 8)) * where;
+        else {
+          if(type_of(h2) == 6)
+            h2 = hr::createStep(h2, bd_down), where = xpush(-log(2)) * where;
+          else if(mapside(h2) == 1)
+            h2 = hr::createStep(h2, bd_left), where = parabolic(+1) * where;
+          else if(mapside(h2) == -1)
+            h2 = hr::createStep(h2, bd_right), where = parabolic(-1) * where;
+          }
         }
       else {
-        #if DIM==2
-        if(type_of(h1) == 6)
-          h1 = hr::createStep(h1, bd_down), gm = gm * xpush(log(2));
-        else if(mapside(h1) == 1)
-          h1 = hr::createStep(h1, bd_left), gm = gm * parabolic(-1);
-        else if(mapside(h1) == -1)
-          h1 = hr::createStep(h1, bd_right), gm = gm * parabolic(+1);
-        #else
-        h1 = hr::createStep(h1, 8), where = where * tmatrix(h1, 8);
-        #endif
+        if(DIM == 3)
+          h1 = hr::createStep(h1, 8), where = where * tmatrix(h1, 8);
+        else {
+          if(type_of(h1) == 6)
+            h1 = hr::createStep(h1, bd_down), gm = gm * xpush(log(2));
+          else if(mapside(h1) == 1)
+            h1 = hr::createStep(h1, bd_left), gm = gm * parabolic(-1);
+          else if(mapside(h1) == -1)
+            h1 = hr::createStep(h1, bd_right), gm = gm * parabolic(+1);
+          }
         }
       }
     return gm * where;
@@ -360,9 +353,8 @@ auto bt_config = addHook(hooks_args, 0, [] () {
   return 1;
   });
 #endif
-#endif
 
-int celldistance(cell *c1, cell *c2) {
+int celldistance3(cell *c1, cell *c2) {
   int steps = 0;
   while(c1 != c2) {
     int d1 = celldistAlt(c1), d2 = celldistAlt(c2);
@@ -371,6 +363,8 @@ int celldistance(cell *c1, cell *c2) {
     }
   return steps;
   }
+#endif
+
 
   }
 }
diff --git a/cell.cpp b/cell.cpp
index baa91bcc..226c9bd5 100644
--- a/cell.cpp
+++ b/cell.cpp
@@ -1868,9 +1868,9 @@ int celldistance(cell *c1, cell *c2) {
     return 64;
     }
 
-  #if DIM == 3
-  if(binarytiling) return binary::celldistance(c1, c2);
-  #endif
+  if(binarytiling && DIM == 3) 
+    return binary::celldistance3(c1, c2);
+
   return hyperbolic_celldistance(c1, c2);
   }
 
diff --git a/classes.cpp b/classes.cpp
index 8250d944..4c8d4b11 100644
--- a/classes.cpp
+++ b/classes.cpp
@@ -1757,11 +1757,7 @@ vector<geometryinfo> ginf = {
   {"{8,3}", "Bolza",    "Bolza Surface",                              "Bolza",    8, 3, qsDOCKS,   gcHyperbolic, 0x18200, {{6, 4}}, eVariation::bitruncated},
   {"{8,3}", "Bolza2",   "Bolza Surface x2",                           "Bolza2",   8, 3, qsDOCKS,   gcHyperbolic, 0x18400, {{6, 4}}, eVariation::bitruncated},
   {"{7,3}", "minimal",  "minimal quotient",                           "minimal",  7, 3, qsSMALLN,  gcHyperbolic, 0x18600, {{7, 5}}, eVariation::bitruncated},
-#if DIM == 2
   {"binary","none",     "variant of the binary tiling",               "binary",   7, 3, 0,gcHyperbolic,       0, {{7, 5}}, eVariation::pure},
-#else
-  {"binary","none",     "variant of the binary tiling",               "binary",   9, 3, 0,gcHyperbolic,       0, {{7, 3}}, eVariation::pure},
-#endif
   {"Arch",  "none",     "Archimedean",                                "A",        7, 3, 0,         gcHyperbolic,       0, {{7, 5}}, eVariation::pure},
   {"{7,3}", "Macbeath", "Macbeath Surface",                           "Macbeath", 7, 3, qsSMALL,   gcHyperbolic, 0x20000, {{7, 5}}, eVariation::bitruncated},
   {"{5,4}", "Bring",    "Bring's Surface",                            "Bring",    5, 4, qsSMALL,   gcHyperbolic, 0x20200, {{6, 4}}, eVariation::bitruncated},
@@ -1769,6 +1765,7 @@ vector<geometryinfo> ginf = {
   {"{12,3}","M4",       "Schmutz's M(4)",                             "M4",      12, 3, qsSMALL,   gcHyperbolic, 0x20600, {{4, 2}}, eVariation::bitruncated},
   {"{6,4}", "Crystal",  "dimensional crystal",                        "Crystal",  6, 4, qANYQ,     gcHyperbolic, 0x28000, {{5, 3}}, eVariation::pure},
   {"{3,4}", "none",     "{3,4} (octahedron)",                         "4x3",      3, 4, qsSMALLB,  gcSphere,     0x28200, {{SEE_ALL, SEE_ALL}}, eVariation::bitruncated},
+  {"bin3",  "none",     "3D binary tiling",                           "binary3",  9, 3, 0,         gcHyperbolic,       0, {{7, 3}}, eVariation::pure},
   };
 
 // remember to match the following mask when specifying codes for extra geometries: 0x78600
diff --git a/classes.h b/classes.h
index 5471d213..778cc64d 100644
--- a/classes.h
+++ b/classes.h
@@ -214,7 +214,7 @@ enum eLand { laNone, laBarrier, laCrossroads, laDesert, laIce, laCaves, laJungle
 enum eGeometry {
   gNormal, gEuclid, gSphere, gElliptic, gZebraQuotient, gFieldQuotient, gTorus, gOctagon, g45, g46, g47, gSmallSphere, gTinySphere, gEuclidSquare, gSmallElliptic, 
   gKleinQuartic, gBolza, gBolza2, gMinimal, gBinaryTiling, gArchimedean, 
-  gMacbeath, gBring, gSchmutzM2, gSchmutzM3, gCrystal, gOctahedron,
+  gMacbeath, gBring, gSchmutzM2, gSchmutzM3, gCrystal, gOctahedron, gBinary3,
   gGUARD};
 
 enum eGeometryClass { gcHyperbolic, gcEuclid, gcSphere };
diff --git a/control.cpp b/control.cpp
index 483bff18..2a14065f 100644
--- a/control.cpp
+++ b/control.cpp
@@ -177,7 +177,6 @@ void closeJoysticks() {
   }
 
 void checkjoy() {
-  #if DIM == 2
   DEBB(DF_GRAPH, (debugfile,"check joy\n"));
   if(!DEFAULTCONTROL) return;
   ld joyvalue1 = sqr(vid.joyvalue);
@@ -200,11 +199,9 @@ void checkjoy() {
     }
   
   joydir = vectodir(hpxy(jx, jy));
-  #endif
   }
 
 void checkpanjoy(double t) {
-  #if DIM == 2
   if(shmup::on) return;
   
   if(vid.joypanspeed < 1e-7) return;
@@ -217,7 +214,6 @@ void checkpanjoy(double t) {
   
   playermoved = false;
   View = gpushxto0(hpxy(jx, jy)) * View;
-  #endif
   }
 
 #endif
@@ -305,7 +301,6 @@ void handlePanning(int sym, int uni) {
   if(sym == SDLK_PAGEUP || sym == SDLK_PAGEDOWN) 
     if(isGravityLand(cwt.at->land)) playermoved = false;
 
-  #if DIM == 2
   if(sym == PSEUDOKEY_WHEELUP) {
     ld jx = (mousex - current_display->xcenter - .0) / current_display->radius / 10;
     ld jy = (mousey - current_display->ycenter - .0) / current_display->radius / 10;
@@ -313,7 +308,6 @@ void handlePanning(int sym, int uni) {
     View = gpushxto0(hpxy(jx, jy)) * View;
     sym = 1;
     }
-  #endif
   }
 
 #ifdef SCALETUNER
@@ -587,7 +581,7 @@ void mainloopiter() {
   SDL_Event ev;
   DEBB(DF_GRAPH, (debugfile,"polling for events\n"));
   
-  if(smooth_movement && DEFAULTNOR) {
+  if(smooth_movement) {
     static int lastticks;
     ld t = (ticks - lastticks) * shiftmul / 1000.;
     lastticks = ticks;
@@ -940,11 +934,11 @@ bool haveMobileCompass() {
 #else
   if(forcetarget) return false;
 #endif
+  if(DIM == 3) return false;
   return canmove && !shmup::on && vid.mobilecompasssize > 0 && isize(screens) == 1;
   }
   
 bool handleCompass() {
-#if DIM == 2
   if(!haveMobileCompass()) return false;
   
   using namespace shmupballs;
@@ -966,9 +960,6 @@ bool handleCompass() {
     }
 
   return false;
-#else
-  return false;
-#endif
   }
 
 // orientation sensitivity
diff --git a/floorshapes.cpp b/floorshapes.cpp
index 15979ff7..7efb8978 100644
--- a/floorshapes.cpp
+++ b/floorshapes.cpp
@@ -477,9 +477,7 @@ void generate_floorshapes_for(int id, cell *c, int siid, int sidir) {
 
 void generate_floorshapes() {
 
-  if(0);
-  
-  else if(binarytiling && DIM == 3) ;
+  if(DIM == 3) ;
   
   #if CAP_IRR
   else if(IRREGULAR) {
@@ -625,9 +623,7 @@ void set_floor(const transmatrix& spin, hpcshape& sh) {
 
 void draw_shapevec(cell *c, const transmatrix& V, const vector<hpcshape> &shv, color_t col, PPR prio = PPR::DEFAULT) {
   if(!c) queuepolyat(V, shv[0], col, prio);
-  #if CAP_BT
   else if(DIM == 3) ;
-  #endif
   #if CAP_GP
   else if(GOLDBERG) {
     int id = gp::get_plainshape_id(c);
diff --git a/geom-exp.cpp b/geom-exp.cpp
index a52af220..87ff31f2 100644
--- a/geom-exp.cpp
+++ b/geom-exp.cpp
@@ -371,7 +371,7 @@ void ge_land_selection() {
 vector<eGeometry> tilinglist = {
   gTinySphere, gSmallSphere, gSphere, gEuclid, gNormal, gOctagon,
   gOctahedron, gEuclidSquare, g45, g46, g47,
-  gArchimedean, gBinaryTiling
+  gArchimedean, gBinaryTiling, gBinary3
   };
 
 vector<eGeometry> quotientlist = {
diff --git a/geometry2.cpp b/geometry2.cpp
index 6307d4b3..d5e80e89 100644
--- a/geometry2.cpp
+++ b/geometry2.cpp
@@ -457,7 +457,10 @@ hyperpoint get_corner_position(cell *c, int cid, ld cf) {
   #endif
   #if CAP_BT
   if(binarytiling) {
-    #if DIM == 2
+    if(DIM == 3) {
+      println(hlog, "get_corner_position called");
+      return C0;
+      }
     ld yx = log(2) / 2;
     ld yy = yx;
     ld xx = 1 / sqrt(2)/2;
@@ -470,10 +473,6 @@ hyperpoint get_corner_position(cell *c, int cid, ld cf) {
     vertices[5] = get_horopoint(-yy, -xx);
     vertices[6] = get_horopoint(-yy, 0);
     return mid_at_actual(vertices[cid], 3/cf);
-    #else
-    println(hlog, "get_corner_position called");
-    return C0;
-    #endif
     }
   #endif
   #if CAP_ARCM
@@ -547,7 +546,10 @@ hyperpoint nearcorner(cell *c, int i) {
   #endif
   #if CAP_BT
   if(binarytiling) {
-    #if DIM == 2
+    if(DIM == 3) {
+      println(hlog, "nearcorner called");
+      return Hypc;
+      }
     ld yx = log(2) / 2;
     ld yy = yx;
     // ld xx = 1 / sqrt(2)/2;
@@ -563,10 +565,6 @@ hyperpoint nearcorner(cell *c, int i) {
     else
       neis[5] = get_horopoint(-yy*2, 0);
     return neis[i];
-    #else
-    println(hlog, "nearcorner called");
-    return Hypc;
-    #endif
     }
   #endif
   double d = cellgfxdist(c, i);
diff --git a/graph.cpp b/graph.cpp
index a11fa44e..1bac2d92 100644
--- a/graph.cpp
+++ b/graph.cpp
@@ -5128,8 +5128,7 @@ void drawcell(cell *c, transmatrix V, int spinv, bool mirrored) {
       
       if(0);
       #if CAP_BT
-      else if(binarytiling) {
-        #if DIM == 2
+      else if(binarytiling && DIM == 2) {
         ld yx = log(2) / 2;
         ld yy = yx;
         ld xx = 1 / sqrt(2)/2;
@@ -5143,9 +5142,11 @@ void drawcell(cell *c, transmatrix V, int spinv, bool mirrored) {
         horizontal(yy, 2*xx, xx, 4, binary::bd_up_right);
         horizontal(yy, xx, -xx, 8, binary::bd_up);
         horizontal(yy, -xx, -2*xx, 4, binary::bd_up_left);
-        #else
+        }
+      #endif
+      #if CAP_BT && MAXDIM == 4
+      else if(binarytiling && DIM == 3) {
         binary::queuecube(V, 1, 0xC0C0C080, 0);
-        #endif
         }
       #endif
       else if(isWarped(c) && has_nice_dual()) {
diff --git a/heptagon.cpp b/heptagon.cpp
index 33eeea53..89b6467c 100644
--- a/heptagon.cpp
+++ b/heptagon.cpp
@@ -224,9 +224,13 @@ heptagon *createStep(heptagon *h, int d) {
     crystal::create_step(h, d);
   #endif
   #if CAP_BT
-  if(!h->move(d) && binarytiling) 
+  if(!h->move(d) && binarytiling && DIM == 2)
     return binary::createStep(h, d);
   #endif
+  #if CAP_BT && MAXDIM == 4
+  if(!h->move(d) && binarytiling && DIM == 3) 
+    return binary::createStep3(h, d);
+  #endif
   #if CAP_ARCM
   if(!h->move(d) && archimedean) {
     arcm::create_adjacent(h, d); 
diff --git a/hyper.h b/hyper.h
index 4c7ac6e7..6dec4e90 100644
--- a/hyper.h
+++ b/hyper.h
@@ -94,7 +94,7 @@ void addMessage(string s, char spamtype = 0);
 #define weirdhyperbolic ((S7 > 7 || S3 > 3 || !STDVAR || binarytiling || archimedean) && hyperbolic)
 #define stdhyperbolic (S7 == 7 && S3 == 3 && STDVAR && !binarytiling && !archimedean)
 
-#define binarytiling (geometry == gBinaryTiling)
+#define binarytiling (geometry == gBinaryTiling || geometry == gBinary3)
 #define archimedean (geometry == gArchimedean)
 #define eubinary (euclid || binarytiling || geometry == gCrystal)
 
@@ -189,18 +189,24 @@ typedef complex<ld> cld;
 
 #define DEBSM(x) 
 
+#if MAXDIM == 3
+#define DIM 2
+#else
+#define DIM (geometry == gBinary3 ? 3 : 2)
+#endif
 #define MDIM (DIM+1)
 
-struct hyperpoint : array<ld, MDIM> {
+struct hyperpoint : array<ld, MAXDIM> {
   hyperpoint() {}
+  
   hyperpoint(ld x, ld y, ld z, ld w) { 
     (*this)[0] = x; (*this)[1] = y; (*this)[2] = z; 
-    if(DIM == 3) (*this)[3] = w;
+    if(MAXDIM == 4) (*this)[3] = w;
     }
   };
 
 struct transmatrix {
-  ld tab[MDIM][MDIM];
+  ld tab[MAXDIM][MAXDIM];
   ld * operator [] (int i) { return tab[i]; }
   const ld * operator [] (int i) const { return tab[i]; }
   };
@@ -225,7 +231,7 @@ inline transmatrix operator * (const transmatrix& T, const transmatrix& U) {
   }
 
 constexpr transmatrix diag(ld a, ld b, ld c, ld d) {
-  #if DIM==2
+  #if MAXDIM==3
   return transmatrix{{{a,0,0}, {0,b,0}, {0,0,c}}};
   #else
   return transmatrix{{{a,0,0,0}, {0,b,0,0}, {0,0,c,0}, {0,0,0,d}}};
@@ -254,18 +260,10 @@ const static transmatrix pispin = diag(-1,-1,1,1);
 // central symmetry
 const static transmatrix centralsym = diag(-1,-1,-1,-1);
 
-#if DIM == 3
-inline hyperpoint hpxyz(ld x, ld y, ld z) { return hyperpoint(x,y,0,z); }
-inline hyperpoint hpxyz3(ld x, ld y, ld z, ld w) { return hyperpoint(x,y,z,w); }
+inline hyperpoint hpxyz(ld x, ld y, ld z) { return DIM == 2 ? hyperpoint(x,y,z,0) : hyperpoint(x,y,0,z); }
+inline hyperpoint hpxyz3(ld x, ld y, ld z, ld w) { return DIM == 2 ? hyperpoint(x,y,w,0) : hyperpoint(x,y,z,w); }
 inline hyperpoint point3(ld x, ld y, ld z) { return hyperpoint(x,y,z,0); }
 inline hyperpoint point2(ld x, ld y) { return hyperpoint(x,y,0,0); }
-#else
-#define point3 hpxyz
-inline hyperpoint hpxyz(ld x, ld y, ld z) { return hyperpoint(x,y,z); }
-inline hyperpoint hpxyz3(ld x, ld y, ld z, ld w) { return hyperpoint(x,y,w); }
-inline hyperpoint point2(ld x, ld y) { return hyperpoint(x,y,0); }
-inline hyperpoint point3(ld x, ld y, ld z) { return hyperpoint(x,y,z); }
-#endif
 
 namespace hyperpoint_vec {
 
@@ -2785,7 +2783,14 @@ typedef float GLfloat;
 typedef array<GLfloat, 2> glvec2;
 typedef array<GLfloat, 3> glvec3;
 typedef array<GLfloat, 4> glvec4;
+
+#if MAXDIM == 4
+#define SHDIM 4
 typedef glvec4 glvertex;
+#else
+#define SHDIM 3
+typedef glvec3 glvertex;
+#endif
 
 struct texture_triangle {
   array<hyperpoint, 3> v;
@@ -3787,7 +3792,6 @@ hyperpoint normalize(hyperpoint H);
 extern ld hrandf();
 
 namespace glhr {
-  static const int SHDIM = 4;
 
   struct glmatrix {
     GLfloat a[4][4];
@@ -3811,11 +3815,15 @@ namespace glhr {
   glvertex pointtogl(const hyperpoint& t);
   
   inline glvertex makevertex(GLfloat x, GLfloat y, GLfloat z) {
+    #if SHDIM == 3
+    return glvertex({x,y,z});
+    #else
     return glvertex({x,y,z,1});
+    #endif
     }
   
   struct colored_vertex {
-    glvec3 coords;
+    glvertex coords;
     glvec4 color;
     colored_vertex(GLfloat x, GLfloat y, GLfloat r, GLfloat g, GLfloat b) {
       coords[0] = x;
@@ -3829,12 +3837,12 @@ namespace glhr {
     };
   
   struct textured_vertex {
-    glvec4 coords;
+    glvertex coords;
     glvec2 texture;
     };
   
   struct ct_vertex {
-    glvec4 coords;
+    glvertex coords;
     glvec4 color;
     glvec2 texture;
     ct_vertex(const hyperpoint& h, ld x1, ld y1, ld col) {
@@ -4221,6 +4229,9 @@ hyperpoint get_horopoint3(ld y, ld x, ld z);
 
 namespace binary {
   heptagon *createStep(heptagon *parent, int d);
+  #if MAXDIM == 4
+  heptagon *createStep3(heptagon *parent, int d);
+  #endif
   transmatrix parabolic(ld u);
   transmatrix parabolic3(ld u, ld v);
   extern ld btrange, btrange_cosh;
@@ -4480,7 +4491,7 @@ struct comma_printer {
 template<class T, size_t X> void print(hstream& hs, const array<T, X>& a) { print(hs, "("); comma_printer c(hs); for(const T& t: a) c(t); print(hs, ")"); }
 template<class T> void print(hstream& hs, const vector<T>& a) { print(hs, "("); comma_printer c(hs); for(const T& t: a) c(t); print(hs, ")"); }
 
-inline void print(hstream& hs, const hyperpoint h) { print(hs, (const array<ld, MDIM>&)h); }
+inline void print(hstream& hs, const hyperpoint h) { print(hs, (const array<ld, MAXDIM>&)h); }
 inline void print(hstream& hs, const transmatrix T) { 
   print(hs, "("); comma_printer c(hs);
   for(int i=0; i<MDIM; i++)
diff --git a/hyperpoint.cpp b/hyperpoint.cpp
index 5cc22532..1fea76da 100644
--- a/hyperpoint.cpp
+++ b/hyperpoint.cpp
@@ -4,13 +4,8 @@
 
 namespace hr {
 
-#if DIM == 3
-eGeometry geometry = gBinaryTiling;
-eVariation variation = eVariation::pure;
-#else
 eGeometry geometry;
 eVariation variation;
-#endif
 
 
 // hyperbolic points and matrices
@@ -304,7 +299,7 @@ inline hyperpoint xspinpush0(ld alpha, ld x) {
 transmatrix ypush(ld alpha) { return cpush(1, alpha); }
 
 transmatrix matrix3(ld a, ld b, ld c, ld d, ld e, ld f, ld g, ld h, ld i) {
-  #if DIM==2
+  #if MAXDIM==3
   return transmatrix {{{a,b,c},{d,e,f},{g,h,i}}};
   #else
   return transmatrix {{{a,b,0,c},{d,e,0,f},{0,0,1,0},{g,h,0,i}}};
@@ -312,7 +307,7 @@ transmatrix matrix3(ld a, ld b, ld c, ld d, ld e, ld f, ld g, ld h, ld i) {
   }
 
 transmatrix matrix4(ld a, ld b, ld c, ld d, ld e, ld f, ld g, ld h, ld i, ld j, ld k, ld l, ld m, ld n, ld o, ld p) {
-  #if DIM==2
+  #if MAXDIM==3
   return transmatrix {{{a,b,d},{e,f,h},{m,n,p}}};
   #else
   return transmatrix {{{a,b,c,d},{e,f,g,h},{i,j,k,l},{m,n,o,p}}};
@@ -368,12 +363,9 @@ transmatrix rspintoc(const hyperpoint& H, int t, int f) {
 
 // rotate the hyperbolic plane around C0 such that H[1] == 0 and H[0] >= 0
 transmatrix spintox(const hyperpoint& H) { 
-  #if DIM==2
-  return spintoc(H, 0, 1); 
-  #else
+  if(DIM == 2) return spintoc(H, 0, 1); 
   transmatrix T1 = spintoc(H, 0, 1);
   return spintoc(T1*H, 0, 2) * T1;
-  #endif
   }
 
 void set_column(transmatrix& T, int i, const hyperpoint& H) {
@@ -392,12 +384,9 @@ transmatrix build_matrix(hyperpoint h1, hyperpoint h2, hyperpoint h3) {
 
 // reverse of spintox(H)
 transmatrix rspintox(const hyperpoint& H) { 
-  #if DIM==2
-  return rspintoc(H, 0, 1); 
-  #else
+  if(DIM == 2) return rspintoc(H, 0, 1); 
   transmatrix T1 = spintoc(H, 0, 1);
   return rspintoc(H, 0, 1) * rspintoc(T1*H, 0, 2);
-  #endif
   }
 
 // for H such that H[1] == 0, this matrix pushes H to C0
@@ -474,31 +463,33 @@ void fixmatrix(transmatrix& T) {
 // show the matrix on screen
 
 ld det(const transmatrix& T) {
-  #if DIM == 2
-  ld det = 0;
-  for(int i=0; i<3; i++) 
-    det += T[0][i] * T[1][(i+1)%3] * T[2][(i+2)%3];
-  for(int i=0; i<3; i++) 
-    det -= T[0][i] * T[1][(i+2)%3] * T[2][(i+1)%3];
-  #else
-  ld det = 1;
-  transmatrix M = T;
-  for(int a=0; a<MDIM; a++) {
-    for(int b=a; b<=DIM; b++)
-      if(M[b][a]) {
-        if(b != a)
-          for(int c=a; c<MDIM; c++) tie(M[b][c], M[a][c]) = make_pair(-M[a][c], M[b][c]);
-        break;
-        }
-    if(!M[a][a]) return 0;
-    for(int b=a+1; b<=DIM; b++) {
-      ld co = -M[b][a] / M[a][a];
-      for(int c=a; c<MDIM; c++) M[b][c] += M[a][c] * co;
-      }
-    det *= M[a][a];
+  if(DIM == 2) {
+    ld det = 0;
+    for(int i=0; i<3; i++) 
+      det += T[0][i] * T[1][(i+1)%3] * T[2][(i+2)%3];
+    for(int i=0; i<3; i++) 
+      det -= T[0][i] * T[1][(i+2)%3] * T[2][(i+1)%3];
+    return det;
+    }
+  else {
+    ld det = 1;
+    transmatrix M = T;
+    for(int a=0; a<MDIM; a++) {
+      for(int b=a; b<=DIM; b++)
+        if(M[b][a]) {
+          if(b != a)
+            for(int c=a; c<MDIM; c++) tie(M[b][c], M[a][c]) = make_pair(-M[a][c], M[b][c]);
+          break;
+          }
+      if(!M[a][a]) return 0;
+      for(int b=a+1; b<=DIM; b++) {
+        ld co = -M[b][a] / M[a][a];
+        for(int c=a; c<MDIM; c++) M[b][c] += M[a][c] * co;
+        }
+      det *= M[a][a];
+      }
+    return det;
     }
-  #endif
-  return det;
   }
 
 void inverse_error(const transmatrix& T) {
@@ -506,51 +497,48 @@ void inverse_error(const transmatrix& T) {
   }
 
 transmatrix inverse(const transmatrix& T) {
-  profile_start(7);
-  
-  #if DIM == 2
-  ld d = det(T);
-  transmatrix T2;
-  if(d == 0) {
-    inverse_error(T); 
-    return Id;
+  if(DIM == 2) {
+    ld d = det(T);
+    transmatrix T2;
+    if(d == 0) {
+      inverse_error(T); 
+      return Id;
+      }
+    
+    for(int i=0; i<3; i++) 
+    for(int j=0; j<3; j++) 
+      T2[j][i] = (T[(i+1)%3][(j+1)%3] * T[(i+2)%3][(j+2)%3] - T[(i+1)%3][(j+2)%3] * T[(i+2)%3][(j+1)%3]) / d;
+    return T2;
     }
+  else {
+    transmatrix T1 = T;
+    transmatrix T2 = Id;
   
-  for(int i=0; i<3; i++) 
-  for(int j=0; j<3; j++) 
-    T2[j][i] = (T[(i+1)%3][(j+1)%3] * T[(i+2)%3][(j+2)%3] - T[(i+1)%3][(j+2)%3] * T[(i+2)%3][(j+1)%3]) / d;
-  
-  #else
-  transmatrix T1 = T;
-  transmatrix T2 = Id;
-
-  for(int a=0; a<MDIM; a++) {
-    for(int b=a; b<=DIM; b++)
-      if(T1[b][a]) {
-        if(b != a)
-          for(int c=0; c<MDIM; c++) 
-            swap(T1[b][c], T1[a][c]), swap(T2[b][c], T2[a][c]);
-        break;
+    for(int a=0; a<MDIM; a++) {
+      for(int b=a; b<=DIM; b++)
+        if(T1[b][a]) {
+          if(b != a)
+            for(int c=0; c<MDIM; c++) 
+              swap(T1[b][c], T1[a][c]), swap(T2[b][c], T2[a][c]);
+          break;
+          }
+      if(!T1[a][a]) { inverse_error(T); return Id; }
+      for(int b=a+1; b<=DIM; b++) {
+        ld co = -T1[b][a] / T1[a][a];
+        for(int c=0; c<MDIM; c++) T1[b][c] += T1[a][c] * co, T2[b][c] += T2[a][c] * co;
         }
-    if(!T1[a][a]) { inverse_error(T); return Id; }
-    for(int b=a+1; b<=DIM; b++) {
-      ld co = -T1[b][a] / T1[a][a];
-      for(int c=0; c<MDIM; c++) T1[b][c] += T1[a][c] * co, T2[b][c] += T2[a][c] * co;
       }
-    }
-  
-  for(int a=MDIM-1; a>=0; a--) {
-    for(int b=0; b<a; b++) {
-      ld co = -T1[b][a] / T1[a][a];
-      for(int c=0; c<MDIM; c++) T1[b][c] += T1[a][c] * co, T2[b][c] += T2[a][c] * co;
+    
+    for(int a=MDIM-1; a>=0; a--) {
+      for(int b=0; b<a; b++) {
+        ld co = -T1[b][a] / T1[a][a];
+        for(int c=0; c<MDIM; c++) T1[b][c] += T1[a][c] * co, T2[b][c] += T2[a][c] * co;
+        }
+      ld co = 1 / T1[a][a];
+      for(int c=0; c<MDIM; c++) T1[a][c] *= co, T2[a][c] *= co;
       }
-    ld co = 1 / T1[a][a];
-    for(int c=0; c<MDIM; c++) T1[a][c] *= co, T2[a][c] *= co;
+    return T2;
     }
-  #endif
-
-  profile_stop(7);
-  return T2;
   }
 
 // distance between mh and 0
diff --git a/polygons.cpp b/polygons.cpp
index 3f637749..c59299d1 100644
--- a/polygons.cpp
+++ b/polygons.cpp
@@ -642,16 +642,13 @@ void dqi_poly::gldraw() {
 #endif
 
 double scale_at(const transmatrix& T) {
-  #if DIM==3
-  return 1 / (tC0(T))[2];
-  #else
+  if(DIM == 3) return 1 / (tC0(T))[2];
   using namespace hyperpoint_vec;
   hyperpoint h1, h2, h3;
   applymodel(tC0(T), h1);
   applymodel(T * xpush0(.01), h2);
   applymodel(T * ypush(.01) * C0, h3);
-  return sqrt(hypot2(h2-h1) * hypot2(h3-h1) / .0001);
-  #endif
+  return sqrt(hypot_d(h2-h1, 2) * hypot_d(h3-h1, 2) / .0001);
   }
 
 double linewidthat(const hyperpoint& h) {
diff --git a/shaders.cpp b/shaders.cpp
index 32a3fa8e..f25babcf 100644
--- a/shaders.cpp
+++ b/shaders.cpp
@@ -599,14 +599,15 @@ void init() {
 
 hyperpoint gltopoint(const glvertex& t) {
   hyperpoint h;
-  h[0] = t[0]; h[1] = t[1]; h[2] = t[2]; h[3] = t[3];
-  // if(DIM == 3) h[3] = 0;
+  h[0] = t[0]; h[1] = t[1]; h[2] = t[2]; 
+  if(SHDIM == 4 && MAXDIM == 4) h[3] = t[3];
   return h;
   }
 
 glvertex pointtogl(const hyperpoint& t) {
   glvertex h;
-  h[0] = t[0]; h[1] = t[1]; h[2] = t[2]; h[3] = t[3];
+  h[0] = t[0]; h[1] = t[1]; h[2] = t[2]; 
+  if(SHDIM == 4) h[3] = (DIM == 3) ? t[3] : 1;
   return h;
   }
 
diff --git a/shmup.cpp b/shmup.cpp
index 2fa206e9..ffa0016a 100644
--- a/shmup.cpp
+++ b/shmup.cpp
@@ -845,15 +845,13 @@ void handleInput(int delta) {
       multi::mdx[i] = multi::mdx[i] * (1 - delta / 1000.) + mdx * delta / 2000.;
       multi::mdy[i] = multi::mdy[i] * (1 - delta / 1000.) + mdy * delta / 2000.;
   
-      #if DIM == 2
-      if(mdx != 0 || mdy != 0) if(!multi::combo[i]) {
-        cwtV = multi::whereis[i]; cwt = multi::player[i];
-        flipplayer = multi::flipped[i];
-        multi::whereto[i] = vectodir(hpxy(multi::mdx[i], multi::mdy[i]));
+      if(DIM == 2) {
+        if(mdx != 0 || mdy != 0) if(!multi::combo[i]) {
+          cwtV = multi::whereis[i]; cwt = multi::player[i];
+          flipplayer = multi::flipped[i];
+          multi::whereto[i] = vectodir(hpxy(multi::mdx[i], multi::mdy[i]));
+          }
         }
-      #else
-      ignore(mdx); ignore(mdy);
-      #endif
       
       if(multi::actionspressed[b+pcFire] || 
         (multi::actionspressed[b+pcMoveLeft] && multi::actionspressed[b+pcMoveRight]))
@@ -1626,21 +1624,21 @@ void movePlayer(monster *m, int delta) {
   
   playerturn[cpid] = mturn * delta / 150.0;
 
-  #if DIM == 2
-  double mdd = hypot(mdx, mdy);
-  
-  if(mdd > 1e-6) {
-    hyperpoint jh = hpxy(mdx/100.0, mdy/100.0);
-    hyperpoint ctr = m->pat * C0;
-
-    if(sphere && vid.alpha > 1.001) for(int i=0; i<3; i++) ctr[i] = -ctr[i];
-
-    hyperpoint h = inverse(m->pat) * rgpushxto0(ctr) * jh;
+  if(DIM == 2) {
+    double mdd = hypot(mdx, mdy);
     
-    playerturn[cpid] = -atan2(h[1], h[0]);
-    mgo += mdd;
+    if(mdd > 1e-6) {
+      hyperpoint jh = hpxy(mdx/100.0, mdy/100.0);
+      hyperpoint ctr = m->pat * C0;
+  
+      if(sphere && vid.alpha > 1.001) for(int i=0; i<3; i++) ctr[i] = -ctr[i];
+  
+      hyperpoint h = inverse(m->pat) * rgpushxto0(ctr) * jh;
+      
+      playerturn[cpid] = -atan2(h[1], h[0]);
+      mgo += mdd;
+      }
     }
-  #endif
 
 #if CAP_SDL
   Uint8 *keystate = SDL_GetKeyState(NULL);
diff --git a/sysconfig.h b/sysconfig.h
index b6706d9c..e218a888 100644
--- a/sysconfig.h
+++ b/sysconfig.h
@@ -443,7 +443,7 @@ union SDL_Event;
 #endif
 #endif
 
-#define DIM 3
+#define MAXDIM 4
 
 #ifndef CAP_GEOMETRY
 #define CAP_GEOMETRY (!(ISMINI))
diff --git a/system.cpp b/system.cpp
index d5d34048..924b52be 100644
--- a/system.cpp
+++ b/system.cpp
@@ -1194,7 +1194,7 @@ void set_geometry(eGeometry target) {
     if(DUAL && geometry != gArchimedean) 
       variation = ginf[geometry].default_variation;
     #if CAP_BT
-    if(geometry == gBinaryTiling) variation = eVariation::pure;
+    if(among(geometry, gBinaryTiling, gBinary3)) variation = eVariation::pure;
     #endif
    
     need_reset_geometry = true;