# HG changeset patch
# User Tuomo Valkonen <tuomov@iki.fi>
# Date 1733509668 18000
# Node ID 56d609b70a3ba269fd1f27fed12a6e617a503d7b
# Parent  1d865db9d3e4fff983c3c5ac2a346a5450fcb8ab
Cylinder numerical stability hacks

diff -r 1d865db9d3e4 -r 56d609b70a3b src/cylinder.rs
--- a/src/cylinder.rs	Fri Dec 06 13:10:57 2024 -0500
+++ b/src/cylinder.rs	Fri Dec 06 13:27:48 2024 -0500
@@ -501,10 +501,16 @@
             let sideedge = SidePoint{ angle : φ, z };
             let t1 = cap.log(&capedge);
             let t2 = sideedge.log(r, side);
-            // Either option should give the same result, but the first one avoids division.
-            t1 + capedge.tangent_from_side(top, t2)
-            // let n = t1.norm(L2);
-            // (t1/n)*(n + t2.norm(L2))
+            // Either option should give the same result.
+            // The first one avoids division, but seems numerically more unstable due to
+            // sines and cosines.
+            //t1 + capedge.tangent_from_side(top, t2)
+            let n = t1.norm(L2);
+            if n < f64::EPSILON {
+                capedge.tangent_from_side(top, t2)
+            } else {
+                (t1/n)*(n + t2.norm(L2))
+            }
         }
     }
 
@@ -528,10 +534,16 @@
             let sideedge = SidePoint{ angle : φ, z };
             let t1 = side.log(r, &sideedge);
             let t2 = capedge.log(cap);
-            // Either option should give the same result, but the first one avoids division.
-            t1 + capedge.tangent_to_side(top, t2)
-            // let n = t1.norm(L2);
-            // (t1/n)*(n + t2.norm(L2))
+            // Either option should give the same result.
+            // The first one avoids division, but seems numerically more unstable due to
+            // sines and cosines.
+            // t1 + capedge.tangent_to_side(top, t2)
+            let n = t1.norm(L2);
+            if n < f64::EPSILON {
+                capedge.tangent_to_side(top, t2)
+            } else {
+                (t1/n)*(n + t2.norm(L2))
+            }
         }
     }
 
@@ -563,18 +575,24 @@
             let t1 = side1.log(r, &sideedge1);
             let t2 = capedge1.log(&capedge2);
             let t3 = sideedge2.log(r, &side2);
-            // Any option should give the same result, but the first one avoids division.
-            // t1 + capedge1.tangent_to_side(top, t2 + capedge2.tangent_from_side(top, t3))
-            //
-            // let n = t2.norm(L2);
-            // t1 + capedge1.tangent_to_side(top, (t2/n)*(n + t3.norm(L2)))
+            // Either option should give the same result.
+            // The first one avoids division, but seems numerically more unstable due to
+            // sines and cosines.
             //
+            // t1 + capedge1.tangent_to_side(top, t2 + capedge2.tangent_from_side(top, t3))
             let n = t1.norm(L2);
-            (t1/n)*(n + t2.norm(L2) + t3.norm(L2))
-            //
-            // let n = t1.norm(L2);
-            // let t23 = t2 + capedge2.tangent_from_side(top, t3);
-            // (t1/n)*(n + t23.norm(L2))
+            if n < f64::EPSILON {
+                let n2 = t2.norm(L2);
+                capedge1.tangent_to_side(top,
+                    if n2 < f64::EPSILON {
+                        capedge2.tangent_from_side(top, t3)
+                    } else {
+                        (t2/n2)*(n2 + t3.norm(L2))
+                    }
+                )
+            } else {
+                (t1/n)*(n + t2.norm(L2) + t3.norm(L2))
+            }
         }
     }
 
@@ -606,10 +624,23 @@
             let t1 = cap1.log(&capedge1);
             let t2 = sideedge1.log(r, &sideedge2);
             let t3 = capedge2.log(cap2);
-            // Either option should give the same result, but the first one avoids division.
-            t1 + capedge1.tangent_from_side(top1, t2 + capedge2.tangent_to_side(top2, t3))
-            //let n = t1.norm(L2);
-            //(t1/n)*(n + t2.norm(L2) + t3.norm(L2))
+            // Either option should give the same result.
+            // The first one avoids division, but seems numerically more unstable due to
+            // sines and cosines.
+            // t1 + capedge1.tangent_from_side(top1, t2 + capedge2.tangent_to_side(top2, t3))
+            let n = t1.norm(L2);
+            if n < f64::EPSILON {
+                let n2 = t2.norm(L2);
+                capedge1.tangent_from_side(top1,
+                    if n2 < f64::EPSILON {
+                        capedge2.tangent_to_side(top2, t3)
+                    } else {
+                        (t2/n2)*(n2 + t3.norm(L2))
+                    }
+                )
+            } else {
+                (t1/n)*(n + t2.norm(L2) + t3.norm(L2))
+            }
         }
     }
 
@@ -851,7 +882,7 @@
 mod tests {
     use super::*;
 
-    static TOL : f64 = 1e-8;
+    static TOL : f64 = 1e-7;
 
     macro_rules! check_distance {
         ($distance : expr, $expected : expr) => {