// Hyperbolic Rogue -- regular honeycombs
// works with spherical and hyperbolic ones -- Euclidean cubic tiling implemented in euclid.cpp
// hyperbolic honeycombs rely on binary:: to deal with floating point errors (just like archimedean)
// Copyright (C) 2011-2019 Zeno Rogue, see 'hyper.cpp' for details

namespace hr {
#if MAXMDIM >= 4
transmatrix cpush(int cid, ld alpha);
transmatrix cspin(int a, int b, ld alpha);
extern
vector<hpcshape> shWall3D, shMiniWall3D;
namespace binary {   
  void build_tmatrix(); 
  void virtualRebaseSimple(heptagon*& base, transmatrix& at);
  int celldistance3(heptagon *c1, heptagon *c2);
  hyperpoint deparabolic3(hyperpoint h);
  }

EX namespace reg3 {

  #if HDR
  inline short& altdist(heptagon *h) { return h->emeraldval; }
  #endif

  map<int, int> close_distances;

  EX int bucketer(ld x) {
    return int(x * 10 + 100000.5) - 100000;
    }
  
  int bucketer(hyperpoint h) {
    return bucketer(h[0]) + 1000 * bucketer(h[1]) + 1000000 * bucketer(h[2]);
    }
  
  int loop, face;

  EX vector<hyperpoint> cellshape;
  vector<hyperpoint> vertices_only;
  
  EX transmatrix spins[12], adjmoves[12];

  ld adjcheck, strafedist;
  EX bool dirs_adjacent[16][16];

  template<class T> ld binsearch(ld dmin, ld dmax, const T& f) {
    for(int i=0; i<200; i++) {
      ld d = (dmin + dmax) / 2;
      if(f(d)) dmax = d;
      else dmin = d;
      }
    return dmin;
    } 

  EX void generate() {
  
    if(S7 == 4) face = 3;
    if(S7 == 6) face = 4;
    if(S7 == 12) face = 5;
    if(S7 == 8) face = 3;
    /* icosahedron not implemented */
    loop = ginf[geometry].tiling_name[5] - '0';
    println(hlog, "face = ", face, " loop = ", loop, " S7 = ", S7);
    
    ld dual_angle = binsearch(0, M_PI, [&] (ld d) {
      hyperpoint h0 = cpush(0, 1) * C0;
      hyperpoint h1 = cspin(0, 1, d) * h0;
      hyperpoint h2 = cspin(1, 2, 2*M_PI/loop) * h1;
      return hdist(h0, h1) > hdist(h1, h2);
      });

    ld dodecahedron_angle = binsearch(0, M_PI, [&] (ld d) {
      hyperpoint h0 = cpush(0, 1) * C0;
      hyperpoint h1 = cspin(0, 1, d) * h0;
      hyperpoint h2 = cspin(1, 2, 2*M_PI/face) * h1;
      return hdist(h0, h1) > hdist(h1, h2);
      });
    
    if(S7 == 8) {
      /* 24-cell is a special case because it is the only one with '4' in the middle of the Schlaefli symbol. */
      /* The computations above assume 3 */
      hyperpoint h1 = hpxy3(.5,.5,.5);
      hyperpoint h2 = hpxy3(.5,.5,-.5);
      dual_angle = hdist(h1, h2);
      }
    
    println(hlog, "dodecahedron angle = ", dodecahedron_angle);
    println(hlog, "dual angle = ", dual_angle);
    
    ld inp_length = binsearch(0, 1.55, [&] (ld d) {
      hyperpoint h = xpush(-d) * spin(2*M_PI/face) * xpush0(d);
      ld alpha = M_PI - atan2(-h[1], h[0]);
      return (alpha < dual_angle / 2) ? hyperbolic : sphere;
      });
    
    println(hlog, "inp length = ", inp_length);
    
    ld edge_length = hdist(xpush0(inp_length), spin(2*M_PI/face) * xpush0(inp_length));
    if(S7 == 8) edge_length = hdist(normalize(hpxyz3(1,1,0,0)), normalize(hpxyz3(1,0,1,0))); 
    println(hlog, "edge length = ", edge_length);
    
    hyperpoint h0 = cpush(0, 1) * C0;
    hyperpoint h1 = cspin(0, 1, dodecahedron_angle) * h0;
    hyperpoint h2 = cspin(1, 2, 2*M_PI/face) * h1;
    hyperpoint h3 = cspin(1, 2, -2*M_PI/face) * h1;
    hyperpoint a2 = S7 == 8 ? normalize(h1 + h2) : normalize(h0 + h1 + h2);
    hyperpoint a3 = S7 == 8 ? normalize(h1 + h3) : normalize(h0 + h1 + h3);
    
    println(hlog, "S7 = ", S7);

    ld whereonline = binsearch(0, 5, [&] (ld d) {
      // sometimes breaks in elliptic
      dynamicval<eGeometry> g(geometry, elliptic ? gCell120 : geometry);
      hyperpoint z2 = a2 * d + C0 * (1-d);
      if(hyperbolic && intval(z2, Hypc) >= 0) return true;
      hyperpoint b2 = normalize(z2);
      hyperpoint z3 = a3 * d + C0 * (1-d);
      hyperpoint b3 = normalize(z3);
      return hdist(b2, b3) >= edge_length;
      });
    
    println(hlog, "whereonline = ", whereonline);
    a2 = normalize(a2 * whereonline + C0 * (1-whereonline));
    a3 = normalize(a3 * whereonline + C0 * (1-whereonline));

    hyperpoint mid = Hypc;
    for(int i=0; i<face; i++) mid += cspin(1, 2, 2*i*M_PI/face) * a2;
    mid = normalize(mid);
    ld between_centers = 2 * hdist0(mid);
    println(hlog, "between_centers = ", between_centers);

    if(S7 == 12 || S7 == 8) {
      spins[0] = Id;
      spins[1] = cspin(0, 1, dodecahedron_angle) * cspin(1, 2, M_PI);
      for(int a=2; a<face+1; a++) spins[a] = cspin(1, 2, 2*M_PI*(a-1)/face) * spins[1];      
      for(int a=S7/2; a<S7; a++) spins[a] = cspin(0, 1, M_PI) * spins[a-S7/2];
      }
    
    if(S7 == 6) {
      spins[0] = Id;
      spins[1] = cspin(0, 1, dodecahedron_angle) * cspin(1, 2, M_PI);
      spins[2] = cspin(1, 2, M_PI/2) * spins[1];
      for(int a=S7/2; a<S7; a++) spins[a] = spins[a-S7/2] * cspin(0, 1, M_PI);
      }
    
    if(S7 == 4) {
      spins[0] = Id;
      spins[1] = cspin(0, 1, dodecahedron_angle) * cspin(1, 2, M_PI);
      for(int a=2; a<face+1; a++) spins[a] = cspin(1, 2, 2*M_PI*(a-1)/face) * spins[1];      
      }
    
    cellshape.clear();
    for(int a=0; a<S7; a++)
    for(int b=0; b<face; b++)
      cellshape.push_back(spins[a] * cspin(1, 2, 2*M_PI*b/face) * a2);
    
    adjmoves[0] = cpush(0, between_centers) * cspin(0, 2, M_PI);
    for(int i=1; i<S7; i++) adjmoves[i] = spins[i] * adjmoves[0];

    for(int a=0; a<S7; a++)
      println(hlog, "center of ", a, " is ", tC0(adjmoves[a]));    
    
    println(hlog, "doublemove = ", tC0(adjmoves[0] * adjmoves[0]));

    adjcheck = hdist(tC0(adjmoves[0]), tC0(adjmoves[1])) * 1.0001;

    int numedges = 0;
    for(int a=0; a<S7; a++) for(int b=0; b<S7; b++) {
      dirs_adjacent[a][b] = a != b && hdist(tC0(adjmoves[a]), tC0(adjmoves[b])) < adjcheck;
      if(dirs_adjacent[a][b]) numedges++;
      }
    println(hlog, "numedges = ", numedges);
    
    if(loop == 4) strafedist = adjcheck;
    else strafedist = hdist(adjmoves[0] * C0, adjmoves[1] * C0);

    vertices_only.clear();
    for(hyperpoint h: cellshape) {
      bool found = false;
      for(hyperpoint h2: vertices_only) if(hdist(h, h2) < 1e-6) found = true;
      if(!found) vertices_only.push_back(h);
      }
    }

  void binary_rebase(heptagon *h, const transmatrix& V) {
    }
  
  void test();
    
  struct hrmap_field3 : hrmap {
    vector<heptagon*> allh;
    vector<cell*> acells;
    
    vector<vector<transmatrix>> tmatrices;
    
    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;
  
    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 = reg3::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);
  
      println(hlog, "full_P = ", full_P, " / ", R);
      println(hlog, "full_R = ", full_R, " / ", P);
      println(hlog, "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]) {
        println(hlog, i, ". ", fullmatrices[i]);
        }
      
      // find cav such that:
      // cav * Id          * C0 = corner0
      // cav * adjmoves[0] * C0 = corner1
      // cav * adjmoves[1] * C0 = corner3
      // cav * adjmoves[2] * C0 = cornerx
      
      hyperpoint corner0 = reg3::cellshape[0];
      hyperpoint corner1 = reg3::cellshape[1];
      hyperpoint corner3 = reg3::cellshape[3];
      hyperpoint cornerx;
  
      for(hyperpoint h: reg3::cellshape) println(hlog, "some corner ", h);
  
      for(hyperpoint h: reg3::cellshape)
        if(hdist(h, corner1) > .1 && hdist(h, corner3) > .1 && abs(hdist(h, corner0)-hdist(corner0, corner1)) < .1)
          cornerx = h;
      println(hlog, "corner0 = ", corner0);
      println(hlog, "corner1 = ", corner1);
      println(hlog, "corner3 = ", corner3);
      println(hlog, "cornerx = ", cornerx);
      
      transmatrix adj = Id, iadj = Id;
  
      geometry = g;
      reg3::generate();
      
      cyclers.clear();
      println(hlog, "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);
        println(hlog, "cav = ", cav);
        println(hlog, "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;
            println(hlog, i, " = ", make_tuple(a,P,b), " = ", code, " T = ", T);
            }
          }
        
      println(hlog, "field_adjmoves = ", field_adjmoves);
      
      println(hlog, "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);
        }
      
      println(hlog, "building allh...");
      int cells = N / perm_group;
      allh.resize(cells);
      for(int i=0; i<cells; i++) {
        allh[i] = tailored_alloc<heptagon> (S7);
        allh[i]->c7 = newCell(S7, allh[i]);
        allh[i]->fieldval = i;
        allh[i]->zebraval = 0;
        allh[i]->alt = NULL;
        acells.push_back(allh[i]->c7);
        }
      
      println(hlog, "finding tmatrices...");
      tmatrices.resize(cells);
      
      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) println(hlog, "bad found: ", i, "/", d, "/", found);
          // println(hlog, "tmatrix(",i,",",d,") = ", tmatrices[i][d]);
          }
        }
  
      println(hlog, "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);
          }
  
      create_patterns();
      }
    
    set<cellwalker> plane;

    void make_plane(cellwalker cw) {
      if(plane.count(cw)) return;
      plane.insert(cw);
      for(int i=0; i<S7; i++)
        if(reg3::dirs_adjacent[i][cw.spin])
          make_plane(reg3::strafe(cw, i));
      }
    
    
    void create_patterns() {
      // change the geometry to make sure that the correct celldistance is used
      dynamicval<eGeometry> g(geometry, S7 == 12 ? gField534 : gField435);
      // also, strafe needs currentmap
      dynamicval<hrmap*> c(currentmap, this);

      if(S7 == 12) {
        // Emerald in 534
        cell *a = gamestart();
        cell *b = a;
        for(cell *c: allcells())
          if(hr::celldistance(a, c) == 5) {
            b = c;
            break;
            }
        for(cell *c: allcells())
          if(hr::celldistance(a, c) > hr::celldistance(b, c))
            c->master->zebraval |= 1;
          
        // Vineyard in 534
        b = (cellwalker(a, 0) + wstep + rev + wstep).at;
        for(cell *c: allcells())
          if(hr::celldistance(a, c) == hr::celldistance(b, c))
            c->master->zebraval |= 2;
        }
      
      if(S7 == 6) {
        // Emerald in 534
        cell *a = gamestart();
        for(cell *c: allcells())
          if(hr::celldistance(a, c) > 3)
            c->master->zebraval |= 1;
        
        // Vineyard in 435
        make_plane(cellwalker(gamestart(), 0));
        println(hlog, "plane size = ", isize(plane));
        
        set<int> plane_indices;
        for(auto cw: plane) plane_indices.insert(cw.at->master->fieldval);

        set<int> nwi;
        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])];
            if(plane_indices.count(j)) ok = false;
            forCellEx(c1, allcells()[j]) if(plane_indices.count(c1->master->fieldval)) ok = false;
            }
          if(ok) nwi.insert(i);
          }

        int gpow = 0;
        
        for(int i: nwi) {
          int pw = 1;
          int at = i;
          while(true) {
            at = currfp_gmul(at, i);
            if(!nwi.count(at)) break;
            pw++;
            }
          if(pw == 4) gpow = i;
          }
        
        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])];
            allcells()[j]->master->zebraval |= 2;
            }
          u = currfp_gmul(u, gpow);
          }
        }
      }
    
    void draw() override {
      sphereflip = Id;
      
      // for(int i=0; i<S6; i++) queuepoly(ggmatrix(cwt.at), shWall3D[i], 0xFF0000FF);
      
      dq::visited_by_matrix.clear();
      dq::enqueue_by_matrix(viewctr.at, cview());
      
      while(!dq::drawqueue.empty()) {
        auto& p = dq::drawqueue.front();
        heptagon *h = get<0>(p);
        transmatrix V = get<1>(p);
        dynamicval<ld> b(band_shift, get<2>(p));
        bandfixer bf(V);
        dq::drawqueue.pop();
              
        cell *c = h->c7;
        if(!do_draw(c, V)) continue;
        drawcell(c, V, 0, false);
    
        for(int d=0; d<S7; d++)
          dq::enqueue_by_matrix(h->move(d), V * tmatrices[h->fieldval][d]);
        }
      }
  
    transmatrix relative_matrix(heptagon *h2, heptagon *h1) override {
      if(h1 == h2) return Id;
      int d = hr::celldistance(h2->c7, h1->c7);

      for(int a=0; a<S7; a++) if(hr::celldistance(h1->move(a)->c7, h2->c7) < d)
        return tmatrices[h1->fieldval][a] * relative_matrix(h2, h1->move(a));
      
      println(hlog, "error in hrmap_field3:::relative_matrix");
      return Id;
      }
  
    heptagon *getOrigin() override { return allh[0]; }
    
    vector<cell*>& allcells() override { return acells; }    

    vector<hyperpoint> get_vertices(cell* c) override {
      return vertices_only;
      }
    };

  struct hrmap_reg3 : hrmap {
  
    heptagon *origin;
    hrmap *binary_map;
    hrmap_field3 *quotient_map;
    
    unordered_map<heptagon*, pair<heptagon*, transmatrix>> reg_gmatrix;
    unordered_map<heptagon*, vector<pair<heptagon*, transmatrix> > > altmap;

    vector<cell*> spherecells;  

    vector<cell*>& allcells() override { 
      if(sphere) return spherecells;
      return hrmap::allcells();
      }

    hrmap_reg3() {
      generate();
      origin = tailored_alloc<heptagon> (S7);
      heptagon& h = *origin;
      h.s = hsOrigin;
      h.cdata = NULL;
      h.alt = NULL;
      h.distance = 0;
      h.fieldval = 0;
      h.c7 = newCell(S7, origin);
      if(sphere) spherecells.push_back(h.c7);
      worst_error1 = 0, worst_error2 = 0;
      
      dynamicval<hrmap*> cr(currentmap, this);
      
      heptagon *alt = NULL;
      transmatrix T = Id;
      
      if(hyperbolic) {
        #if CAP_FIELD
        quotient_map = new hrmap_field3;
        h.zebraval = quotient_map->allh[0]->zebraval;
        #endif

        dynamicval<eGeometry> g(geometry, gBinary3); 
        binary::build_tmatrix();
        alt = tailored_alloc<heptagon> (S7);
        alt->s = hsOrigin;
        alt->emeraldval = 0;
        alt->zebraval = 0;
        alt->distance = 0;
        alt->alt = alt;
        alt->cdata = NULL;
        alt->c7 = NULL;
        binary_map = binary::new_alt_map(alt);
        T = xpush(.01241) * spin(1.4117) * xpush(0.1241) * cspin(0, 2, 1.1249) * xpush(0.07) * Id;
        }
      else binary_map = NULL, quotient_map = NULL;
      
      reg_gmatrix[origin] = make_pair(alt, T);
      altmap[alt].emplace_back(origin, T);
      
      celllister cl(origin->c7, 4, 100000, NULL);
      for(cell *c: cl.lst) {
        hyperpoint h = tC0(relative_matrix(c->master, origin));
        close_distances[bucketer(h)] = cl.getdist(c);
        }
      }
    
    ld worst_error1, worst_error2;

    heptagon *getOrigin() override {
      return origin;
      }
    
    void fix_distances(heptagon *h, heptagon *h2) {
      vector<heptagon*> to_fix;

      auto fix_pair = [&] (heptagon *h, heptagon *h2) {
        if(!h2) return;
        if(h->distance > h2->distance+1) {
          h->distance = h2->distance + 1;
          to_fix.push_back(h);
          }
        else if(h2->distance > h->distance+1) {
          h2->distance = h->distance + 1;
          to_fix.push_back(h2);
          }
        if(h->alt && h->alt == h2->alt) {
          if(altdist(h) > altdist(h2) + 1) {
            altdist(h) = altdist(h2) + 1;
            to_fix.push_back(h);
            }
          else if (altdist(h2) > altdist(h) + 1) {
            altdist(h2) = altdist(h) + 1;
            to_fix.push_back(h2);
            }
          }
        };

      if(!h2) to_fix = {h};
      else fix_pair(h, h2);
      
      for(int i=0; i<isize(to_fix); i++) {
        h = to_fix[i];
        for(int j=0; j<S7; j++) fix_pair(h, h->move(j));
        }
      }
    
    #define DEB 0
    
    heptagon *counterpart(heptagon *h) {
      return quotient_map->allh[h->fieldval];
      }

    heptagon *create_step(heptagon *parent, int d) override {
      auto& p1 = reg_gmatrix[parent];
      if(DEB) println(hlog, "creating step ", parent, ":", d, ", at ", p1.first, tC0(p1.second));
      heptagon *alt = p1.first;
      #if CAP_FIELD
      transmatrix T = p1.second * (hyperbolic ? quotient_map->tmatrices[parent->fieldval][d] : adjmoves[d]);
      #else
      transmatrix T = p1.second * adjmoves[d];
      #endif
      transmatrix T1 = T;
      if(hyperbolic) {
        dynamicval<eGeometry> g(geometry, gBinary3);
        dynamicval<hrmap*> cm(currentmap, binary_map);
        binary::virtualRebaseSimple(alt, T);
        }
      
      fixmatrix(T);
      auto hT = tC0(T);
      
      if(DEB) println(hlog, "searching at ", alt, ":", hT);

      if(DEB) for(auto& p2: altmap[alt]) println(hlog, "for ", tC0(p2.second), " intval is ", intval(tC0(p2.second), hT));
      
      ld err;
      
      for(auto& p2: altmap[alt]) if((err = intval(tC0(p2.second), hT)) < 1e-3) {
        if(err > worst_error1) println(hlog, format("worst_error1 = %lg", double(worst_error1 = err)));
        // println(hlog, "YES found in ", isize(altmap[alt]));
        if(DEB) println(hlog, "-> found ", p2.first);
        int fb = 0;
        hyperpoint old = T * (inverse(T1) * tC0(p1.second));
        #if CAP_FIELD
        if(hyperbolic) {
          p2.first->c.connect(counterpart(parent)->c.spin(d), parent, d, false);
          fix_distances(p2.first, parent);
          return p2.first;
          }
        #endif
        for(int d2=0; d2<S7; d2++) {
          hyperpoint back = p2.second * tC0(adjmoves[d2]);
          if((err = intval(back, old)) < 1e-3) {
            if(err > worst_error2) println(hlog, format("worst_error2 = %lg", double(worst_error2 = err)));
            if(p2.first->move(d2)) println(hlog, "error: repeated edge");
            p2.first->c.connect(d2, parent, d, false);
            fix_distances(p2.first, parent);
            fb++;
            }
          }
        if(fb != 1) { 
          println(hlog, "found fb = ", fb); 
          println(hlog, old);
          for(int d2=0; d2<S7; d2++) {
            println(hlog, p2.second * tC0(adjmoves[d2]), " in distance ", intval(p2.second * tC0(adjmoves[d2]), old));
            }
          parent->c.connect(d, parent, d, false);
          return parent;
          }
        return p2.first;
        }
      
      if(DEB) println(hlog, "-> not found");
      int d2 = 0, fv = isize(reg_gmatrix);
      #if CAP_FIELD
      if(hyperbolic) {
        auto cp = counterpart(parent);
        d2 = cp->c.spin(d);
        fv = cp->c.move(d)->fieldval;
        }
      #endif
      heptagon *created = tailored_alloc<heptagon> (S7);
      created->c7 = newCell(S7, created);
      if(sphere) spherecells.push_back(created->c7);
      created->alt = NULL;
      created->cdata = NULL;
      #if CAP_FIELD
      if(hyperbolic) {
        created->zebraval = quotient_map->allh[fv]->zebraval;
        }
      else
      #endif
      created->zebraval = hrand(10);
      created->fieldval = fv;
      created->distance = parent->distance + 1;
      fixmatrix(T);
      reg_gmatrix[created] = make_pair(alt, T);
      altmap[alt].emplace_back(created, T);
      created->c.connect(d2, parent, d, false);
      return created;
      }

    ~hrmap_reg3() {
      if(binary_map) {        
        dynamicval<eGeometry> g(geometry, gBinary3);
        delete binary_map;
        }
      if(quotient_map) delete quotient_map;
      clearfrom(origin);
      }
    
    map<heptagon*, int> reducers;

    void link_alt(const cellwalker& hs) override {
      auto h = hs.at->master;
      altdist(h) = 0;
      if(h->alt->s != hsOrigin) reducers[h] = hs.spin;
      }
    
    void generateAlts(heptagon* h, int levs, bool link_cdata) override {
      if(reducers.count(h)) {
        heptspin hs(h, reducers[h]);
        reducers.erase(h);
        hs += wstep;
        hs += rev;
        altdist(hs.at) = altdist(h) - 1;
        hs.at->alt = h->alt;
        reducers[hs.at] = hs.spin;
        fix_distances(hs.at, NULL);
        }
      for(int i=0; i<S7; i++) {
        auto h2 = h->cmove(i);
        if(h2->alt == NULL) {
          h2->alt = h->alt;
          altdist(h2) = altdist(h) + 1;
          fix_distances(h2, NULL);
          }
        }
      }

    void draw() override {
      sphereflip = Id;
      
      // for(int i=0; i<S6; i++) queuepoly(ggmatrix(cwt.at), shWall3D[i], 0xFF0000FF);
      
      dq::visited.clear();
      dq::enqueue(viewctr.at, cview());
      
      while(!dq::drawqueue.empty()) {
        auto& p = dq::drawqueue.front();
        heptagon *h = get<0>(p);
        transmatrix V = get<1>(p);
        dynamicval<ld> b(band_shift, get<2>(p));
        bandfixer bf(V);
        dq::drawqueue.pop();
        
        
        cell *c = h->c7;
        if(!do_draw(c, V)) continue;
        drawcell(c, V, 0, false);
    
        for(int i=0; i<S7; i++) if(h->move(i)) {
          #if CAP_FIELD
          if(hyperbolic) dq::enqueue(h->move(i), V * quotient_map->tmatrices[h->fieldval][i]);
          else
          #endif
          dq::enqueue(h->move(i), V * relative_matrix(h->move(i), h));
          }
        }
      }
    
    transmatrix relative_matrix(heptagon *h2, heptagon *h1) override {
      auto p1 = reg_gmatrix[h1];
      auto p2 = reg_gmatrix[h2];
      transmatrix T = Id;
      if(hyperbolic) { 
        dynamicval<eGeometry> g(geometry, gBinary3);
        dynamicval<hrmap*> cm(currentmap, binary_map);
        T = binary_map->relative_matrix(p2.first, p1.first);
        }
      return inverse(p1.second) * T * p2.second;
      }
    
    vector<hyperpoint> get_vertices(cell* c) override {
      return vertices_only;
      }
    };
  
