From 3f8194746a9c58ab7ab2dd0b96c84ed1dd94b0ee Mon Sep 17 00:00:00 2001
From: Zeno Rogue <zeno@attnam.com>
Date: Sat, 10 Nov 2018 11:51:47 +0100
Subject: [PATCH] made Archimedean more stable

---
 archimedean.cpp | 50 ++++++++++++++++++++++++++++++++++++++-----------
 1 file changed, 39 insertions(+), 11 deletions(-)

diff --git a/archimedean.cpp b/archimedean.cpp
index 1a9ad916..292ce474 100644
--- a/archimedean.cpp
+++ b/archimedean.cpp
@@ -520,6 +520,20 @@ void connectHeptagons(heptspin hi, heptspin hs) {
   // heptagon *hnew = build_child(h, d, get_adj(h, d).first, get_adj(h, d).second);
   }
 
+// T and X are supposed to be equal -- move T so that it is closer to X
+void fixup_matrix(transmatrix& T, const transmatrix& X, ld step) {
+  for(int i=0; i<3; i++)
+  for(int j=0; j<3; j++)
+    T[i][j] = (T[i][j] * (1-step) + X[i][j] * step);
+
+  /*
+  for(int i=0; i<3; i++)
+  for(int j=0; j<3; j++)
+    if(T[i][j] - X[i][j] > 1e-3) exit(1);
+  */
+  fixmatrix(T);
+  }
+
 void create_adjacent(heptagon *h, int d) {
 
   SDEBUG( printf("%p.%d ~ ?\n", h, d); )
@@ -532,31 +546,45 @@ void create_adjacent(heptagon *h, int d) {
 
   // * spin(-tri[id][pi+i].first) * xpush(t.second) * pispin * spin(tri[id'][p'+d'].first)
   
-  auto& p = archimedean_gmatrix[h];
+  auto& p1 = archimedean_gmatrix[h];
   
-  heptagon *alt = p.first;
+  heptagon *alt = p1.first;
 
-  transmatrix T = p.second * spin(-t1.first) * xpush(t1.second);
+  transmatrix T = p1.second * spin(-t1.first) * xpush(t1.second);
+  transmatrix U = Id;
   
   if(hyperbolic) {
     dynamicval<eGeometry> g(geometry, gNormal); 
+    U = T;
     virtualRebaseSimple(alt, T);
+    U = U * inverse(T);
     }
   
   if(euclid) 
     alt = encodeId(pair_to_vec(int(T[0][2]), int(T[1][2])));
     
   SDEBUG( printf("look for: %p / %s\n", alt, display(T * C0)); )
-  
-  for(auto& p: altmap[alt]) if(intval(p.second * C0, T * C0) < 1e-4) {
-    SDEBUG( printf("cell found: %p\n", p.first); )
-    for(int d2=0; d2<p.first->degree(); d2++) {
-      heptspin hs(p.first, d2);
-      auto& t2 = current.get_triangle(p.first, d2);
+
+  for(auto& p2: altmap[alt]) if(intval(p2.second * C0, T * C0) < 1e-4) {
+    SDEBUG( printf("cell found: %p\n", p2.first); )
+    for(int d2=0; d2<p2.first->degree(); d2++) {
+      heptspin hs(p2.first, d2);
+      auto& t2 = current.get_triangle(p2.first, d2);
       transmatrix T1 = T * spin(M_PI + t2.first);
       SDEBUG( printf("compare: %s", display(T1 * xpush0(1)));  )
-      SDEBUG( printf(":: %s\n", display(p.second * xpush0(1)));  )
-      if(intval(T1 * xpush0(1), p.second * xpush0(1)) < 1e-4) {
+      SDEBUG( printf(":: %s\n", display(p2.second * xpush0(1)));  )
+      if(intval(T1 * xpush0(1), p2.second * xpush0(1)) < 1e-4) {
+      
+        // T1 = p2.second
+        // T * spin(pi+t2.first) == p2.second
+        // p1.second * spinm(-t1.first) * xpush(t1.second) * spin(pi+t2.first) == p2.second
+        
+        // bring p1 and p2 closer, to prevent floating point errors
+        if(hyperbolic) {
+          fixup_matrix(p1.second, U * p2.second * spin(-M_PI - t2.first) * xpush(-t1.second) * spin(t1.first), 0.25);
+          fixup_matrix(p2.second, T1, 0.25);
+          }
+
         while(skip_digons(hs, -1)) hs--;
         connectHeptagons(hi, hs);
         connect_digons_too(hi, hs);