From a1a7f216e02d4a7033527c07d5974ab885a03a65 Mon Sep 17 00:00:00 2001 From: Tarik Moussa Date: Sun, 31 May 2026 16:07:35 +0200 Subject: [PATCH] fix(num): N5 stable triangle area (Kahan) in euclidean_cot_weights MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit The cotangent weights divided by 2·√(t12·t23·t31·l123). For a needle/cap (sliver) triangle one of the t-values is the difference of near-equal edge lengths → catastrophic cancellation, and the area under the sqrt loses precision, feeding large relative error into every cotangent weight, the Hessian, and the linear solve (numerical-stability audit N5). Replace the area computation with Kahan's stable side-length formula (sort a≥b≥c, evaluate ¼·√[(a+(b+c))(c−(a−b))(c+(a−b))(a+(b−c))]). The denominator is still exactly 8·Area for well-shaped triangles but accurate for slivers. The triangle-inequality guard and the cotangent numerators are unchanged. Test: CotWeights_SliverMatchesHighPrecisionReference cross-checks a thin triangle (apex ≈ 0.01 rad) against a long-double law-of-cosines reference. 292/292 CGAL tests pass. Co-Authored-By: Claude Opus 4.8 --- code/include/euclidean_hessian.hpp | 28 +++++++++++---- code/tests/cgal/test_euclidean_hessian.cpp | 41 ++++++++++++++++++++++ 2 files changed, 63 insertions(+), 6 deletions(-) diff --git a/code/include/euclidean_hessian.hpp b/code/include/euclidean_hessian.hpp index 83f473d..b4706ac 100644 --- a/code/include/euclidean_hessian.hpp +++ b/code/include/euclidean_hessian.hpp @@ -17,7 +17,8 @@ // │ side lengths lij = exp(Λ̃ij/2): │ // │ │ // │ t12 = −l12+l23+l31, t23 = l12−l23+l31, t31 = l12+l23−l31 │ -// │ l123 = l12+l23+l31, denom2 = 2·sqrt(t12·t23·t31·l123) = 8·Area │ +// │ l123 = l12+l23+l31, denom2 = 8·Area (Area via Kahan's stable │ +// │ side-length formula — see euclidean_cot_weights)│ // │ │ // │ Cotangent at vertex k (opposite t_opp, adjacent t_a and t_b): │ // │ cot_k = (t_opp · l123 − t_a · t_b) / denom2 │ @@ -48,6 +49,7 @@ #include #include #include +#include #include namespace conformallab { @@ -85,12 +87,26 @@ inline EuclCotWeights euclidean_cot_weights(double l12, double l23, double l31) if (t12 <= 0.0 || t23 <= 0.0 || t31 <= 0.0) return {0.0, 0.0, 0.0, false}; - const double l123 = l12 + l23 + l31; - const double denom2_sq = t12 * t23 * t31 * l123; - if (denom2_sq <= 0.0) return {0.0, 0.0, 0.0, false}; + const double l123 = l12 + l23 + l31; - // denom2 = 2·sqrt(t12·t23·t31·l123) = 8·Area - const double denom2 = 2.0 * std::sqrt(denom2_sq); + // Triangle area via Kahan's numerically stable formula. Sort the side + // lengths a ≥ b ≥ c and evaluate with the cancellation-avoiding grouping + // Area = ¼·√[ (a+(b+c))·(c−(a−b))·(c+(a−b))·(a+(b−c)) ]. + // This replaces the naive 2·√(t12·t23·t31·l123): that product of + // near-equal-length differences loses precision for needle/cap triangles, + // where it would feed large relative error straight into the cotangent + // weights and the linear solve. The result denom2 = 8·Area is identical + // up to rounding for well-shaped triangles, but accurate for slivers. + double a = l12, b = l23, c = l31; + if (a < b) std::swap(a, b); + if (a < c) std::swap(a, c); + if (b < c) std::swap(b, c); // now a ≥ b ≥ c + + const double kahan = (a + (b + c)) * (c - (a - b)) + * (c + (a - b)) * (a + (b - c)); + if (kahan <= 0.0) return {0.0, 0.0, 0.0, false}; + + const double denom2 = 8.0 * (0.25 * std::sqrt(kahan)); // = 8·Area // Formula: cot_k = (t_opp · l123 − t_a · t_b) / denom2 // cot1: t_opp=t23, t_a=t12, t_b=t31 (v1 opposite l23) diff --git a/code/tests/cgal/test_euclidean_hessian.cpp b/code/tests/cgal/test_euclidean_hessian.cpp index 80305c6..72e7ca6 100644 --- a/code/tests/cgal/test_euclidean_hessian.cpp +++ b/code/tests/cgal/test_euclidean_hessian.cpp @@ -61,6 +61,47 @@ TEST(EuclideanHessian, CotWeights_RightIsoscelesTriangle) EXPECT_NEAR(cw.cot3, 1.0, 1e-12); // 45° at v3 } +// ════════════════════════════════════════════════════════════════════════════ +// N5: cotangent weights on a SLIVER triangle match a high-precision reference. +// +// The area is now computed by Kahan's stable side-length formula instead of +// the naive 2·√(t12·t23·t31·l123). On a thin (sliver) triangle the naive +// product of near-equal-length differences loses precision, biasing every +// cotangent weight. Here we cross-check against an independent law-of-cosines +// reference in long double (80-bit on x86 CI — a genuine high-precision oracle; +// equal to double on ARM64, where it still serves as a regression cross-check). +// ════════════════════════════════════════════════════════════════════════════ + +TEST(EuclideanHessian, CotWeights_SliverMatchesHighPrecisionReference) +{ + // Thin isosceles: short base l12, two near-unit legs (apex angle ≈ 0.01 rad). + const double l12 = 0.01, l23 = 1.0, l31 = 1.0; + auto cw = euclidean_cot_weights(l12, l23, l31); + ASSERT_TRUE(cw.valid); + + // cot at the vertex opposite side `opp`, with adjacent sides s1, s2: + // cos = (s1² + s2² − opp²) / (2·s1·s2), sin = √((1−cos)(1+cos)). + auto cot_ref = [](long double opp, long double s1, long double s2) { + long double cosA = (s1 * s1 + s2 * s2 - opp * opp) / (2.0L * s1 * s2); + long double sinA = std::sqrt((1.0L - cosA) * (1.0L + cosA)); + return static_cast(cosA / sinA); + }; + const double c1 = cot_ref(l23, l12, l31); // v1 opposite l23 + const double c2 = cot_ref(l31, l12, l23); // v2 opposite l31 + const double c3 = cot_ref(l12, l23, l31); // v3 opposite l12 (needle angle) + + auto rel = [](double a, double b) { + return std::abs(a - b) / std::max(1.0, std::abs(b)); + }; + EXPECT_LT(rel(cw.cot1, c1), 1e-9); + EXPECT_LT(rel(cw.cot2, c2), 1e-9); + EXPECT_LT(rel(cw.cot3, c3), 1e-9); + + EXPECT_TRUE(std::isfinite(cw.cot1) && + std::isfinite(cw.cot2) && + std::isfinite(cw.cot3)); +} + // ════════════════════════════════════════════════════════════════════════════ // Hessian is symmetric: H[i,j] == H[j,i] // ════════════════════════════════════════════════════════════════════════════