EX hrmap* new_map() {
  if(quotient) return new hrmap_field3;
  return new hrmap_reg3;
  }

hrmap_reg3* regmap() {
  return ((hrmap_reg3*) currentmap);
  }

EX int celldistance(cell *c1, cell *c2) {
  if(c1 == c2) return 0;
  if(c1 == currentmap->gamestart()) return c2->master->distance;
  if(c2 == currentmap->gamestart()) return c1->master->distance;

  auto r = regmap();

  hyperpoint h = tC0(r->relative_matrix(c1->master, c2->master));
  int b = bucketer(h);
  if(close_distances.count(b)) return close_distances[b];

  dynamicval<eGeometry> g(geometry, gBinary3);  
  return 20 + binary::celldistance3(r->reg_gmatrix[c1->master].first, r->reg_gmatrix[c2->master].first);
  }

EX bool pseudohept(cell *c) {
  auto m = regmap();
  if(sphere) {
    hyperpoint h = tC0(m->relative_matrix(c->master, regmap()->origin));
    if(S7 == 12) {
      hyperpoint h1 = cspin(0, 1, atan2(16, 69) + M_PI/4) * h;
      for(int i=0; i<4; i++) if(abs(abs(h1[i]) - .5) > .01) return false;
      return true;
      }
    if(S7 == 8)
      return h[3] >= .99 || h[3] <= -.99 || abs(h[3]) < .01; 
    if(loop == 3 && face == 3 && S7 == 4)
      return c == m->gamestart();
    if(loop == 4 && face == 3)
      return abs(h[3]) > .9;
    if(loop == 3 && face == 4)
      return abs(h[3]) > .9;
    if(loop == 5 && face == 3)
      return abs(h[3]) > .99 || abs(h[0]) > .99 || abs(h[1]) > .99 || abs(h[2]) > .99;
    }
  // chessboard pattern in 534
  if(geometry == gSpace534) 
    return c->master->distance & 1;
  if(geometry == gField534) 
    return hr::celldistance(c, currentmap->gamestart()) & 1;
  if(hyperbolic) {
    heptagon *h = m->reg_gmatrix[c->master].first;
    return (h->zebraval == 1) && (h->distance & 1);
    }    
  return false;
  }
