From 56d628d4d8a2bb3c52069158fba2558e4d86102e Mon Sep 17 00:00:00 2001
From: ? <zeno@attnam.com>
Date: Fri, 22 Feb 2019 20:58:40 +0100
Subject: [PATCH] adjusting for 3D part II

---
 basegraph.cpp  |   8 ++--
 geometry.cpp   |   2 +-
 geometry2.cpp  |  10 +++--
 graph.cpp      |   3 +-
 hyper.h        |  86 ++++++++++++++++++-----------------
 hyperpoint.cpp | 118 +++++++++++++++++++++++++++++++------------------
 hypgraph.cpp   |  34 +++++++-------
 polygons.cpp   |  30 ++++++-------
 shaders.cpp    |  43 +++++++++---------
 9 files changed, 188 insertions(+), 146 deletions(-)

diff --git a/basegraph.cpp b/basegraph.cpp
index ee78a6ee..75af207c 100644
--- a/basegraph.cpp
+++ b/basegraph.cpp
@@ -529,9 +529,9 @@ namespace glhr { void texture_vertices(GLfloat *f, int qty, int stride = 2) {
   } 
  void oldvertices(GLfloat *f, int qty) {
   #if CAP_SHADER
-  glVertexAttribPointer(aPosition, 3, GL_FLOAT, GL_FALSE, 3 * sizeof(GLfloat), f);
+  glVertexAttribPointer(aPosition, SHDIM, GL_FLOAT, GL_FALSE, SHDIM * sizeof(GLfloat), f);
   #else
-  glVertexPointer(3, GL_FLOAT, 0, f);
+  glVertexPointer(SHDIM, GL_FLOAT, 0, f);
   #endif
   // #endif
   }
@@ -544,9 +544,9 @@ glhr::textured_vertex charvertex(int x1, int y1, ld tx, ld ty) {
   res.coords[0] = x1;
   res.coords[1] = y1;
   res.coords[2] = 0;
+  res.coords[3] = 1;
   res.texture[0] = tx;
   res.texture[1] = ty;
-  res.texture[2] = 0;
   return res;
   }
 
@@ -965,7 +965,7 @@ void drawCircle(int x, int y, int size, color_t color, color_t fillcolor) {
     if(ISMOBILE && pts > 72) pts = 72;
     for(int r=0; r<pts; r++) {
       float rr = (M_PI * 2 * r) / pts;
-      glcoords.push_back(make_array<GLfloat>(x + size * sin(rr), y + size * vid.stretch * cos(rr), current_display->scrdist));
+      glcoords.push_back(glhr::makevertex(x + size * sin(rr), y + size * vid.stretch * cos(rr), current_display->scrdist));
       }
     glhr::vertices(glcoords);
     glhr::set_depthtest(false);
diff --git a/geometry.cpp b/geometry.cpp
index 59911e56..a0a85434 100644
--- a/geometry.cpp
+++ b/geometry.cpp
@@ -197,7 +197,7 @@ void precalc() {
 
 transmatrix xspinpush(ld dir, ld dist) {
   if(euclid)
-    return eupush(cos(dir) * dist, -sin(dir) * dist DC(,0));
+    return eupush(cos(dir) * dist, -sin(dir));
   else
     return spin(dir) * xpush(dist) * spin(-dir);
   }
diff --git a/geometry2.cpp b/geometry2.cpp
index ea33dac2..6307d4b3 100644
--- a/geometry2.cpp
+++ b/geometry2.cpp
@@ -362,7 +362,11 @@ void virtualRebase(cell*& base, transmatrix& at, bool tohex) {
 
 void virtualRebase(cell*& base, hyperpoint& h, bool tohex) {
   // we perform fixing in check, so that it works with larger range
-  virtualRebase(base, h, tohex, [] (const hyperpoint& h) { return hyperbolic ? hpxy(h[0], h[1] DC(,h[2])) :h; });
+  virtualRebase(base, h, tohex, [] (const hyperpoint& h) { 
+    if(hyperbolic && DIM == 2) return hpxy(h[0], h[1]);
+    if(hyperbolic && DIM == 3) return hpxy3(h[0], h[1], h[2]);
+    return h; 
+    });
   }
 
 // works only in geometries similar to the standard one, and only on heptagons
@@ -436,8 +440,8 @@ hyperpoint randomPointIn(int t) {
   }
 
 #if CAP_BT
-hyperpoint get_horopoint(ld y, ld x DC(,ld z)) {
-  return xpush(-y) * binary::parabolic(x,z) * C0;
+hyperpoint get_horopoint(ld y, ld x) {
+  return xpush(-y) * binary::parabolic(x) * C0;
   }
 #endif
 
diff --git a/graph.cpp b/graph.cpp
index 85115f52..83a49bbd 100644
--- a/graph.cpp
+++ b/graph.cpp
@@ -2286,6 +2286,7 @@ array<array<int,4>,AURA+1> aurac;
 
 bool haveaura() {
   if(!(vid.aurastr>0 && !svg::in && (auraNOGL || vid.usingGL))) return false;
+  if(DIM == 3) return false;
   if(sphere && mdAzimuthalEqui()) return true;
   if(among(pmodel, mdJoukowsky, mdJoukowskyInverted) && hyperbolic && conformal::model_transition < 1) 
     return true;
@@ -5779,7 +5780,7 @@ void drawmovestar(double dx, double dy) {
     Centered = eupush(H);
   else if(R > 1e-9) Centered = rgpushxto0(H);
   
-  Centered = Centered * rgpushxto0(hpxy0(dx*5, dy*5));
+  Centered = Centered * rgpushxto0(hpxy(dx*5, dy*5));
   if(multi::cpid >= 0) multi::crosscenter[multi::cpid] = Centered;
   
   int rax = vid.axes;
diff --git a/hyper.h b/hyper.h
index 048db6a7..d30ec9fb 100644
--- a/hyper.h
+++ b/hyper.h
@@ -191,17 +191,12 @@ typedef complex<ld> cld;
 
 #define MDIM (DIM+1)
 
-#if DIM == 2
-#define D3(x) // x
-#define DC(x,y) // x,y
-#else
-#define D3(x) x
-#define DC(x,y) x,y
-#endif
-
 struct hyperpoint : array<ld, MDIM> {
   hyperpoint() {}
-  hyperpoint(ld x, ld y, ld z DC(, ld w)) { (*this)[0] = x; (*this)[1] = y; (*this)[2] = z; D3((*this)[3] = w;) }
+  hyperpoint(ld x, ld y, ld z, ld w) { 
+    (*this)[0] = x; (*this)[1] = y; (*this)[2] = z; 
+    if(DIM == 3) (*this)[3] = w;
+    }
   };
 
 struct transmatrix {
@@ -221,9 +216,11 @@ inline hyperpoint operator * (const transmatrix& T, const hyperpoint& H) {
 
 inline transmatrix operator * (const transmatrix& T, const transmatrix& U) {
   transmatrix R;
-  // for(int i=0; i<3; i++) for(int j=0; j<3; j++) R[i][j] = 0;
-  for(int i=0; i<MDIM; i++) for(int j=0; j<MDIM; j++) // for(int k=0; k<3; k++)
-    R[i][j] = T[i][0] * U[0][j] + T[i][1] * U[1][j] + T[i][2] * U[2][j] D3(+ T[i][3] * U[3][j]);
+  for(int i=0; i<MDIM; i++) for(int j=0; j<MDIM; j++) {
+    R[i][j] = 0;
+    for(int k=0; k<MDIM; k++)
+      R[i][j] += T[i][k] * U[k][j];
+    }
   return R;
   }
 
@@ -239,6 +236,9 @@ constexpr transmatrix diag(ld a, ld b, ld c, ld d) {
 // identity matrix
 const static transmatrix Id = diag(1,1,1,1);
 
+// zero matrix
+const static transmatrix Zero = diag(0,0,0,0);
+
 // mirror image
 const static transmatrix Mirror = diag(1,-1,1,1);
 
@@ -254,14 +254,17 @@ const static transmatrix pispin = diag(-1,-1,1,1);
 // central symmetry
 const static transmatrix centralsym = diag(-1,-1,-1,-1);
 
-#define hpxyz hyperpoint
-
 #if DIM == 3
-static hyperpoint point3(ld x, ld y, ld z) { return hpxyz(x,y,z,0); }
-static hyperpoint point2(ld x, ld y) { return hpxyz(x,y,0,0); }
+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 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
-static hyperpoint point2(ld x, ld y) { return hpxyz(x,y,0); }
+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 {
@@ -298,13 +301,14 @@ namespace hyperpoint_vec {
       h1[1] * h2[2] - h1[2] * h2[1],
       h1[2] * h2[0] - h1[0] * h2[2],
       h1[0] * h2[1] - h1[1] * h2[0]
-      DC(,0)
       );
     }
 
   // inner product (in R^3)
   inline ld operator | (hyperpoint h1, hyperpoint h2) {
-    return h1[0] * h2[0] + h1[1] * h2[1] + h1[2] * h2[2] D3(+ h1[3] * h2[3]);
+    ld sum = 0;
+    for(int i=0; i<MDIM; i++) sum += h1[i] * h2[i];
+    return sum;
     }    
   }
 
@@ -2778,9 +2782,10 @@ void drawqueue();
 typedef float GLfloat;
 #endif
 
+typedef array<GLfloat, 2> glvec2;
 typedef array<GLfloat, 3> glvec3;
 typedef array<GLfloat, 4> glvec4;
-typedef glvec3 glvertex;
+typedef glvec4 glvertex;
 
 struct texture_triangle {
   array<hyperpoint, 3> v;
@@ -3564,8 +3569,6 @@ template<class T> array<T, 2> make_array(T a, T b) { array<T,2> x; x[0] = a; x[1
 
 extern cell *lastmountpos[MAXPLAYER];
 
-ld hypot3(const hyperpoint& h);
-
 extern const hyperpoint Hypc;
 ld det(const transmatrix& T);
 void queuechr(const hyperpoint& h, int size, char chr, color_t col, int frame = 0);
@@ -3573,20 +3576,15 @@ void queuechr(const hyperpoint& h, int size, char chr, color_t col, int frame =
 string fts(float x);
 bool model_needs_depth();
 
-hyperpoint hpxy(ld x, ld y DC(, ld z));
-ld sqhypot2(const hyperpoint& h);
-ld hypot2(const hyperpoint& h);
+hyperpoint hpxy(ld x, ld y);
+hyperpoint hpxy3(ld x, ld y, ld z);
+ld sqhypot_d(const hyperpoint& h, int d);
+ld hypot_d(const hyperpoint& h, int d);
 transmatrix pushxto0(const hyperpoint& H);
 transmatrix rpushxto0(const hyperpoint& H);
 transmatrix spintox(const hyperpoint& H);
 transmatrix ypush(ld alpha);
 
-#if DIM == 3
-static hyperpoint hpxy0(ld x, ld y) { return hpxy(x, y, 0); }
-#else
-static hyperpoint hpxy0(ld x, ld y) { return hpxy(x, y); }
-#endif
-
 #if CAP_SURFACE
 namespace surface {
 
@@ -3776,10 +3774,10 @@ template<> struct saver<ld> : dsaver<ld> {
 #endif
 extern vector<unique_ptr<drawqueueitem>> ptds;
 extern ld intval(const hyperpoint &h1, const hyperpoint &h2);
-extern ld intvalxy(const hyperpoint &h1, const hyperpoint &h2);
 transmatrix euscalezoom(hyperpoint h);
 transmatrix euaffine(hyperpoint h);
-transmatrix eupush(ld x, ld y DC(, ld z));
+transmatrix eupush(ld x, ld y);
+transmatrix eupush3(ld x, ld y, ld z);
 transmatrix eupush(hyperpoint h);
 transmatrix rspintox(const hyperpoint& H);
 transmatrix gpushxto0(const hyperpoint& H);
@@ -3789,6 +3787,7 @@ hyperpoint normalize(hyperpoint H);
 extern ld hrandf();
 
 namespace glhr {
+  static const int SHDIM = 4;
 
   struct glmatrix {
     GLfloat a[4][4];
@@ -3810,7 +3809,11 @@ namespace glhr {
   void set_modelview(const glmatrix& m);  
   hyperpoint gltopoint(const glvertex& t);
   glvertex pointtogl(const hyperpoint& t);
-
+  
+  inline glvertex makevertex(GLfloat x, GLfloat y, GLfloat z) {
+    return glvertex({x,y,z,1});
+    }
+  
   struct colored_vertex {
     glvec3 coords;
     glvec4 color;
@@ -3826,14 +3829,14 @@ namespace glhr {
     };
   
   struct textured_vertex {
-    glvec3 coords;
-    glvec3 texture;
+    glvec4 coords;
+    glvec2 texture;
     };
   
   struct ct_vertex {
-    glvec3 coords;
+    glvec4 coords;
     glvec4 color;
-    glvec3 texture;
+    glvec2 texture;
     ct_vertex(const hyperpoint& h, ld x1, ld y1, ld col) {
       coords = pointtogl(h);
       texture[0] = x1;
@@ -4213,11 +4216,14 @@ bool saved_tortoise_on(cell *c);
 
 #if CAP_BT
 void horopoint(ld y, ld x);
-hyperpoint get_horopoint(ld y, ld x DC(,ld z));
+hyperpoint get_horopoint(ld y, ld x);
+hyperpoint get_horopoint3(ld y, ld x, ld z);
 
 namespace binary {
   heptagon *createStep(heptagon *parent, int d);
-  transmatrix parabolic(ld u DC(, ld v));
+  transmatrix parabolic(ld u);
+  transmatrix parabolic3(ld u, ld v);
+  extern ld btrange, btrange_cosh;
   }
 #endif
 
diff --git a/hyperpoint.cpp b/hyperpoint.cpp
index 3bb6a45e..0d27a59e 100644
--- a/hyperpoint.cpp
+++ b/hyperpoint.cpp
@@ -124,27 +124,32 @@ ld atan2_auto(ld y, ld x) {
 // by points in 3D space (Minkowski space) such that x^2+y^2-z^2 == -1, z > 0
 // (this is analogous to representing a sphere with points such that x^2+y^2+z^2 == 1)
 
-hyperpoint hpxy(ld x, ld y DC(, ld z)) { 
-  return hpxyz(x,y,DC(z,) euclid ? 1 : sphere ? sqrt(1-x*x-y*y D3(-z*z)) : sqrt(1+x*x+y*y D3(+z*z)));
+hyperpoint hpxy(ld x, ld y) { 
+  return hpxyz(x,y,euclid ? 1 : sphere ? sqrt(1-x*x-y*y) : sqrt(1+x*x+y*y));
+  }
+
+hyperpoint hpxy3(ld x, ld y, ld z) { 
+  return hpxyz3(x,y,z, euclid ? 1 : sphere ? sqrt(1-x*x-y*y-z*z) : sqrt(1+x*x+y*y+z*z));
   }
 
 // center of the pseudosphere
-const hyperpoint Hypc(0,0,0 DC(,0));
+const hyperpoint Hypc = hpxyz(0,0,0);
 
 // origin of the hyperbolic plane
-const hyperpoint C0(0,0 DC(,0) ,1);
+const hyperpoint C0 = hpxyz(0,0,1);
 
 // a point (I hope this number needs no comments ;) )
-const hyperpoint Cx1(1,0 DC(,0) ,1.41421356237);
+const hyperpoint Cx1 = hpxyz(1,0,1.41421356237);
 
 // this function returns approximate square of distance between two points
 // (in the spherical analogy, this would be the distance in the 3D space,
 // through the interior, not on the surface)
 // also used to verify whether a point h1 is on the hyperbolic plane by using Hypc for h2
 
-bool zero2(hyperpoint h) { return h[0] == 0 && h[1] == 0 D3(&& h[2] == 0); }
-
-bool zero3(hyperpoint h) { return h[0] == 0 && h[1] == 0 && h[2] == 0 D3(&& h[3] == 0); }
+bool zero_d(hyperpoint h, int d) { 
+  for(int i=0; i<d; i++) if(h[i]) return false;
+  return true;
+  }
 
 ld intval(const hyperpoint &h1, const hyperpoint &h2) {
   ld res = 0;
@@ -157,28 +162,14 @@ ld intval(const hyperpoint &h1, const hyperpoint &h2) {
   return res;
   }
 
-ld intvalxy(const hyperpoint &h1, const hyperpoint &h2) {
-  return squar(h1[0]-h2[0]) + squar(h1[1]-h2[1]) D3(+ squar(h1[2]-h2[2]));
-  }
-
-ld intvalxyz(const hyperpoint &h1, const hyperpoint &h2) {
-  return squar(h1[0]-h2[0]) + squar(h1[1]-h2[1]) + squar(h1[2]-h2[2]) D3(+ squar(h1[3]-h2[3]));
-  }
-
-ld sqhypot2(const hyperpoint& h) {
-  return h[0]*h[0]+h[1]*h[1] D3(+h[2]*h[2]);
+ld sqhypot_d(const hyperpoint& h, int d) {
+  ld sum = 0;
+  for(int i=0; i<d; i++) sum += h[i]*h[i];
+  return sum;
   }
   
-ld sqhypot3(const hyperpoint& h) {
-  return h[0]*h[0]+h[1]*h[1]+h[2]*h[2] D3(+h[3]*h[3]);
-  }
-  
-ld hypot2(const hyperpoint& h) {
-  return sqrt(sqhypot2(h));
-  }
-  
-ld hypot3(const hyperpoint& h) {
-  return sqrt(sqhypot3(h));
+ld hypot_d(const hyperpoint& h, int d) {
+  return sqrt(sqhypot_d(h, d));
   }
   
 ld zlevel(const hyperpoint &h) {
@@ -242,11 +233,18 @@ transmatrix cspin(int a, int b, ld alpha) {
 
 transmatrix spin(ld alpha) { return cspin(0, 1, alpha); }
 
-transmatrix eupush(ld x, ld y DC(, ld z)) {
+transmatrix eupush(ld x, ld y) {
   transmatrix T = Id;
   T[0][DIM] = x;
   T[1][DIM] = y;
-  D3( T[2][DIM] = z );
+  return T;
+  }
+
+transmatrix eupush3(ld x, ld y, ld z) {
+  transmatrix T = Id;
+  T[0][DIM] = x;
+  T[1][DIM] = y;
+  if(DIM == 3) T[2][DIM] = z;
   return T;
   }
 
@@ -305,16 +303,46 @@ inline hyperpoint xspinpush0(ld alpha, ld x) {
 // push alpha units vertically
 transmatrix ypush(ld alpha) { return cpush(1, alpha); }
 
-transmatrix parabolic1(ld u DC(, ld v)) {
+transmatrix matrix3(ld a, ld b, ld c, ld d, ld e, ld f, ld g, ld h, ld i) {
+  #if DIM==2
+  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}}};
+  #endif
+  }
+
+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
+  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}}};
+  #endif
+  }
+
+transmatrix parabolic1(ld u) {
   if(euclid)
     return ypush(u);
   else {
-    ld diag = (u*u D3(+v*v))/2;
-    #if DIM==2
-    return transmatrix {{{-diag+1, u, diag}, {-u, 1, u}, {-diag, u, diag+1}}};
-    #else
-    return transmatrix {{{-diag+1, u, v, diag}, {-u, 1, 0, u}, {-v, 0, 1, v}, {-diag, u, v, diag+1}}};
-    #endif
+    ld diag = u*u/2;
+    return matrix3(
+      -diag+1, u, diag,
+      -u, 1, u,
+      -diag, u, diag+1
+      );
+    }
+  }
+
+transmatrix parabolic13(ld u, ld v) {
+  if(euclid)
+    return ypush(u);
+  else {
+    ld diag = (u*u+v*v)/2;
+    return matrix4(
+      -diag+1, u, v, diag,
+      -u, 1, 0, u,
+      -v, 0, 1, v,
+      -diag, u, v, diag+1
+      );
     }
   }
 
@@ -354,12 +382,11 @@ void set_column(transmatrix& T, int i, const hyperpoint& H) {
   }
 
 transmatrix build_matrix(hyperpoint h1, hyperpoint h2, hyperpoint h3) {
-  transmatrix T;
+  transmatrix T = Id;
   for(int i=0; i<MDIM; i++)
     T[i][0] = h1[i],
     T[i][1] = h2[i],
-    T[i][2] = h3[i] 
-    DC(, T[i][3] = 0);
+    T[i][2] = h3[i];
   return T;
   }
 
@@ -390,10 +417,13 @@ transmatrix rpushxto0(const hyperpoint& H) {
   }
 
 transmatrix ggpushxto0(const hyperpoint& H, ld co) {
-  if(euclid) return eupush(co * H[0], co * H[1] DC(, co * H[2]));
+  if(euclid) {
+    using namespace hyperpoint_vec;
+    return eupush(co * H);
+    }
   transmatrix res = Id;
-  if(sqhypot2(H) < 1e-12) return res;
-  ld fac = (H[DIM]-1) / sqhypot2(H);
+  if(sqhypot_d(H, DIM) < 1e-12) return res;
+  ld fac = (H[DIM]-1) / sqhypot_d(H, DIM);
   for(int i=0; i<DIM; i++)
   for(int j=0; j<DIM; j++)
     res[i][j] += H[i] * H[j] * fac;
@@ -530,7 +560,7 @@ double hdist0(const hyperpoint& mh) {
       if(mh[DIM] < 1) return 0;
       return acosh(mh[DIM]);
     case gcEuclid: {
-      ld d = sqhypot2(mh);
+      ld d = sqhypot_d(mh, DIM);
       return d;
       }
     case gcSphere: {
diff --git a/hypgraph.cpp b/hypgraph.cpp
index a45b306f..99d05b5e 100644
--- a/hypgraph.cpp
+++ b/hypgraph.cpp
@@ -34,7 +34,7 @@ hyperpoint perspective_to_space(hyperpoint h, ld alpha, eGeometryClass gc) {
   ld hx = h[0], hy = h[1];
   
   if(gc == gcEuclid)
-    return hpxy0(hx * (1 + alpha), hy * (1 + alpha));
+    return hpxy(hx * (1 + alpha), hy * (1 + alpha));
     
   ld hr = hx*hx+hy*hy;
   
@@ -85,7 +85,7 @@ hyperpoint gethyper(ld x, ld y) {
   
   if(vid.camera_angle) camrotate(hx, hy);
   
-  return perspective_to_space(hpxyz(hx, hy, 0 DC(,0)));
+  return perspective_to_space(hpxyz(hx, hy, 0));
   }
 
 void ballmodel(hyperpoint& ret, double alpha, double d, double zl) {
@@ -180,7 +180,7 @@ template<class T> void makeband(hyperpoint H, hyperpoint& ret, const T& f) {
   
   ld yzf = y * zf; y *= yf;
   conformal::apply_orientation(y, x);
-  ret = hpxyz(x / M_PI, y / M_PI, 0 DC(,0));
+  ret = hpxyz(x / M_PI, y / M_PI, 0);
   if(zlev != 1 && current_display->stereo_active()) 
     apply_depth(ret, yzf / M_PI);
   return;
@@ -288,7 +288,7 @@ void applymodel(hyperpoint H, hyperpoint& ret) {
       conformal::apply_orientation(H[0], H[1]);
   
       H[1] += 1;
-      double rad = sqhypot2(H);
+      double rad = sqhypot_d(H, 2);
       H /= -rad;
       H[1] += .5;
       
@@ -318,7 +318,7 @@ void applymodel(hyperpoint H, hyperpoint& ret) {
         case gcHyperbolic: {
           ld zl = zlevel(H);
           ret = H / H[2];
-          ret[2] = sqrt(1 - sqhypot2(ret));
+          ret[2] = sqrt(1 - sqhypot_d(ret, 2));
           ret = ret * (1 + (zl - 1) * ret[2]);
           break;
           }
@@ -326,7 +326,7 @@ void applymodel(hyperpoint H, hyperpoint& ret) {
         case gcEuclid: {
           // stereographic projection to a sphere
           auto hd = hdist0(H) / vid.euclid_to_sphere;
-          if(hd == 0) ret = hpxyz(0, 0, -1 DC(,0));
+          if(hd == 0) ret = hpxyz(0, 0, -1);
           else {
             ld x = 2 * hd / (1 + hd * hd);
             ld y = x / hd;
@@ -356,7 +356,7 @@ void applymodel(hyperpoint H, hyperpoint& ret) {
       if(pmodel == mdHyperboloid) {
         ld& topz = conformal::top_z;
         if(H[2] > topz) {
-          ld scale = sqrt(topz*topz-1) / hypot2(H);
+          ld scale = sqrt(topz*topz-1) / hypot_d(H, 2);
           H *= scale;
           H[2] = topz;
           }
@@ -378,7 +378,7 @@ void applymodel(hyperpoint H, hyperpoint& ret) {
       ld zlev = find_zlev(H);
       H = space_to_perspective(H);
       H[2] = zlev;
-      ret = H / sqrt(1 + sqhypot3(H));
+      ret = H / sqrt(1 + sqhypot_d(H, 3));
       break;
       }
     
@@ -398,7 +398,7 @@ void applymodel(hyperpoint H, hyperpoint& ret) {
           cld w(ret[0], ret[1]);
           cld z = sqrt(c4*w*w-c1) + c2*w;
           if(abs(z) > 1) z = c1 / z;
-          hyperpoint zr = hpxyz(real(z), imag(z), 0 DC(,0));
+          hyperpoint zr = hpxyz(real(z), imag(z), 0);
           
           hyperpoint inhyp = perspective_to_space(zr, 1, gcHyperbolic);
           last_skiprope = vid.skiprope;
@@ -408,7 +408,7 @@ void applymodel(hyperpoint H, hyperpoint& ret) {
         }
   
       H = space_to_perspective(H);
-      ld r = hypot2(H);
+      ld r = hypot_d(H, 2);
       ld c = H[0] / r;
       ld s = H[1] / r;
       ld& mt = conformal::model_transition;
@@ -423,7 +423,7 @@ void applymodel(hyperpoint H, hyperpoint& ret) {
         ret = mobius(ret, vid.skiprope, 2);
       
       if(pmodel == mdJoukowskyInverted) {
-        ld r2 = sqhypot2(ret);
+        ld r2 = sqhypot_d(ret, 2);
         ret[0] = ret[0] / r2;
         ret[1] = -ret[1] / r2;
         conformal::apply_orientation(ret[1], ret[0]);
@@ -473,7 +473,7 @@ void applymodel(hyperpoint H, hyperpoint& ret) {
         H[0] -= .5;
         
         ld phi = atan2(H);
-        ld r = hypot2(H);
+        ld r = hypot_d(H, 2);
         
         r = pow(r, 1 - mt);
         phi *= (1 - mt);
@@ -510,7 +510,7 @@ void applymodel(hyperpoint H, hyperpoint& ret) {
     case mdEquidistant: case mdEquiarea: {
       ld zlev = find_zlev(H);
 
-      ld rad = hypot2(H);
+      ld rad = hypot_d(H, 2);
       if(rad == 0) rad = 1;
       ld d = hdist0(H);
       ld df, zf;
@@ -1051,7 +1051,7 @@ void centerpc(ld aspd) {
   ors::unrotate(cwtV); ors::unrotate(View);
   
   hyperpoint H = ypush(-vid.yshift) * sphereflip * tC0(cwtV);
-  ld R = zero2(H) ? 0 : hdist0(H); // = sqrt(H[0] * H[0] + H[1] * H[1]);
+  ld R = zero_d(H, 2) ? 0 : hdist0(H); // = sqrt(H[0] * H[0] + H[1] * H[1]);
   if(R < 1e-9) {
     // either already centered or direction unknown
     /* if(playerfoundL && playerfoundR) {
@@ -1232,7 +1232,7 @@ void circle_around_center(ld radius, color_t linecol, color_t fillcol, PPR prio)
   if(among(pmodel, mdDisk, mdEquiarea, mdEquidistant, mdFisheye) && !(pmodel == mdDisk && hyperbolic && vid.alpha <= -1) && vid.camera_angle == 0) {
     hyperpoint ret;
     applymodel(xpush0(radius), ret);
-    ld r = hypot2(ret);
+    ld r = hypot_d(ret, 2);
     queuecircle(current_display->xcenter, current_display->ycenter, r * current_display->radius, linecol, prio, fillcol);
     return;
     }  
@@ -1328,7 +1328,7 @@ void queuestraight(hyperpoint X, int style, color_t lc, color_t fc, PPR p) {
   hyperpoint H;
   applymodel(X, H);
   H *= current_display->radius;
-  ld mul = hypot(vid.xres, vid.yres) / hypot2(H);
+  ld mul = hypot(vid.xres, vid.yres) / hypot_d(H, 2);
   ld m = style == 1 ? -mul : -1;
   
   queuereset(mdUnchanged, p);
@@ -1385,7 +1385,7 @@ void draw_boundary(int w) {
         ld z = -sqrt(1 - x*x);
         conformal::apply_orientation(y, x);
         hyperpoint h1;
-        applymodel(hpxyz(x,y,DC(0,) z), h1);
+        applymodel(hpxyz(x,y,z), h1);
         
         conformal::apply_orientation(h1[0], h1[1]);      
         h1[1] = abs(h1[1]) * b;
diff --git a/polygons.cpp b/polygons.cpp
index 9c8c3ffd..7733e6b5 100644
--- a/polygons.cpp
+++ b/polygons.cpp
@@ -220,7 +220,7 @@ SDL_Surface *aux;
 #if CAP_POLY
 #define POLYMAX 60000
 
-vector<array<GLfloat, 3>> glcoords, ourshape;
+vector<glvertex> glcoords, ourshape;
 
 void initPolyForGL() {
 
@@ -252,7 +252,7 @@ int polyx[POLYMAX], polyxr[POLYMAX], polyy[POLYMAX];
 int poly_flags;
 
 void add1(const hyperpoint& H) {
-  glcoords.push_back(make_array<GLfloat>(H[0], H[1], H[2])); 
+  glcoords.push_back(glhr::pointtogl(H)); 
   }  
 
 bool is_behind(const hyperpoint& H) {
@@ -314,7 +314,7 @@ bool correct_side(const hyperpoint& H) {
   return get_side(H) == spherespecial;
   }
 
-void fixpoint(array<float, 3>& hscr, hyperpoint H) {
+void fixpoint(glvertex& hscr, hyperpoint H) {
   hyperpoint bad = H, good = goodpoint;
 
   for(int i=0; i<10; i++) {
@@ -326,7 +326,7 @@ void fixpoint(array<float, 3>& hscr, hyperpoint H) {
     }
   hyperpoint Hscr;
   applymodel(good, Hscr); 
-  hscr = make_array<GLfloat>(Hscr[0]*current_display->radius, Hscr[1]*current_display->radius*vid.stretch, Hscr[2]*current_display->radius); 
+  hscr = glhr::makevertex(Hscr[0]*current_display->radius, Hscr[1]*current_display->radius*vid.stretch, Hscr[2]*current_display->radius); 
   }
 
 void addpoint(const hyperpoint& H) {
@@ -598,10 +598,10 @@ void dqi_poly::gldraw() {
           GLfloat yy = vid.yres;
           GLfloat dist = shaderside_projection ? current_display->scrdist : 0;
           vector<glvertex> scr = {
-            make_array<GLfloat>(-xx, -yy, dist), 
-            make_array<GLfloat>(+xx, -yy, dist), 
-            make_array<GLfloat>(+xx, +yy, dist), 
-            make_array<GLfloat>(-xx, +yy, dist)
+            glhr::makevertex(-xx, -yy, dist), 
+            glhr::makevertex(+xx, -yy, dist), 
+            glhr::makevertex(+xx, +yy, dist), 
+            glhr::makevertex(-xx, +yy, dist)
             };
           glhr::vertices(scr);
           glhr::id_modelview();
@@ -834,11 +834,11 @@ void compute_side_by_centerin(dqi_poly *p, bool& nofill) {
   
   /*
   if(poly_flags & POLY_BADCENTERIN) {
-    glcoords.push_back(make_array<GLfloat>(hscr[0]+10, hscr[1]*vid.stretch, hscr[2]));
-    glcoords.push_back(make_array<GLfloat>(hscr[0], hscr[1]*vid.stretch+10, hscr[2]));
-    glcoords.push_back(make_array<GLfloat>(hscr[0]-10, hscr[1]*vid.stretch, hscr[2]));
-    glcoords.push_back(make_array<GLfloat>(hscr[0], hscr[1]*vid.stretch-10, hscr[2]));
-    glcoords.push_back(make_array<GLfloat>(hscr[0]+10, hscr[1]*vid.stretch, hscr[2]));
+    glcoords.push_back(glhr::makevertex(hscr[0]+10, hscr[1]*vid.stretch, hscr[2]));
+    glcoords.push_back(glhr::makevertex(hscr[0], hscr[1]*vid.stretch+10, hscr[2]));
+    glcoords.push_back(glhr::makevertex(hscr[0]-10, hscr[1]*vid.stretch, hscr[2]));
+    glcoords.push_back(glhr::makevertex(hscr[0], hscr[1]*vid.stretch-10, hscr[2]));
+    glcoords.push_back(glhr::makevertex(hscr[0]+10, hscr[1]*vid.stretch, hscr[2]));
     } */
   }
 
@@ -949,7 +949,7 @@ void dqi_poly::draw() {
           ld c1 = ah1[1], c2 = -ah2[1];
           if(c1 < 0) c1 = -c1, c2 = -c2;
           hyperpoint h = ah1 * c1 + ah2 * c2;
-          h /= hypot3(h);
+          h /= hypot_d(h, 3);
           if(h[2] < 0 && abs(h[0]) < sin(vid.twopoint_param)) cpha = 1-cpha, pha = 2;
           }
         if(cpha == 1) pha = 0;
@@ -1098,7 +1098,7 @@ void dqi_poly::draw() {
         ld h = atan2(glcoords[0][0], glcoords[0][1]);
         for(int i=0; i<=360; i++) {
           ld a = i * degree + h;
-          glcoords.push_back(make_array<GLfloat>(current_display->radius * sin(a), current_display->radius * vid.stretch * cos(a), current_display->scrdist));
+          glcoords.push_back(glhr::makevertex(current_display->radius * sin(a), current_display->radius * vid.stretch * cos(a), current_display->scrdist));
           }
         poly_flags ^= POLY_INVERSE;
         }
diff --git a/shaders.cpp b/shaders.cpp
index e1797e94..5085b549 100644
--- a/shaders.cpp
+++ b/shaders.cpp
@@ -590,13 +590,14 @@ void init() {
 
 hyperpoint gltopoint(const glvertex& t) {
   hyperpoint h;
-  h[0] = t[0]; h[1] = t[1]; h[2] = t[2];
+  h[0] = t[0]; h[1] = t[1]; h[2] = t[2]; h[3] = t[3];
+  // if(DIM == 3) h[3] = 0;
   return h;
   }
 
 glvertex pointtogl(const hyperpoint& t) {
   glvertex h;
-  h[0] = t[0]; h[1] = t[1]; h[2] = t[2];
+  h[0] = t[0]; h[1] = t[1]; h[2] = t[2]; h[3] = t[3];
   return h;
   }
 
@@ -620,16 +621,16 @@ void vertices(const vector<glvertex>& v) {
     if(&v[0] == current_vertices) return;
     current_vertices = buffered_vertices;
     glBindBuffer(GL_ARRAY_BUFFER, buf_buffered);
-    glVertexAttribPointer(glhr::aPosition, 3, GL_FLOAT, GL_FALSE, sizeof(glvertex), 0);
+    glVertexAttribPointer(glhr::aPosition, SHDIM, GL_FLOAT, GL_FALSE, sizeof(glvertex), 0);
     return;
     }
   bindbuffer(v);
-  glVertexAttribPointer(glhr::aPosition, 3, GL_FLOAT, GL_FALSE, sizeof(glvertex), 0);
+  glVertexAttribPointer(glhr::aPosition, SHDIM, GL_FLOAT, GL_FALSE, sizeof(glvertex), 0);
   #else
   if(current_vertices == &v[0]) return;
   current_vertices = &v[0];
   #if CAP_SHADER
-  glVertexAttribPointer(aPosition, 3, GL_FLOAT, GL_FALSE, sizeof(glvertex), &v[0]);
+  glVertexAttribPointer(aPosition, SHDIM, GL_FLOAT, GL_FALSE, sizeof(glvertex), &v[0]);
   #else
   glVertexPointer(3, GL_FLOAT, sizeof(glvertex), &v[0]);
   #endif
@@ -642,9 +643,9 @@ void vertices_texture(const vector<glvertex>& v, const vector<glvertex>& t) {
   #else
   vertices(v);
   #if CAP_SHADER
-  glVertexAttribPointer(aTexture, 3, GL_FLOAT, GL_FALSE, sizeof(glvertex), &t[0]);
+  glVertexAttribPointer(aTexture, SHDIM, GL_FLOAT, GL_FALSE, sizeof(glvertex), &t[0]);
   #else
-  glTexCoordPointer(3, GL_FLOAT, 0, &v[0]);
+  glTexCoordPointer(SHDIM, GL_FLOAT, 0, &v[0]);
   #endif
   #endif
   }
@@ -652,16 +653,16 @@ void vertices_texture(const vector<glvertex>& v, const vector<glvertex>& t) {
 void prepare(vector<colored_vertex>& v) {
   #if CAP_VERTEXBUFFER
   bindbuffer(v);
-  PTR(glhr::aPosition, 3, coords);
+  PTR(glhr::aPosition, SHDIM, coords);
   PTR(glhr::aColor, 4, color);
   #else
   if(current_vertices == &v[0]) return;
   current_vertices = &v[0];
   #if CAP_SHADER
-  glVertexAttribPointer(aPosition, 3, GL_FLOAT, GL_FALSE, sizeof(colored_vertex), &v[0].coords);
+  glVertexAttribPointer(aPosition, SHDIM, GL_FLOAT, GL_FALSE, sizeof(colored_vertex), &v[0].coords);
   glVertexAttribPointer(aColor, 4, GL_FLOAT, GL_FALSE, sizeof(colored_vertex), &v[0].color);
   #else
-  glVertexPointer(3, GL_FLOAT, sizeof(colored_vertex), &v[0].coords);
+  glVertexPointer(SHDIM, GL_FLOAT, sizeof(colored_vertex), &v[0].coords);
   glColorPointer(4, GL_FLOAT, sizeof(colored_vertex), &v[0].color);
   #endif
   #endif
@@ -670,17 +671,17 @@ void prepare(vector<colored_vertex>& v) {
 void prepare(vector<textured_vertex>& v) {
   #if CAP_VERTEXBUFFER
   bindbuffer(v);
-  PTR(glhr::aPosition, 3, coords);
+  PTR(glhr::aPosition, SHDIM, coords);
   PTR(glhr::aTexture, 2, texture);
   #else
   if(current_vertices == &v[0]) return;
   current_vertices = &v[0];
   #if CAP_SHADER
-  glVertexAttribPointer(aPosition, 3, GL_FLOAT, GL_FALSE, sizeof(textured_vertex), &v[0].coords);
-  glVertexAttribPointer(aTexture, 3, GL_FLOAT, GL_FALSE, sizeof(textured_vertex), &v[0].texture);
+  glVertexAttribPointer(aPosition, SHDIM, GL_FLOAT, GL_FALSE, sizeof(textured_vertex), &v[0].coords);
+  glVertexAttribPointer(aTexture, SHDIM, GL_FLOAT, GL_FALSE, sizeof(textured_vertex), &v[0].texture);
   #else
-  glVertexPointer(3, GL_FLOAT, sizeof(textured_vertex), &v[0].coords);
-  glTexCoordPointer(3, GL_FLOAT, sizeof(textured_vertex), &v[0].texture);
+  glVertexPointer(SHDIM, GL_FLOAT, sizeof(textured_vertex), &v[0].coords);
+  glTexCoordPointer(2, GL_FLOAT, sizeof(textured_vertex), &v[0].texture);
   #endif
   #endif
   // color2(col);
@@ -689,19 +690,19 @@ void prepare(vector<textured_vertex>& v) {
 void prepare(vector<ct_vertex>& v) {
   #if CAP_VERTEXBUFFER
   bindbuffer(v);
-  PTR(glhr::aPosition, 3, coords);
+  PTR(glhr::aPosition, SHDIM, coords);
   PTR(glhr::aColor, 4, color);
   PTR(glhr::aTexture, 2, texture);
   #else
   if(current_vertices == &v[0]) return;
   current_vertices = &v[0];
   #if CAP_SHADER
-  glVertexAttribPointer(aPosition, 3, GL_FLOAT, GL_FALSE, sizeof(ct_vertex), &v[0].coords);
+  glVertexAttribPointer(aPosition, SHDIM, GL_FLOAT, GL_FALSE, sizeof(ct_vertex), &v[0].coords);
   glVertexAttribPointer(aColor, 4, GL_FLOAT, GL_FALSE, sizeof(ct_vertex), &v[0].color);
-  glVertexAttribPointer(aTexture, 3, GL_FLOAT, GL_FALSE, sizeof(ct_vertex), &v[0].texture);
+  glVertexAttribPointer(aTexture, 2, GL_FLOAT, GL_FALSE, sizeof(ct_vertex), &v[0].texture);
   #else
-  glVertexPointer(3, GL_FLOAT, sizeof(ct_vertex), &v[0].coords);
-  glTexCoordPointer(3, GL_FLOAT, sizeof(ct_vertex), &v[0].texture);
+  glVertexPointer(SHDIM, GL_FLOAT, sizeof(ct_vertex), &v[0].coords);
+  glTexCoordPointer(2, GL_FLOAT, sizeof(ct_vertex), &v[0].texture);
   glColorPointer(4, GL_FLOAT, sizeof(ct_vertex), &v[0].color);
   #endif
   #endif
@@ -716,7 +717,7 @@ void store_in_buffer(vector<glvertex>& v) {
   printf("storing %d in buffer: %p\n", isize(v), &v[0]);
   current_vertices = buffered_vertices = &v[0];
   glBindBuffer(GL_ARRAY_BUFFER, buf_buffered);
-  glVertexAttribPointer(glhr::aPosition, 3, GL_FLOAT, GL_FALSE, sizeof(glvertex), 0);
+  glVertexAttribPointer(glhr::aPosition, SHDIM, GL_FLOAT, GL_FALSE, sizeof(glvertex), 0);
   glBufferData(GL_ARRAY_BUFFER, isize(v) * sizeof(glvertex), &v[0], GL_STATIC_DRAW);
   printf("Stored.\n");
 #endif