From 59a9605b1312c728c6c60f543ae1c883c26f91d5 Mon Sep 17 00:00:00 2001
From: Zeno Rogue <zeno@attnam.com>
Date: Wed, 15 Jan 2020 17:58:41 +0100
Subject: [PATCH] arb::relative_matrix

---
 arbitrile.cpp    |  29 +++++
 cell.cpp         |   2 +-
 fieldpattern.cpp | 256 +++++++++++++++++++++++++++++++++++++-
 geom-exp.cpp     |  24 +++-
 reg3.cpp         | 316 +++++++++++++++++++----------------------------
 5 files changed, 430 insertions(+), 197 deletions(-)

diff --git a/arbitrile.cpp b/arbitrile.cpp
index e7c24cda..8919e239 100644
--- a/arbitrile.cpp
+++ b/arbitrile.cpp
@@ -443,6 +443,35 @@ struct hrmap_arbi : hrmap {
       }
     }
 
+  transmatrix relative_matrix(heptagon *h2, heptagon *h1, const hyperpoint& hint) override {
+    if(gmatrix0.count(h2->c7) && gmatrix0.count(h1->c7))
+      return inverse(gmatrix0[h1->c7]) * gmatrix0[h2->c7];
+    transmatrix gm = Id, where = Id;
+    while(h1 != h2) {
+      for(int i=0; i<h1->type; i++) {
+        if(h1->move(i) == h2) {
+          return gm * adj(h1, i) * where;
+          }
+        }
+      if(h1->distance > h2->distance) {
+        for(int i=0; i<h1->type; i++) if(h1->move(i) && h1->move(i)->distance < h1->distance) {
+          gm = gm * adj(h1, i);
+          h1 = h1->move(i);
+          goto again;
+          }
+        }
+      else {
+        for(int i=0; i<h2->type; i++) if(h2->move(i) && h2->move(i)->distance < h2->distance) {
+          where = iadj(h2, 0) * where;
+          h2 = h2->move(i);
+          goto again;
+          }
+        }
+      again: ;
+      }
+    return gm * where;
+    }
+
   transmatrix adj(cell *c, int dir) override { return adj(c->master, dir); }
   
   ld spin_angle(cell *c, int d) override { return SPIN_NOT_AVAILABLE; }