#endif

#if 0
/* More precise, but very slow distance. Not used/optimized for now */

ld adistance(cell *c) {  
  hyperpoint h = tC0(regmap()->reg_gmatrix[c->master].second);
  h = binary::deparabolic3(h);
  return regmap()->reg_gmatrix[c->master].first->distance * log(2) - h[0];
  }

unordered_map<pair<cell*, cell*>, int> memo;

bool cdd;

int celldistance(cell *c1, cell *c2) {
  if(memo.count(make_pair(c1, c2))) return memo[make_pair(c1, c2)];
  if(c1 == c2) return 0;
  vector<cell*> v[2];
  v[0].push_back(c1);
  v[1].push_back(c2);
  
  int steps = 0;
  
  map<cell*, int> visited;
  visited[c1] = 1;
  visited[c2] = 2;
  
  while(true) {
    if(cdd) {
      println(hlog, "state ", steps, "/",isize(v[0]), "/", isize(v[1]));
      println(hlog, "  A: ", v[0]);
      println(hlog, "  B: ", v[1]); 
      }
    for(int i: {0,1}) {
      vector<cell*> new_v;
      for(cell *c: v[i]) forCellCM(cn, c) if(adistance(cn) < adistance(c)) {
        auto &vi = visited[cn];
        if((vi&3) == 0) {
          vi = 4 * (steps+1);
          vi |= (1<<i);
          new_v.push_back(cn);
          }
        else if((vi&3) == 2-i) {
          vector<pair<cell*, int>> ca1, ca2;
          int b1 = 4*steps-4;
          int b2 = ((vi>>2)<<2) - 4;
          for(auto p: visited) {
            if(cdd) println(hlog, p);
            int ps = p.second & 3;
            if(ps == 1+i && p.second >= b1)
              ca1.emplace_back(p.first, p.second/4);
            if(ps == 2-i && p.second >= b2 && p.second <= b2+8)
              ca2.emplace_back(p.first, p.second/4);
            }
          int bound = 1<<16;
          for(auto p1: ca1) for(auto p2: ca2) {
            hyperpoint h = tC0(relative_matrix(p1.first->master, p2.first->master));
            int b = bucketer(h);
            if(close_distances.count(b)) {
              int d = close_distances[b] + p1.second + p2.second;
              if(cdd) println(hlog, "candidate: close=", close_distances[b], p1, p2, "; h = ", h);
              if(d < bound) bound = d;
              }
            else if(cdd) println(hlog, "bucket missing");
            }
          return memo[make_pair(c1, c2)] = bound;
          return bound;
          }
        }
      v[i] = std::move(new_v);
      }
    steps++;
    }
  }

cellwalker target;
int tsteps;

int dist_alt(cell *c) {
  if(!target.at) {
    target = cellwalker(currentmap->gamestart(), 0);
    tsteps = 0;
    for(int i=0; i<30; i++) target += wstep, target += rev, tsteps++;
    }
  if(specialland == laCamelot) return reg3::celldistance(c, target.at);
  else {
    int d = reg3::celldistance(c, target.at) - tsteps;
    if(d < 10) target += wstep, target += rev, tsteps++;
    return d;
    }
  }
#endif

// Construct a cellwalker in direction j from cw.at, such that its direction is as close
// as possible to cw.spin. Assume that j and cw.spin are adjacent

#if MAXMDIM >= 4
EX cellwalker strafe(cellwalker cw, int j) {
  hyperpoint hfront = tC0(adjmoves[cw.spin]);
  transmatrix T = currentmap->relative_matrix(cw.at->cmove(j)->master, cw.at->master);
  for(int i=0; i<S7; i++) if(i != cw.at->c.spin(j))
    if(hdist(hfront, T * tC0(adjmoves[i])) < strafedist + .01)
      return cellwalker(cw.at->move(j), i);
  println(hlog, "incorrect strafe");
  exit(1);
  }
}
#endif
}