From fd3b030119f566969e4faac09e7495449f7f7366 Mon Sep 17 00:00:00 2001
From: Zeno Rogue <zeno@attnam.com>
Date: Fri, 30 Nov 2018 16:31:55 +0100
Subject: [PATCH] crystal:: added multi-dimensional geometries (referred to as
 crystal), available via CLI for now

---
 cell.cpp       |   6 +
 classes.cpp    |   1 +
 classes.h      |   2 +-
 compileunits.h |   1 +
 crystal.cpp    | 457 +++++++++++++++++++++++++++++++++++++++++++++++++
 geometry2.cpp  |   4 +-
 heptagon.cpp   |   4 +-
 hyper.h        |   7 +
 pattern2.cpp   |   5 +
 9 files changed, 483 insertions(+), 4 deletions(-)
 create mode 100644 crystal.cpp

diff --git a/cell.cpp b/cell.cpp
index 3924641c..d614177f 100644
--- a/cell.cpp
+++ b/cell.cpp
@@ -1170,6 +1170,7 @@ void initcells() {
   
   hrmap* res = callhandlers((hrmap*)nullptr, hooks_newmap);
   if(res) currentmap = res;  
+  else if(geometry == gCrystal) currentmap = crystal::new_map();
   else if(archimedean) currentmap = arcm::new_map();
   else if(fulltorus) currentmap = new hrmap_torus;
   else if(euclid) currentmap = new hrmap_euclidean;
@@ -1333,6 +1334,8 @@ int compdist(int dx[]) {
 int celldist(cell *c) {
   if(fulltorus) 
     return torusmap()->dists[decodeId(c->master)];
+  if(geometry == gCrystal)
+    return crystal::distance(c, currentmap->gamestart());
   if(euwrap)
     return torusconfig::cyldist(decodeId(c->master), 0);
   if(masterless)
@@ -1734,6 +1737,9 @@ int celldistance(cell *c1, cell *c2) {
   
   if(geometry == gFieldQuotient && !GOLDBERG)
     return currfp.getdist(fieldpattern::fieldval(c1), fieldpattern::fieldval(c2));
+  
+  if(geometry == gCrystal)
+    return crystal::distance(c1, c2);
 
   if(bounded) {
     
diff --git a/classes.cpp b/classes.cpp
index a5861ebd..71d3f650 100644
--- a/classes.cpp
+++ b/classes.cpp
@@ -1693,6 +1693,7 @@ vector<geometryinfo> ginf = {
   {"Bring's Surface",     "Bring",    5, 4, qsSMALL,   gcHyperbolic, 0x20200, {{6, 4}}, eVariation::bitruncated},
   {"Schmutz's M(3)",      "M3",      12, 3, qsSMALL,   gcHyperbolic, 0x20400, {{4, 2}}, eVariation::bitruncated},
   {"Schmutz's M(4)",      "M4",      12, 3, qsSMALL,   gcHyperbolic, 0x20600, {{4, 2}}, eVariation::bitruncated},
+  {"Crystal",             "Crystal",  8, 4, qANYQ,     gcHyperbolic,     0x28000, {{3, 2}}, eVariation::pure},
   };
 
 // remember to match the following mask when specifying codes for extra geometries: 0x78600
diff --git a/classes.h b/classes.h
index a960f634..8f09c7b8 100644
--- a/classes.h
+++ b/classes.h
@@ -209,7 +209,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, 
+  gMacbeath, gBring, gSchmutzM2, gSchmutzM3, gCrystal,
   gGUARD};
 
 enum eGeometryClass { gcHyperbolic, gcEuclid, gcSphere };
diff --git a/compileunits.h b/compileunits.h
index 7ce4dbb0..030def83 100644
--- a/compileunits.h
+++ b/compileunits.h
@@ -31,6 +31,7 @@
 #include "heptagon.cpp"
 #include "binary-tiling.cpp"
 #include "archimedean.cpp"
+#include "crystal.cpp"
 #include "language.cpp"
 #include "cell.cpp"
 #include "expansion.cpp"
diff --git a/crystal.cpp b/crystal.cpp
new file mode 100644
index 00000000..07a7f1fe
--- /dev/null
+++ b/crystal.cpp
@@ -0,0 +1,457 @@
+// Hyperbolic Rogue
+// This file implements the multi-dimensional (aka crystal) geometries.
+// Copyright (C) 2011-2018 Zeno Rogue, see 'hyper.cpp' for details
+
+namespace hr {
+
+namespace crystal {
+
+bool add_bitruncation = false;
+bool view_coordinates = false;
+
+const int MAXDIM = 7;
+
+typedef array<int, MAXDIM> coord;
+static const coord c0 = {};
+
+int tocode(int cname) { return (1 << (cname >> 1)); }
+
+void resize2(vector<vector<int>>& v, int a, int b, int z) {
+  v.clear();
+  v.resize(a);
+  for(auto& w: v) w.resize(b, z);
+  }
+
+const int FULLSTEP = 16;
+const int HALFSTEP = 8;
+
+struct crystal_structure {
+  int dir;
+  int dim;
+  
+  vector<vector<int>> cmap;
+  vector<vector<int>> next;
+  vector<vector<int>> prev;
+  vector<vector<int>> order;
+  
+  void coord_to_next() {
+    resize2(next, 1<<dim, 2*dim, -1);
+    for(int a=0; a<(1<<dim); a++) 
+      for(int b=0; b<dir; b++)
+        next[a][cmap[a][b]] = cmap[a][(b+1)%dir];
+    println(hlog, next);
+    }
+  
+  void next_to_coord() {
+    resize2(cmap, 1<<dim, dir, -1);
+    for(int a=0; a<(1<<dim); a++) {
+      int at = 0;
+      for(int b=0; b<dir; b++) {
+        cmap[a][b] = at;
+        at = next[a][at];
+        }
+      }
+    println(hlog, "coordinate map is:\n", cmap);
+    }
+
+  void next_to_prev() {
+    resize2(prev, 1<<dim, 2*dim, -1);
+    for(int a=0; a<(1<<dim); a++)
+      for(int b=0; b<dir; b++) {
+        if(next[a][b] != -1)
+          prev[a][next[a][b]] = b;
+        }
+    }
+
+  void coord_to_order() {
+    println(hlog, dir, dim);
+    resize2(order, 1<<dim, 2*dim, -1);
+    for(int a=0; a<(1<<dim); a++) 
+    for(int b=0; b<dir; b++)
+      order[a][cmap[a][b]] = b;
+    println(hlog, order);
+    }
+  
+  int count_bugs() {
+    int bugcount = 0;
+      
+    for(int a=0; a<(1<<dim); a++) 
+    for(int b=0; b<2*dim; b++) {
+      if(next[a][b] == -1) continue;
+      int qa = a, qb = b;
+      for(int i=0; i<4; i++) {
+        if(i == 2 && (qb != (b^1))) bugcount++;
+        qa ^= tocode(qb);
+        qb ^= 1;
+        qb = next[qa][qb];
+        }
+      if(a != qa || b != qb) bugcount++;
+      }
+    
+    return bugcount;
+    }
+  
+  void next_insert(int a, int at, int val) {
+    int pd = next[a].size();
+
+    next[a].resize(pd + 2);
+    next[a][val] = next[a][at];
+    next[a][at] = val;
+    next[a][val^1] = next[a][at^1];
+    next[a][at^1] = val^1;
+
+    prev[a].resize(pd + 2);
+    prev[a][val] = at;
+    prev[a][next[a][val]] = val;
+    prev[a][val^1] = at^1;
+    prev[a][next[a][val^1]] = val^1;
+    }
+
+  void prev_insert(int a, int at, int val) {
+    next_insert(a, prev[a][at], val);
+    }
+  
+  int errors = 0;
+  
+  bool may_next_insert(int a, int at, int val) {
+    if(isize(next[a]) != dir) {
+      next_insert(a, at, val);
+      return true;
+      }
+    else if(next[a][at] != val) errors++;
+    return false;
+    }
+  
+  bool may_prev_insert(int a, int at, int val) {
+    if(isize(prev[a]) != dir) {
+      prev_insert(a, at, val);
+      return true;
+      }
+    else if(prev[a][at] != val) errors++;
+    return false;
+    }
+  
+  void add_dimension_to(crystal_structure& poor) {
+    dir = poor.dir + 2;
+    dim = poor.dim + 1;
+
+    printf("Building dimension %d\n", dim);
+
+    next.resize(1<<dim);
+    prev.resize(1<<dim);
+    int mask = (1<<poor.dim) - 1;
+    
+    int mm = tocode(poor.dir);
+
+    for(int i=0; i<(1<<dim); i++) {
+      if(i < mm)
+        next[i] = poor.next[i&mask], prev[i] = poor.prev[i&mask];
+      else
+        next[i] = poor.prev[i&mask], prev[i] = poor.next[i&mask];
+      }
+
+    next_insert(0, 0, poor.dir);
+    
+    for(int s=2; s<1<<(dim-2); s+=2) {
+      if(next[s][0] < 4)
+        prev_insert(s, 0, poor.dir);
+      else
+        next_insert(s, 0, poor.dir);
+      }
+    // printf("next[%d][%d] = %d\n", 4, 2, next[4][2]);
+    
+    for(int s=0; s<8; s++) for(int a=0; a<(1<<dim); a++) if(isize(next[a]) > poor.dir) {
+      int which = next[a][poor.dir];
+      int a1 = a ^ tocode(which);
+      
+      may_next_insert(a1, which^1, poor.dir);
+      may_next_insert(a ^ mm, which, poor.dir^1);
+      which = prev[a][poor.dir];
+      a1 = a ^ tocode(which);
+      may_prev_insert(a1, which^1, poor.dir);
+      }
+      
+    // println(hlog, next);
+    
+    if(errors) { printf("errors: %d\n", errors); exit(1);; }
+    
+    int unf = 0;
+    for(int a=0; a<(1<<dim); a++) if(isize(next[a]) == poor.dir) {
+      if(!unf) printf("unf: ");
+      printf("%d ", a);
+      unf ++;
+      }
+    
+    if(unf) { printf("\n"); exit(2); }
+
+    for(int a=0; a<(1<<dim); a++) for(int b=0; b<dir; b++)
+      if(prev[a][next[a][b]] != b) {
+        println(hlog, next[a], prev[a]);
+        printf("next/prev %d\n", a);
+        exit(3);
+        }
+      
+    if(count_bugs()) {
+      printf("bugs reported: %d\n", count_bugs());
+      exit(4);
+      }
+    }
+  
+  void remove_half_dimension() {
+    dir--;
+    for(int i=0; i<(1<<dim); i++) {
+      int take_what = dir;
+      if(i >= (1<<(dim-1))) take_what = dir-1;
+      next[i][prev[i][take_what]] = next[i][take_what],
+      prev[i][next[i][take_what]] = prev[i][take_what],
+      next[i].resize(dir),
+      prev[i].resize(dir);
+      }
+    }
+
+  void build() {
+    dir = 4;
+    dim = 2;
+    next.resize(4, {2,3,1,0});
+    next_to_prev();
+    while(dir < S7) {
+      crystal_structure csx = move(*this);
+      add_dimension_to(csx);
+      }
+    if(dir > S7) remove_half_dimension();
+      
+    next_to_coord();
+    
+    coord_to_order();
+    coord_to_next();
+    if(count_bugs()) {
+      printf("bugs found\n");
+      }
+    if(dir > MAX_EDGE || dim > MAXDIM) {
+      printf("Dimension or directions exceeded -- I have generated it, but won't play");
+      exit(0);
+      }
+    }
+  
+  };
+
+struct lwalker {
+  crystal_structure cs;
+  int id;
+  int spin;
+  };
+
+lwalker operator +(lwalker a, int v) { a.spin = gmod(a.spin + v, a.cs.dir); return a; }
+lwalker operator +(lwalker a, wstep_t) {
+  a.spin = a.cs.cmap[a.id][a.spin];
+  a.id ^= tocode(a.spin);
+  a.spin = a.cs.order[a.id][a.spin^1];
+  return a;
+  }
+
+coord add(coord c, lwalker a, int val) {
+  int code = a.cs.cmap[a.id][a.spin];
+  c[code>>1] += (code&1) ? val : -val;
+  return c;
+  }
+  
+map<heptagon*, coord> hcoords;
+map<coord, heptagon*> heptagon_at;
+
+crystal_structure cs;
+
+coord add(coord c, int cname, int val) {
+  int dim = (cname>>1);
+  c[dim] = (c[dim] + (cname&1?val:-val));
+  return c;
+  }
+
+lwalker makewalker(crystal_structure& cs, coord c, int d) {
+  lwalker a;
+  a.cs = cs;
+  a.id = 0;
+  for(int i=0; i<cs.dim; i++) if(c[i] & FULLSTEP) a.id += (1<<i);
+  a.spin = d;
+  return a;
+  }
+
+void crystalstep(heptagon *h, int d);
+
+heptagon *get_heptagon_at(coord c, int deg) {
+  if(heptagon_at.count(c)) return heptagon_at[c];
+  heptagon*& h = heptagon_at[c];
+  h = tailored_alloc<heptagon> (deg);
+  h->alt = NULL;
+  h->cdata = NULL;
+  h->c7 = newCell(deg, h);
+  h->distance = 0;
+  for(int i=0; i<cs.dim; i++) h->distance += abs(c[i]);
+  hcoords[h] = c;
+  // for(int i=0; i<6; i++) crystalstep(h, i);
+  return h;
+  }
+
+coord get_coord(cell *c) {
+  if(c->master->c7 != c) {
+    coord res = c0;
+    for(int i=0; i<c->type; i+=2) {
+      coord co = hcoords[c->move(i)->master];
+      for(int d=0; d<cs.dim; d++) res[d] += co[d];
+      }
+    for(int d=0; d<cs.dim; d++) res[d] = (2 * res[d] + c->type/2) / c->type;
+    return res;
+    }
+  else
+    return hcoords[c->master];
+  }
+
+struct hrmap_crystal : hrmap {
+  heptagon *getOrigin() { return get_heptagon_at(c0, S7); }
+
+  hrmap_crystal() {
+    cs.build();
+    }
+
+  void verify() { }
+  };
+ 
+hrmap *new_map() {
+  return new hrmap_crystal;
+  }
+
+bool is_bi(coord co) {
+  for(int i=0; i<cs.dim; i++) if(co[i] & HALFSTEP) return true;
+  return false;
+  }
+
+void create_step(heptagon *h, int d) {
+  if(geometry != gCrystal) return;
+  if(!hcoords.count(h)) {
+    printf("not found\n");
+    return;
+    }
+  auto co = hcoords[h];
+  
+  if(is_bi(co)) {
+    heptspin hs(h, d);
+    (hs + 1 + wstep + 1).cpeek();
+    return;
+    }
+
+  auto lw = makewalker(cs, co, d);
+
+  if(!add_bitruncation) {
+    auto c1 = add(co, lw, FULLSTEP);
+    auto lw1 = lw+wstep;
+    
+    h->c.connect(d, heptspin(get_heptagon_at(c1, S7), lw1.spin));
+    }
+  else {
+    auto coc = add(add(co, lw, HALFSTEP), lw+1, HALFSTEP);
+    auto hc = get_heptagon_at(coc, HALFSTEP);
+    for(int a=0; a<8; a+=2) {
+      hc->c.connect(a, heptspin(h, lw.spin));
+      if(h->modmove(lw.spin-1)) {
+        hc->c.connect(a+1, heptspin(h, lw.spin) - 1 + wstep - 1);
+        }
+      co = add(co, lw, FULLSTEP);
+      lw = lw + wstep + (-1);
+      h = get_heptagon_at(co, S7);
+      }
+    }
+  }
+
+array<array<int,2>, MAX_EDGE> distlimit_table = {{
+  {SEE_ALL,SEE_ALL}, {SEE_ALL,SEE_ALL}, {SEE_ALL,SEE_ALL}, {SEE_ALL,SEE_ALL}, {15, 10}, 
+  {6, 4}, {5, 3}, {4, 3}, {4, 3}, {3, 2}, {3, 2}, {3, 2}, {3, 2}, {3, 2}
+  }};
+
+int readArgs() {
+  using namespace arg;
+           
+  if(0) ;
+  else if(argis("-crystal")) {
+    stop_game();
+    geometry = gCrystal; variation = eVariation::pure;
+    shift(); int N = argi();
+    ginf[gCrystal].sides = N;
+    if(N < MAX_EDGE)
+      ginf[gCrystal].distlimit = distlimit_table[N];
+    add_bitruncation = false;
+    }
+  else if(argis("-crystalb")) {
+    stop_game();
+    geometry = gCrystal; variation = eVariation::bitruncated;
+    ginf[gCrystal].sides = 8;
+    add_bitruncation = true;
+    }
+  else if(argis("-cview")) {
+    view_coordinates = true;
+    }
+  else return 1;
+  return 0;
+  }
+
+color_t colorize(cell *c) {
+  coord co = get_coord(c);
+  color_t res;
+  res = 0;
+  for(int i=0; i<3; i++)
+    res |= ((i == 2 && S7 == 5) ? (co[i] ? 255 : 0) : (128 + co[i] * 3)) << (8*i);
+  return res;
+  }
+
+bool crystal_cell(cell *c, transmatrix V) {
+
+  if(geometry != gCrystal) return false;
+
+  if(view_coordinates && cheater) for(int i=0; i<S7; i++) {
+    
+    if(c->master->c7 == c) {    
+      transmatrix V1 = cellrelmatrix(c, i);
+      ld dist = hdist0(V1 * C0);
+      ld alpha = -atan2(V1 * C0);
+      transmatrix T = V * spin(alpha) * xpush(dist*.3);
+
+      auto co = hcoords[c->master];
+      int our_id = 0;
+      for(int a=0; a<MAXDIM; a++) if(co[a] & FULLSTEP) our_id += (1<<a);
+      int cx = cs.cmap[our_id][i];
+      
+      int coordcolors[MAXDIM] = {0x4040D0, 0x40D040, 0xD04040, 0xFFD500, 0xF000F0, 0x00F0F0, 0xF0F0F0 };
+    
+      queuestr(T, 0.3, its(co[cx>>1] / (add_bitruncation ? HALFSTEP : FULLSTEP)), coordcolors[cx>>1], 1);
+      }
+
+    if(PURE) {
+      cellwalker cw(c, i);
+      cellwalker cw2 = cw;
+      for(int i=0; i<(add_bitruncation?3:4); i++) cw2 = cw2 + wstep + 1;
+      if(cw2 != cw) { printf("crystal valence error\n"); cw.at->item = itGold; }
+      }
+    }
+  return false;
+  }
+
+int hypot2(coord co1, coord co2) {
+  int result = 0;
+  for(int a=0; a<cs.dim; a++) result += (co1[a] - co2[a]) * (co1[a] - co2[a]);
+  return result;
+  }
+  
+int distance(cell *c1, cell *c2) {
+  if(true || (PURE && !add_bitruncation)) {
+    coord co1 = hcoords[c1->master];
+    coord co2 = hcoords[c2->master];
+    int result = 0;
+    for(int a=0; a<cs.dim; a++) result += abs(co1[a] - co2[a]);
+    return result / FULLSTEP;
+    }
+  }
+
+auto crystalhook = addHook(hooks_args, 100, readArgs)
+  + addHook(hooks_drawcell, 100, crystal_cell);
+
+}
+
+}
diff --git a/geometry2.cpp b/geometry2.cpp
index a0146c04..651a91e0 100644
--- a/geometry2.cpp
+++ b/geometry2.cpp
@@ -140,10 +140,10 @@ transmatrix calc_relative_matrix(cell *c2, cell *c1, const hyperpoint& point_hin
       int sp = h2->c.spin(d);
       return gm * heptmove[sp] * spin(2*M_PI*d/S7) * where;
       }
-    if(among(geometry, gFieldQuotient, gBring, gMacbeath)) {
+    if(among(geometry, gFieldQuotient, gBring, gMacbeath, gCrystal)) {
       int bestdist = 1000, bestd = 0;
       for(int d=0; d<S7; d++) {
-        int dist = celldistance(h2->move(d)->c7, c1);
+        int dist = celldistance(h2->cmove(d)->c7, c1);
         if(dist < bestdist) bestdist = dist, bestd = d;
         }
       int sp = h2->c.spin(bestd);
diff --git a/heptagon.cpp b/heptagon.cpp
index 60f92c43..20395f7c 100644
--- a/heptagon.cpp
+++ b/heptagon.cpp
@@ -205,13 +205,15 @@ heptagon *createStep(heptagon *h, int d) {
   d = h->c.fix(d);
   if(!h->move(d))
     callhooks(hooks_createStep, h, d);
+  if(!h->move(d) && geometry == gCrystal) 
+    crystal::create_step(h, d);
   if(!h->move(d) && binarytiling) 
     return binary::createStep(h, d);
   if(!h->move(d) && archimedean) {
     arcm::create_adjacent(h, d); 
     return h->move(d);
     }
-  if(!h->move(0) && h->s != hsOrigin && !binarytiling) {
+  if(!h->move(0) && h->s != hsOrigin && !binarytiling && (geometry != gCrystal)) {
     // cheating: 
     int pard=0;
     if(S3 == 3) 
diff --git a/hyper.h b/hyper.h
index efebe819..a50820f0 100644
--- a/hyper.h
+++ b/hyper.h
@@ -4135,6 +4135,13 @@ namespace arcm {
   int fix(heptagon *h, int spin);
   }
 
+namespace crystal {
+  color_t colorize(cell *c);
+  int distance(cell *c1, cell *c2);
+  hrmap *new_map();
+  void create_step(heptagon *h, int d);
+  }
+
 hyperpoint get_warp_corner(cell *c, int cid);
 hyperpoint get_corner_position(cell *c, int cid, ld cf = 3);
 
diff --git a/pattern2.cpp b/pattern2.cpp
index 0a830f92..8e835c5c 100644
--- a/pattern2.cpp
+++ b/pattern2.cpp
@@ -1343,6 +1343,8 @@ namespace patterns {
         int z = currfp.getdist(fieldval(c), make_pair(0,false));
         return 255 * (currfp.maxdist+1-z) / currfp.maxdist;
         }
+      case 'K':
+        return crystal::colorize(c);
       case 'N': {
         if(!hyperbolic) return canvasback;
         using namespace fieldpattern;
@@ -1446,6 +1448,9 @@ namespace patterns {
     if(archimedean)
       dialog::addSelItem(XLAT("Archimedean"), "Archimedean", 'A');
 
+    if(geometry == gCrystal)
+      dialog::addSelItem(XLAT("Crystal coordinates"), "Crystal", 'K');
+
     dialog::addSelItem(XLAT("sides"), "sides", 'B');
 
     dialog::addBreak(100);