conformallab/ConformalLabpp
3
0
Fork 0
Code Issues Pull Requests Actions Packages Projects Releases 1 Wiki Activity

fix(num): N5 stable triangle area (Kahan) in euclidean_cot_weights

Browse Source 
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 <noreply@anthropic.com>
This commit is contained in:
Tarik Moussa
2026-05-31 16:07:35 +02:00
parent b29eb4e0c4
commit b7c7d108b5
2 changed files with 63 additions and 6 deletions
Download Patch File Download Diff File Expand all files Collapse all files

28
code/include/euclidean_hessian.hpp
View File

@@ -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 <vector>
#include <array>
#include <cmath>
#include <algorithm>
#include <stdexcept>
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)

41
code/tests/cgal/test_euclidean_hessian.cpp
View File

@@ -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<double>(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]
// ════════════════════════════════════════════════════════════════════════════