From 6c9f7e71e1aba21ffbdd280615cf89ee7c946c3e Mon Sep 17 00:00:00 2001
From: Zeno Rogue <zeno@attnam.com>
Date: Sun, 22 Aug 2021 23:41:38 +0200
Subject: [PATCH] better handling of precision errors in arcm and arb

---
 arbitrile.cpp   | 10 ++--------
 archimedean.cpp |  4 ++--
 hyperpoint.cpp  | 16 ++++++++++++++++
 3 files changed, 20 insertions(+), 10 deletions(-)

diff --git a/arbitrile.cpp b/arbitrile.cpp
index 171c01b1..1022c043 100644
--- a/arbitrile.cpp
+++ b/arbitrile.cpp
@@ -930,15 +930,9 @@ struct hrmap_arbi : hrmap {
       alt = (heptagon*) s;
       }
 
-    for(auto& p2: altmap[alt]) if(id_of(p2.first) == co.sid && hdist(tC0(p2.second), tC0(T)) < 1e-2) {
+    for(auto& p2: altmap[alt]) if(id_of(p2.first) == co.sid && same_point_may_warn(tC0(p2.second), tC0(T))) {
       for(int oth=0; oth < p2.first->type; oth++) {
-        ld err = hdist(p2.second * xsh.vertices[oth], T * xsh.vertices[co.eid]);
-        if(err < 1e-2) {
-          static ld max_err = 0;
-          if(err > max_err) {
-            println(hlog, "err = ", err);
-            max_err = err;
-            }
+        if(same_point_may_warn(p2.second * xsh.vertices[oth], T * xsh.vertices[co.eid])) {
           h->c.connect(d, p2.first, oth%p2.first->type, co.mirror);
           return p2.first;
           }
diff --git a/archimedean.cpp b/archimedean.cpp
index e59eac96..f3580da6 100644
--- a/archimedean.cpp
+++ b/archimedean.cpp
@@ -649,14 +649,14 @@ struct hrmap_archimedean : hrmap {
       
     DEBB(DF_GEOM, ("look for: ", alt, " / ", T * C0));
   
-    for(auto& p2: altmap[alt]) if(intval(p2.second * C0, T * C0) < 1e-4) {
+    for(auto& p2: altmap[alt]) if(same_point_may_warn(p2.second * C0, T * C0)) {
       DEBB(DF_GEOM, ("cell found: ", 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);
         DEBB(DF_GEOM, ("compare: ", T1 * xpush0(1), ":: ", p2.second * xpush0(1)));
-        if(intval(T1 * xpush0(1), p2.second * xpush0(1)) < 1e-4) {
+        if(same_point_may_warn(T1 * xpush0(1), p2.second * xpush0(1))) {
         
           // T1 = p2.second
           // T * spin(pi+t2.first) == p2.second
diff --git a/hyperpoint.cpp b/hyperpoint.cpp
index 4f6ced20..63fea01d 100644
--- a/hyperpoint.cpp
+++ b/hyperpoint.cpp
@@ -1681,5 +1681,21 @@ EX bool clockwise(hyperpoint h1, hyperpoint h2) {
   return h1[0] * h2[1] > h1[1] * h2[0];
   }
 
+EX ld worst_precision_error;
+
+#if HDR
+struct hr_precision_error : hr_exception { hr_precision_error() : hr_exception("precision error") {} };
+#endif
+
+/** check if a and b are the same, testing for equality. Throw an exception or warning if not sure */
+EX bool same_point_may_warn(hyperpoint a, hyperpoint b) {
+  ld d = hdist(a, b);
+  if(d > 1e-2) return false;
+  if(d > 1e-3) throw hr_precision_error();
+  if(d > 1e-6 && worst_precision_error <= 1e-6)
+    addMessage("warning: precision errors are building up!");
+  if(d > worst_precision_error) worst_precision_error = d;
+  return true;
+  }
 
 }