mesh: portable atan2 and cbrt for the quality measures (#115)

The platform maths libraries' atan2 and cbrt differ in the last bits
between macOS and glibc, so the dihedral angles and the shape measure a
mesh reports did too. The crate's own versions use + - * / and sqrt
only (correctly rounded everywhere) and agree with the platform's to a
few ulps.

Co-Authored-By: Claude Opus 5.5 <[email protected]>
Claude-Session: https://claude.ai/code/session_01UhhZMBTviDQxZvwnyP89fo

crates/mesh/src/quality/shape.rs

@@ -199,6 +199,6 @@ pub(crate) fn pyramid_scaled_jacobian(p: [P; 5]) -> f64 {
 pub(crate) fn jacobian_metrics(a: P, b: P, c: P) -> (f64, f64) {
     let d = det(a, b, c);
     let frob = dot(a, a) + dot(b, b) + dot(c, c);
-    let shape = if d > 0.0 { ratio(3.0 * d.cbrt().powi(2), frob, 0.0) } else { 0.0 };
+    let shape = if d > 0.0 { ratio(3.0 * crate::vec3::cbrt(d).powi(2), frob, 0.0) } else { 0.0 };
     (corner_det(a, b, c), shape)
 }

crates/mesh/src/tet/measure.rs

@@ -104,7 +104,7 @@ pub fn dihedral_range(x: &[P; 4]) -> (f64, f64) {
     let mut lo = f64::INFINITY;
     let mut hi = f64::NEG_INFINITY;
     for (sv, cv) in s.iter().zip(c) {
-        let a = sv.abs().atan2(cv).to_degrees();
+        let a = crate::vec3::atan2(sv.abs(), cv).to_degrees();
         lo = lo.min(a);
         hi = hi.max(a);
     }

crates/mesh/src/vec3.rs

@@ -21,7 +21,59 @@ pub(crate) fn det(a: P, b: P, c: P) -> f64 {
 
 /// The angle between two vectors in degrees (0 when either is zero).
 pub(crate) fn angle_deg(a: P, b: P) -> f64 {
-    norm(cross(a, b)).atan2(dot(a, b)).to_degrees()
+    atan2(norm(cross(a, b)), dot(a, b)).to_degrees()
+}
+
+/// `atan(t)` for `0 ≤ t ≤ 1` from `+ − × ÷ √` only: two half-angle reductions
+/// (`atan t = 2 atan(t / (1 + √(1 + t²)))`) bring `t` below tan(π/16) ≈ 0.199, where the Taylor
+/// series converges to double precision in 20 terms.
+fn atan_unit(t: f64) -> f64 {
+    let mut t = t;
+    for _ in 0..2 {
+        t /= 1.0 + (1.0 + t * t).sqrt();
+    }
+    let (t2, mut term, mut sum) = (t * t, t, 0.0);
+    for k in 0..20 {
+        let s = term / (2 * k + 1) as f64;
+        sum += if k % 2 == 0 { s } else { -s };
+        term *= t2;
+    }
+    4.0 * sum
+}
+
+/// The angle of `(x, y)` in radians, in `(−π, π]`, like [`f64::atan2`] but the same bit for bit
+/// on every platform: the platform maths libraries' `atan2` differ in the last bits (macOS and
+/// glibc), and a mesh or report must not (AGENTS.md, Determinism). 0 for `(0, 0)` and for
+/// non-finite input.
+pub(crate) fn atan2(y: f64, x: f64) -> f64 {
+    use std::f64::consts::{FRAC_PI_2, PI};
+    if !(x.is_finite() && y.is_finite()) || (x == 0.0 && y == 0.0) {
+        return 0.0;
+    }
+    let (ax, ay) = (x.abs(), y.abs());
+    let a = if ay <= ax { atan_unit(ay / ax) } else { FRAC_PI_2 - atan_unit(ax / ay) };
+    let a = if x < 0.0 { PI - a } else { a };
+    if y < 0.0 { -a } else { a }
+}
+
+/// The cube root, from `+ − × ÷` only (Newton's iteration from an exact power-of-two estimate),
+/// the same on every platform (see [`atan2`]).
+pub(crate) fn cbrt(x: f64) -> f64 {
+    if !x.is_finite() || x == 0.0 {
+        return x;
+    }
+    let a = x.abs();
+    // 2^(e/3) for a's binary exponent e: within a factor 2 of the root.
+    let e = ((a.to_bits() >> 52) & 0x7ff) as i64 - 1023;
+    let mut r = f64::from_bits((((e / 3) + 1023).clamp(1, 2046) as u64) << 52);
+    for _ in 0..60 {
+        let next = r - (r * r * r - a) / (3.0 * r * r);
+        if next == r || !next.is_finite() {
+            break;
+        }
+        r = next;
+    }
+    if x < 0.0 { -r } else { r }
 }
 pub(crate) fn add(a: P, b: P) -> P {
     [a[0] + b[0], a[1] + b[1], a[2] + b[2]]
@@ -38,3 +90,30 @@ pub(crate) fn lerp(a: P, b: P, t: f64) -> P {
 pub(crate) fn mid(a: P, b: P) -> P {
     [0.5 * (a[0] + b[0]), 0.5 * (a[1] + b[1]), 0.5 * (a[2] + b[2])]
 }
+
+#[cfg(test)]
+mod tests {
+    use super::*;
+
+    /// The portable atan2 and cbrt agree with the platform's to a few ulps.
+    #[test]
+    fn portable_functions_match_the_platform() {
+        let mut worst: f64 = 0.0;
+        for i in 0..2000 {
+            let a = -std::f64::consts::PI + i as f64 * (2.0 * std::f64::consts::PI / 1999.0);
+            for r in [1e-300, 1e-3, 1.0, 7.5, 1e200] {
+                let (y, x) = (r * a.sin(), r * a.cos());
+                worst = worst.max((atan2(y, x) - y.atan2(x)).abs());
+            }
+        }
+        assert!(worst < 4e-15, "{worst}");
+        assert_eq!(atan2(0.0, 0.0), 0.0);
+        assert_eq!(atan2(1.0, 0.0), std::f64::consts::FRAC_PI_2);
+        assert_eq!(atan2(0.0, -1.0), std::f64::consts::PI);
+        for x in [1e-300, 2.5e-7, 0.3, 1.0, 8.0, 27.0, 1234.5, 1e300, -64.0] {
+            let rel: f64 = (cbrt(x) - x.cbrt()).abs() / x.cbrt().abs();
+            assert!(rel < 4e-16, "{x}: {} vs {}", cbrt(x), x.cbrt());
+        }
+        assert_eq!(cbrt(0.0), 0.0);
+    }
+}