diff --git a/cell.cpp b/cell.cpp
index a76ff566..5ee179a4 100644
--- a/cell.cpp
+++ b/cell.cpp
@@ -983,7 +983,7 @@ EX int celldistance(cell *c1, cell *c2) {
   if(hybri) return hybrid::celldistance(c1, c2);
   
   #if CAP_FIELD
-  if(geometry == gFieldQuotient && !GOLDBERG)
+  if(geometry == gFieldQuotient && !GOLDBERG && WDIM == 2)
     return currfp.getdist(fieldpattern::fieldval(c1), fieldpattern::fieldval(c2));
   #endif
   
diff --git a/fieldpattern.cpp b/fieldpattern.cpp
index 7e5f195b..27ea5c7c 100644
--- a/fieldpattern.cpp
+++ b/fieldpattern.cpp
@@ -29,7 +29,7 @@ struct fgeomextra {
   };
 #endif
 
-bool isprime(int n) {
+EX bool isprime(int n) {
   for(int k=2; k<n; k++) if(n%k == 0) return false;
   return true;
   }
@@ -69,6 +69,8 @@ EX int btspin(int id, int d) {
 
 #if HDR
 struct fpattern {
+
+  unsigned force_hash;
   
   int Prime, wsquare, Field, dual;
   // we perform our computations in the field Z_Prime[w] where w^2 equals wsquare
@@ -273,6 +275,7 @@ struct fpattern {
     }
     
   fpattern(int p) {
+    force_hash = 0;
     if(!p) return;
     init(p);
     }
@@ -282,6 +285,17 @@ struct fpattern {
   vector<matrix> generate_isometries();
   
   bool check_order(matrix M, int req);
+  
+  unsigned compute_hash();
+
+  // general 4D
+  vector<transmatrix> fullv;
+
+  void add1(const matrix& M);
+  void add1(const matrix& M, const transmatrix& Full);
+  vector<matrix> generate_isometries3();
+  int solve3();
+  void generate_all3();
   };
 #endif
 
@@ -323,6 +337,166 @@ vector<matrix> fpattern::generate_isometries() {
   return res;
   }
 
+vector<matrix> fpattern::generate_isometries3() {
+  
+  matrix T = Id;
+  int low = wsquare ? 1-Prime : 0;
+  vector<matrix> res;
+  
+  auto colprod = [&] (int a, int b) {
+    return add(add(mul(T[0][a], T[0][b]), mul(T[1][a], T[1][b])), sub(mul(T[2][a], T[2][b]), mul(T[3][a], T[3][b])));
+    };
+
+  auto rowcol = [&] (int a, int b) {
+    return add(add(mul(T[a][0], T[0][b]), mul(T[a][1], T[1][b])), add(mul(T[a][2], T[2][b]), mul(T[a][3], T[3][b])));
+    };
+
+  for(T[0][0]=low; T[0][0]<Prime; T[0][0]++)
+  for(T[1][0]=low; T[1][0]<Prime; T[1][0]++)
+  for(T[2][0]=low; T[2][0]<Prime; T[2][0]++)
+  for(T[3][0]=low; T[3][0]<Prime; T[3][0]++)
+  if(colprod(0, 0) == 1)
+  for(T[0][1]=low; T[0][1]<Prime; T[0][1]++)
+  for(T[1][1]=low; T[1][1]<Prime; T[1][1]++)
+  for(T[2][1]=low; T[2][1]<Prime; T[2][1]++)
+  for(T[3][1]=low; T[3][1]<Prime; T[3][1]++)
+  if(colprod(1, 1) == 1)
+  if(colprod(1, 0) == 0)
+  for(T[0][2]=low; T[0][2]<Prime; T[0][2]++)
+  for(T[0][3]=low; T[0][3]<Prime; T[0][3]++)
+  if(rowcol(0, 0) == 1)
+  if(rowcol(0, 1) == 0)
+  for(T[1][2]=low; T[1][2]<Prime; T[1][2]++)
+  for(T[1][3]=low; T[1][3]<Prime; T[1][3]++)
+  if(rowcol(1, 0) == 0)
+  if(rowcol(1, 1) == 1)
+  for(T[2][2]=low; T[2][2]<Prime; T[2][2]++)
+  for(T[3][2]=low; T[3][2]<Prime; T[3][2]++)
+  if(colprod(2, 2) == 1)
+  if(colprod(2, 0) == 0)
+  if(colprod(2, 1) == 0)
+  for(T[2][3]=low; T[2][3]<Prime; T[2][3]++)
+  for(T[3][3]=low; T[3][3]<Prime; T[3][3]++)
+  if(rowcol(2, 0) == 0)
+  if(rowcol(2, 1) == 0)
+  if(rowcol(2, 2) == 1)
+  // if(colprod(3, 3) == 1)
+  if(add(colprod(3, 3), 1) == 0)
+  if(colprod(3, 0) == 0)
+  if(colprod(3, 1) == 0)
+  if(colprod(3, 2) == 0)
+  if(rowcol(3, 3) == 1)
+  if(rowcol(3, 0) == 0)
+  if(rowcol(3, 1) == 0)
+  if(rowcol(3, 2) == 0)
+    res.push_back(T);
+
+  return res;
+  }
+
+void fpattern::add1(const matrix& M) {
+  if(!matcode.count(M)) {
+    int i = matrices.size();
+    matcode[M] = i, matrices.push_back(M);
+    }
+  }
+
+void fpattern::add1(const matrix& M, const transmatrix& Full) {
+  if(!matcode.count(M)) {
+    int i = matrices.size();
+    matcode[M] = i, matrices.push_back(M), fullv.push_back(Full);
+    }
+  }
+
+map<unsigned,int> hash_found;
+
+unsigned fpattern::compute_hash() {
+  unsigned hashv = 0;
+  int iR = matcode[R];
+  int iP = matcode[P];
+  int iX = matcode[X];
+  for(int i=0; i<isize(matrices); i++) {
+    hashv = 3 * hashv + gmul(i, iP) + 7 * gmul(i, iR);
+    if(MWDIM == 4) hashv += 11 * gmul(i, iX);
+    }
+  return hashv;
+  }
+
+void fpattern::generate_all3() {
+  matrices.clear();
+  matcode.clear();
+  add1(Id);
+  fullv = {hr::Id};
+  for(int i=0; i<isize(matrices); i++) {
+    add1(mmul(matrices[i], R), fullv[i] * reg3::full_R);
+    add1(mmul(matrices[i], X), fullv[i] * reg3::full_X);
+    }
+  local_group = isize(matrices);
+  for(int i=0; i<(int)matrices.size(); i++) {
+    matrix E = mmul(matrices[i], P);
+    if(!matcode.count(E))
+      for(int j=0; j<local_group; j++) add1(mmul(E, matrices[j]));
+    }
+  unsigned hashv = compute_hash();
+  DEBB(DF_FIELD, ("all = ", isize(matrices), "/", local_group, " = ", isize(matrices) / local_group, " hash = ", hashv, " count = ", ++hash_found[hashv]));
+  }
+
+int fpattern::solve3() {
+  if(!wsquare) return 0;
+  reg3::construct_relations();
+  
+  DEBB(DF_FIELD, ("generating isometries for ", Field));
+  
+  auto iso3 = generate_isometries();
+  auto iso4 = generate_isometries3();
+
+  int cmb = 0;
+  
+  int N = isize(reg3::rels);
+  
+  vector<int> fails(N);
+  
+  vector<matrix> possible_P, possible_X, possible_R;
+  
+  for(auto& M: iso3) {
+    if(check_order(M, 2)) 
+      possible_X.push_back(M);
+    if(check_order(M, reg3::r_order)) 
+      possible_R.push_back(M);
+    }
+  for(auto& M: iso4) 
+    if(check_order(M, 2)) 
+      possible_P.push_back(M);
+    
+  DEBB(DF_FIELD, ("field = ", Field, " #P = ", isize(possible_P), " #X = ", isize(possible_X), " #R = ", isize(possible_R), " r_order = ", reg3::r_order, " xp_order = ", reg3::xp_order));
+                                                                                                                               
+  for(auto& xX: possible_X)
+  for(auto& xP: possible_P) if(check_order(mmul(xP, xX), reg3::xp_order))
+  for(auto& xR: possible_R) if(check_order(mmul(xR, xX), reg3::rx_order)) { // if(xR[0][0] == 1 && xR[0][1] == 0) {
+    auto by = [&] (char ch) -> matrix& { return ch == 'X' ? xX : ch == 'R' ? xR : xP; };
+    for(int i=0; i<N; i++) {
+      matrix ml = Id;
+      for(char c: reg3::rels[i].first) { ml = mmul(ml, by(c)); if(ml == Id) { fails[i]++; goto bad; }}
+      matrix mr = Id;
+      for(char c: reg3::rels[i].second) { mr = mmul(mr, by(c)); if(mr == Id) { fails[i]++; goto bad; }}
+      if(ml != mr) { fails[i]++; goto bad;}
+      }
+    P = xP; R = xR; X = xX;
+    generate_all3();
+    if(force_hash && compute_hash() != force_hash) continue;
+    cmb++;
+    goto ok;
+    bad: ;
+    }
+  
+  ok:
+
+  DEBB(DF_FIELD, ("cmb = ", cmb, " for field = ", Field));
+  for(int i=0; i<N; i++) if(fails[i]) DEBB(DF_FIELD, (reg3::rels[i], " fails = ", fails[i]));
+  
+  return cmb;
+  }
+
 int fpattern::solve() {
   
   for(int a=0; a<MWDIM; a++) for(int b=0; b<MWDIM; b++) Id[a][b] = a==b?1:0;
@@ -347,6 +521,14 @@ int fpattern::solve() {
         }
       } else wsquare = 0;
 
+    if(WDIM == 3) {
+      if(dual == 0) {
+        int s = solve3();
+        if(s) return 0;
+        }
+      continue;
+      }
+
     if(dual == 2) {
       if(Field <= 10) {
         vector<matrix> all_isometries = generate_isometries();
@@ -419,6 +601,8 @@ int fpattern::order(const matrix& M) {
 
 void fpattern::build() {
 
+  if(WDIM == 3) return;
+  
   for(int i=0; i<isize(qpaths); i++) {
     matrix M = strtomatrix(qpaths[i]);
     qcoords.push_back(M);
@@ -789,6 +973,7 @@ EX void info() {
       #ifndef EASY
       neasy = 0; 
       #endif
+      if(WDIM == 3) continue;
       fp.build();
       #ifndef EASY
       printf("Not easy: %d\n", neasy);
@@ -812,12 +997,61 @@ EX bool quotient_field_changed;
 
 EX struct fpattern& getcurrfp() {
   if(geometry == gFieldQuotient && quotient_field_changed)
-    return current_quotient_field;  
-  if(WDIM == 3) {
-    dynamicval<eGeometry> g(geometry, gSpace435);
+    return current_quotient_field;
+  if(geometry == gSpace535) {
+    // 120 cells, hash = 9EF7A9C4
     static fpattern fp(5);
     return fp;
     }
+  if(geometry == gSpace534) {
+    // 260 cells, hash = 72414D0C (not 0C62E214)
+    static fpattern fp(0);
+    if(fp.Prime) return fp;
+    fp.Prime = 5; fp.force_hash = 0x72414D0C; fp.solve();
+    return fp;
+    }
+  if(geometry == gSpace435) {
+    // 650 cells, hash = EB201050
+    static fpattern fp(0);
+    if(fp.Prime) return fp;
+    fp.Prime = 5; fp.force_hash = 0x72414D0C; fp.solve();
+    return fp;
+    }
+  if(geometry == gSpace336) {
+    // 672 cells in E3F6B7BC
+    // 672 cells in 885F1184
+    // 9408 cells in C4089F34
+    static fpattern fp(0);
+    if(fp.Prime) return fp;
+    fp.Prime = 7; fp.force_hash = 0xE3F6B7BCu; fp.solve();
+    return fp;
+    }
+  if(geometry == gSpace344) {
+    // 600 cells in 558C8ED0
+    // 2400 cells in AF042EA8 
+    // 2600 cells in D26948E0 
+    // 2600 cells in EC29DCEC 
+    static fpattern fp(0);
+    if(fp.Prime) return fp;
+    fp.Prime = 5; fp.force_hash = 0x558C8ED0u; fp.solve();
+    return fp;
+    // 4900 cells in CDCC7860 (7)
+    }
+  if(geometry == gSpace536) {
+    static fpattern fp(0);
+    if(fp.Prime) return fp;
+    // 130 cells in 3BA5C5A4
+    // 260 cells in 9FDE7B38
+    fp.Prime = 7; fp.force_hash = 0x9FDE7B38u; fp.solve();
+    return fp;
+    }
+  if(WDIM == 3) {
+    static fpattern fp(0);
+    if(fp.Prime) return fp;
+    for(int p=2; p<8; p++) { fp.Prime = p; if(!fp.solve()) break; }
+    DEBB(DF_FIELD, ("set prime = ", fp.Prime));
+    return fp;
+    }
   if(S7 == 8 && S3 == 3) {
     static fpattern fp(17);
     return fp;
@@ -900,6 +1134,20 @@ EX void enableFieldChange() {
   fieldpattern::current_quotient_field.init(gxcur.primes[gxcur.current_prime_id].p);
   }
 
+EX void field_from_current() {
+  auto& go = ginf[geometry];
+  dynamicval<eGeometry> g(geometry, gFieldQuotient);
+  auto& gg = ginf[geometry];
+  gg.sides = go.sides;
+  gg.vertex = go.vertex;
+  gg.distlimit = go.distlimit;
+  gg.tiling_name = go.tiling_name;
+  gg.flags = go.flags | qANYQ | qFIELD | qBOUNDED;
+  gg.g = go.g;
+  gg.default_variation = go.default_variation;
+  fieldpattern::quotient_field_changed = true;
+  }
+
 EX }
 
 #define currfp fieldpattern::getcurrfp()
diff --git a/geom-exp.cpp b/geom-exp.cpp
index 9e75b21e..55d164da 100644
--- a/geom-exp.cpp
+++ b/geom-exp.cpp
@@ -316,6 +316,8 @@ void set_or_configure_geometry(eGeometry g) {
 bool same_tiling(eGeometry g2) {
   if(g2 == gCrystal)
     return S3 == 4;
+  if(g2 == gFieldQuotient && hyperbolic && standard_tiling())
+    return true;
   if(g2 == gFieldQuotient && geometry != gFieldQuotient) {
     int ce = 0;
     for(auto& ge: fieldpattern::fgeomextras) {
@@ -456,7 +458,15 @@ EX void select_quotient_screen() {
            "no quotient",
         g == geometry, key++);
       dialog::add_action([g] {
-        if(g == gFieldQuotient) 
+        if(g == gFieldQuotient && WDIM == 3) {
+          stop_game();
+          fieldpattern::field_from_current();
+          set_geometry(gFieldQuotient);
+          for(int p=2;; p++) { currfp.Prime = p; currfp.force_hash = 0; if(!currfp.solve()) break; }
+          println(hlog, "set prime = ", currfp.Prime);
+          start_game();
+          }
+        else if(g == gFieldQuotient) 
           pushScreen(showQuotientConfig);
         else if(g == gCrystal)
           pushScreen(crystal::show);
@@ -900,6 +910,14 @@ int read_geom_args() {
     enableFieldChange();
     set_geometry(gFieldQuotient);
     }
+  else if(argis("-to-fq")) {
+    shift(); unsigned hash = arghex();
+    stop_game_and_switch_mode(rg::nothing);
+    fieldpattern::field_from_current();
+    set_geometry(gFieldQuotient);
+    for(int p=2;; p++) { currfp.Prime = p; currfp.force_hash = hash; if(!currfp.solve()) break; }
+    println(hlog, "set prime = ", currfp.Prime);
+    }
   else if(argis("-cs")) {
     shift(); cheat();
     fieldpattern::matrix M = currfp.strtomatrix(args());
@@ -940,6 +958,10 @@ int read_geom_args() {
     fieldpattern::info();
     exit(0);
     } 
+  else if(argis("-fi-geo")) {
+    fieldpattern::info();
+    exit(0);
+    } 
   else if(argis("-qs")) {
     cheat();
     shift(); currfp.qpaths.push_back(args());
diff --git a/reg3.cpp b/reg3.cpp
index f9b2429c..87aeab3e 100644
--- a/reg3.cpp
+++ b/reg3.cpp
@@ -32,7 +32,7 @@ EX namespace reg3 {
   EX int face;
 
   EX vector<hyperpoint> cellshape;
-  vector<hyperpoint> vertices_only;
+  EX vector<hyperpoint> vertices_only;
   
   EX transmatrix spins[12], adjmoves[12];
 
@@ -237,6 +237,8 @@ EX namespace reg3 {
     
     void initialize(int cell_count);
     vector<cell*>& allcells() override { return acells; }
+
+    vector<hyperpoint> get_vertices(cell* c) override { return vertices_only; }
     };
   #endif
   
@@ -324,196 +326,37 @@ EX namespace reg3 {
         }
       }        
     };
-    
-  struct hrmap_field3 : hrmap_quotient3 {
-    
-    int mgmul(std::initializer_list<int> v) { 
-      int a = 0;
-      for(int b: v) a = a ? currfp_gmul(a, b) : b;
-      return a;
-      }
-    
-    vector<transmatrix> fullmatrices;
-    
-    int P, R, X;
-    transmatrix full_P, full_R, full_X;
-    
-    vector<int> field_adjmoves;
-    vector<int> cyclers;
-    int perm_group;
+
+  struct hrmap_field3 : reg3::hrmap_quotient3 {
   
-    vector<int> cell_to_code;
-    vector<int> code_to_cell;
-    
-    void seek(set<int>& seen_matrices, set<int>& seen_codes, const transmatrix& at, int ccode, const hyperpoint checker) {
-      if(hdist0(tC0(at)) > 4) return;
-      int b = bucketer(tC0(at));
-      if(seen_matrices.count(b)) return;
-      seen_matrices.insert(b);
-      for(int a=0; a<perm_group; a++) {
-        transmatrix T = at * fullmatrices[a];
-        if(hdist(T * checker, checker) < 1e-2) {
-          int co = mgmul({ccode, a});
-          seen_codes.insert(co);
-          fullmatrices[co] = T;
-          }
-        }
-      for(int a=0; a<perm_group; a++) seek(seen_matrices, seen_codes, at * fullmatrices[a] * full_P, mgmul({ccode, a, P}), checker);
-      }
-    
     hrmap_field3() {
-      eGeometry g = geometry;
-      geometry = gSpace435;
-      reg3::generate();
-      R = currfp_get_R();
-      P = currfp_get_P();
-      X = currfp_get_X();
-      full_P = reg3::adjmoves[0] * cspin(0, 2, M_PI) * cspin(0, 1, M_PI);
-      full_R = spin(-2 * M_PI / 4);
-      full_X = cspin(1, 2, M_PI / 2);
-  
-      DEBB(DF_FIELD, ("full_P = ", full_P, " / ", R));
-      DEBB(DF_FIELD, ("full_R = ", full_R, " / ", P));
-      DEBB(DF_FIELD, ("full_X = ", full_X, " / ", X));
-  
-      int N = currfp_n();
-  
-      perm_group = 24;
-      fullmatrices.resize(N);
-      fullmatrices[0] = Id;
-      vector<bool> known(perm_group, false);
-      known[0] = true;
-      for(int a=0; a<perm_group; a++) 
-      for(int i=0; i<perm_group; i++) if(known[i]) {
-        int iR = currfp_gmul(i, R);
-        fullmatrices[iR] = fullmatrices[i] * full_R;
-        known[iR] = true;
-        int iX = currfp_gmul(i, X);
-        fullmatrices[iX] = fullmatrices[i] * full_X;
-        known[iX] = true;
-        }
-      for(int i=0; i<perm_group; i++) if(known[i]) {
-        DEBB(DF_FIELD, (i, ". ", fullmatrices[i]));
-        }
+
+      auto& f = currfp;
+      auto lgr = f.local_group;
       
-      // find cav such that:
-      // cav * Id          * C0 = corner0
-      // cav * adjmoves[0] * C0 = corner1
-      // cav * adjmoves[1] * C0 = corner3
-      // cav * adjmoves[2] * C0 = cornerx
+      int N = isize(f.matrices) / lgr;
+      initialize(N);
       
-      hyperpoint corner0 = reg3::cellshape[0];
-      hyperpoint corner1 = reg3::cellshape[1];
-      hyperpoint corner3 = reg3::cellshape[3];
-      hyperpoint cornerx;
+      vector<int> moveid(S7), movedir(lgr);
+      for(int s=0; s<lgr; s++) 
+      for(int i=0; i<S7; i++) if(eqmatrix(f.fullv[s] * reg3::adjmoves[0], reg3::adjmoves[i]))
+        moveid[i] = s;
   
-      for(hyperpoint h: reg3::cellshape) DEBB(DF_FIELD, ("some corner ", h));
+      for(int s=0; s<lgr; s++) 
+      for(int i=0; i<S7; i++) if(hdist(tC0(inverse(f.fullv[s]) * reg3::adjmoves[0]), tC0(reg3::adjmoves[i])) < 1e-4)
+        movedir[s] = i;
   
-      for(hyperpoint h: reg3::cellshape)
-        if(hdist(h, corner1) > .1 && hdist(h, corner3) > .1 && abs(hdist(h, corner0)-hdist(corner0, corner1)) < .1)
-          cornerx = h;
-      DEBB(DF_FIELD, ("corner0 = ", corner0));
-      DEBB(DF_FIELD, ("corner1 = ", corner1));
-      DEBB(DF_FIELD, ("corner3 = ", corner3));
-      DEBB(DF_FIELD, ("cornerx = ", cornerx));
-      
-      transmatrix adj = Id, iadj = Id;
-  
-      geometry = g;
-      reg3::generate();
-      
-      cyclers.clear();
-      DEBB(DF_FIELD, ("S7 = ", S7));
-      if(S7 == 12) {
-      
-        transmatrix resmatrix;
-        set_column(resmatrix, 0, corner0);
-        set_column(resmatrix, 1, corner1);
-        set_column(resmatrix, 2, corner3);
-        set_column(resmatrix, 3, cornerx);
-        
-        transmatrix transformer;
-        set_column(transformer, 0, C0);
-        set_column(transformer, 1, tC0(reg3::adjmoves[0]));
-        set_column(transformer, 2, tC0(reg3::adjmoves[1]));
-        set_column(transformer, 3, tC0(reg3::adjmoves[2]));
-        
-        transmatrix cav = resmatrix * inverse(transformer);
-        DEBB(DF_FIELD, ("cav = ", cav));
-        DEBB(DF_FIELD, ("cav * C0 = ", cav * C0));
-  
-        set<int> seen_matrices;
-        set<int> seen_codes;
-        seek(seen_matrices, seen_codes, Id, 0, corner0);
-        
-        for(int x: seen_codes) cyclers.push_back(x);
-        perm_group = isize(cyclers);
-        adj = cav;
-        iadj = inverse(cav);
-        }
-      else {
-        for(int i=0; i<perm_group; i++) cyclers.push_back(i);
-        }
-      
-      field_adjmoves.resize(S7);
-      for(int i=0; i<S7; i++) field_adjmoves[i] = -1;
-  
-      for(int i=0; i<S7; i++) 
-        for(int a: cyclers)
-        for(int b: cyclers) {
-          transmatrix T = iadj * fullmatrices[a] * full_P * fullmatrices[b] * adj;
-          if(eqmatrix(T, reg3::adjmoves[i])) {
-            int code = mgmul({a,P,b});
-            field_adjmoves[i] = code;
-            DEBB(DF_FIELD, (i, " = ", make_tuple(a,P,b), " = ", code, " T = ", T));
+      for(int a=0; a<N; a++) {
+        tmatrices[a].resize(S7);
+        for(int b=0; b<S7; b++) {
+          int k = lgr*a;
+          k = f.gmul(f.gmul(k, moveid[b]), lgr);
+          for(int l=0; l<lgr; l++) if(f.gmul(k, l) % lgr == 0) {
+            tmatrices[a][b] = reg3::adjmoves[b] * f.fullv[l];
+            allh[a]->c.connect(b, allh[k/lgr], movedir[l], false);
             }
           }
-        
-      DEBB(DF_FIELD, ("field_adjmoves = ", field_adjmoves));
-      
-      DEBB(DF_FIELD, ("finding code_to_cell/cell_to_code..."));
-  
-      cell_to_code.clear();
-      code_to_cell.resize(N);
-      for(int i=0; i<N; i++) code_to_cell[i] = -1;
-      for(int i=0; i<N; i++) if(code_to_cell[i] == -1) {
-        for(int j: cyclers) code_to_cell[currfp_gmul(i, j)] = isize(cell_to_code);
-        cell_to_code.push_back(i);
-        }
-      
-      DEBB(DF_FIELD, ("building allh..."));
-      int cells = N / perm_group;
-      initialize(cells);
-      
-      DEBB(DF_FIELD, ("finding tmatrices..."));      
-      for(int i=0; i<cells; i++) {
-        for(int d=0; d<S7; d++) {
-          int found = 0;
-          int tmul = currfp_gmul(cell_to_code[i], field_adjmoves[d]);
-          for(int s: cyclers) {
-            int tmul2 = currfp_gmul(tmul, s);
-            if(cell_to_code[code_to_cell[tmul2]] == tmul2) {
-              allh[i]->move(d) = allh[code_to_cell[tmul2]];
-              allh[i]->c7->move(d) = allh[i]->move(d)->c7;
-              tmatrices[i].push_back(reg3::adjmoves[d] * iadj * fullmatrices[s] * adj);
-              found++;
-              }
-            }
-          if(found != 1) DEBB(DF_FIELD, ("bad found: ", i, "/", d, "/", found));
-          // println(hlog, "tmatrix(",i,",",d,") = ", tmatrices[i][d]);
-          }
-        }
-  
-      DEBB(DF_FIELD, ("setting spin..."));
-      for(int i=0; i<cells; i++) 
-      for(int d=0; d<S7; d++)
-      for(int e=0; e<S7; e++) 
-        if(allh[i]->move(d)->move(e) == allh[i]) {
-          allh[i]->c.setspin(d, e, false);
-          allh[i]->c7->c.setspin(d, e, false);
-          }
-  
-      DEBB(DF_FIELD, ("creating patterns..."));
+        }      
       create_patterns();
       }
     
@@ -529,8 +372,9 @@ EX namespace reg3 {
     
     
     void create_patterns() {
+      auto& f = currfp;
       // change the geometry to make sure that the correct celldistance is used
-      dynamicval<eGeometry> g(geometry, S7 == 12 ? gField534 : gField435);
+      dynamicval<eGeometry> g(geometry, gFieldQuotient);
       // also, strafe needs currentmap
       dynamicval<hrmap*> c(currentmap, this);
 
@@ -572,7 +416,7 @@ EX namespace reg3 {
         for(int i=0; i<currfp_n(); i++) {
           bool ok = true;
           for(auto o: plane_indices) {
-            int j = code_to_cell[currfp_gmul(i, cell_to_code[o])];
+            int j = currfp_gmul(i, o * f.local_group) / f.local_group;
             if(plane_indices.count(j)) ok = false;
             forCellEx(c1, allcells()[j]) if(plane_indices.count(c1->master->fieldval)) ok = false;
             }
@@ -595,17 +439,13 @@ EX namespace reg3 {
         int u = 0;
         for(int a=0; a<5; a++) {
           for(int o: plane_indices) {
-            int j = code_to_cell[currfp_gmul(u, cell_to_code[o])];
+            int j = currfp_gmul(u, o * f.local_group) / f.local_group;
             allcells()[j]->master->zebraval |= 2;
             }
           u = currfp_gmul(u, gpow);
           }
         }
       }
-    
-    vector<hyperpoint> get_vertices(cell* c) override {
-      return vertices_only;
-      }
     };
 
   /** homology cover of the Seifert-Weber space */
@@ -770,7 +610,7 @@ EX namespace reg3 {
         quotient_map = new seifert_weber::hrmap_seifert_cover;
         h.zebraval = quotient_map->allh[0]->zebraval;
         }
-      else if(hyperbolic && !(cgflags & qIDEAL)) {
+      else if(hyperbolic) {
         quotient_map = new hrmap_field3;
         h.zebraval = quotient_map->allh[0]->zebraval;
         }
@@ -1085,6 +925,7 @@ EX bool pseudohept(cell *c) {
     return hr::celldistance(c, currentmap->gamestart()) & 1;
   if(geometry == gCrystal344 || geometry == gCrystal534 || geometry == gSeifertCover)
     return false;
+  if(quotient) return false; /* added */
   if(hyperbolic) {
     heptagon *h = m->reg_gmatrix[c->master].first;
     return (h->zebraval == 1) && (h->distance & 1);
@@ -1225,6 +1066,99 @@ EX cellwalker strafe(cellwalker cw, int j) {
   println(hlog, "incorrect strafe");
   exit(1);
   }
+
+EX vector<pair<string, string> > rels;
+EX int xp_order, r_order, rx_order;
+
+EX transmatrix full_X, full_R, full_P;
+geometry_information *for_cgi;
+
+EX int matrix_order(const transmatrix A) {
+  transmatrix T = A;
+  int res = 1;
+  while(!eqmatrix(T, Id)) {
+    res++; T = T * A;
+    }
+  return res;
+  }
+
+EX void construct_relations() {
+  if(for_cgi == &cgi) return;
+  for_cgi = &cgi;
+  rels.clear();
+
+  reg3::generate();
+  reg3::generate_cellrotations();
+  vector<transmatrix> all;
+
+  vector<string> formulas;
+  
+  formulas.push_back("");
+
+  all.push_back(Id);
+  hyperpoint v = reg3::cellshape[0];
+  auto add = [&] (transmatrix T) {
+    for(int i=0; i<isize(all); i++) if(eqmatrix(all[i], T)) return i;
+    int S = isize(all);
+    all.push_back(T);
+    return S;
+    };
+  
+  auto cons = [&] (int i0, int i1, int i2) {
+    using reg3::adjmoves;
+    transmatrix T = build_matrix(adjmoves[ 0]*C0, adjmoves[ 1]*C0, adjmoves[ 2]*C0, C0);
+    transmatrix U = build_matrix(adjmoves[i0]*C0, adjmoves[i1]*C0, adjmoves[i2]*C0, C0);
+    return U * inverse(T);
+    };
+  
+  full_P = reg3::adjmoves[0];
+  full_R = S7 == 8 ? cons(1, 7, 0) : cons(1, 2, 0);
+  full_X = S7 == 8 ? cons(1, 0, 6) : S7 == 6 ? cons(1, 0, 5) : cons(1, 0, reg3::face);
+  
+  println(hlog, reg3::cellshape);
+
+  println(hlog, "cellshape = ", isize(reg3::cellshape));
+  bool ok = true;
+  int last_i = -1;
+  for(hyperpoint h: reg3::cellshape) {
+    int i = 0, j = 0;
+    for(hyperpoint u: reg3::cellshape) if(hdist(h, full_X*u) < 1e-4) i++;
+    for(hyperpoint u: reg3::cellshape) if(hdist(h, full_R*u) < 1e-4) j++;
+    if(last_i == -1) last_i = i;
+    if(i != j || i != last_i) ok = false;
+    }
+  
+  if(!ok) { println(hlog, "something wrong"); exit(1); }
+  
+  add(Id);
+  
+  auto work = [&] (transmatrix T, int p, char c) {
+    if(hdist0(tC0(T)) > 5) return;
+    for(hyperpoint h: reg3::cellshape) if(hdist(T * h, v) < 1e-4) goto ok;
+    return;
+    ok:
+    int id = add(T);
+    // println(hlog, p, " x ", (s0+c), " = ", id);
+
+    if(id >= isize(formulas)) formulas.push_back(formulas[p] + c);
+    else if(id == 0) println(hlog, "reached identity: ", formulas[p]+c);
+    else if(formulas[p][0] != formulas[id][0])
+      rels.emplace_back(formulas[p] + c, formulas[id]);
+    };
+  
+  for(int i=0; i<isize(all); i++) {
+    transmatrix T = all[i];
+    work(T * full_R, i, 'R');
+    work(T * full_X, i, 'X');
+    work(T * full_P, i, 'P');
+    }
+  
+  xp_order = matrix_order(full_X * full_P);
+  r_order = matrix_order(full_R);
+  rx_order = matrix_order(full_R * full_X);
+  println(hlog, "orders = ", tie(rx_order, r_order, xp_order));
+  }
+
 EX }
 #endif
 }