mesh: platform-independent maths (libm) so meshes are bit-identical on macOS and Linux

The system maths library's sin, cos, atan2, log … are not correctly rounded and differ
between platforms: mesh.hex.generate gave the plate of the strong-relax test a different
mesh on Linux than on macOS (nodes from 2809 on), so the test passed on a Mac and failed
in CI. Every transcendental call in fieldcraft-mesh now goes through the libm crate (a
pure-Rust port of musl's), which gives the same bits everywhere.

Part of #438.

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

Cargo.lock

@@ -1418,6 +1418,7 @@ name = "fieldcraft-mesh"
 version = "0.6.0"
 dependencies = [
  "fieldcraft-model",
+ "libm",
  "num-bigint",
  "num-rational",
  "num-traits",

crates/mesh/Cargo.toml

@@ -10,6 +10,10 @@ workspace = true
 
 [dependencies]
 fieldcraft-model = { workspace = true }
+# Pure-Rust maths (a port of musl's libm) for sin, cos, log, exp, pow, hypot, acos: the same
+# bits on every platform, where the system library's are not correctly rounded and differ (the
+# plate's hex mesh differed between macOS and Linux, #438). atan2 and cbrt are crate::vec3's.
+libm = "0.2"
 thiserror = { workspace = true }
 
 [dev-dependencies]

crates/mesh/src/block.rs

@@ -146,7 +146,7 @@ impl Path {
             if norm(sub(back, p)) > 1e-9 * la.max(lb) || ((c * c + s * s).sqrt() - 1.0).abs() > 1e-9 {
                 return Err(Error::Invalid(format!("corner {p:?} is not on the edge's ellipse or arc about {center:?}")));
             }
-            Ok(s.atan2(c))
+            Ok(crate::vec3::atan2(s, c))
         };
         let start = angle(from)?;
         let mut sweep = angle(to)? - start;
@@ -170,7 +170,7 @@ impl Path {
             Path::Line { from, to } => add_scaled(from, s, sub(to, from)),
             Path::Ellipse { center, a, b, start, sweep, .. } => {
                 let t = start + s * sweep;
-                add_scaled(add_scaled(center, t.cos(), a), t.sin(), b)
+                add_scaled(add_scaled(center, libm::cos(t), a), libm::sin(t), b)
             }
         }
     }
@@ -179,7 +179,7 @@ impl Path {
 /// The parameters of `n` elements along one direction with last/first size ratio `bias`, at
 /// `order` points per element (`order · n + 1` values from 0 to 1; mid-side values halfway).
 fn parameters(n: usize, bias: f64, order: usize) -> Vec<f64> {
-    let q = if n > 1 { bias.powf(1.0 / (n - 1) as f64) } else { 1.0 };
+    let q = if n > 1 { libm::pow(bias, 1.0 / (n - 1) as f64) } else { 1.0 };
     let mut ends = Vec::with_capacity(n + 1);
     let (mut sum, mut size) = (0.0, 1.0);
     ends.push(0.0);
@@ -645,7 +645,7 @@ mod tests {
         ];
         let m = block(&Block { kind: ElementKind::Hex20, corners, divisions: [3, 6, 2], bias: [1.0; 3], curves }).unwrap();
         for p in &m.points {
-            let r = p[0].hypot(p[1]);
+            let r = libm::hypot(p[0], p[1]);
             assert!(r > a - 1e-12 && r < b + 1e-12, "{p:?}");
         }
         let exact = std::f64::consts::PI / 4.0 * (b * b - a * a) * 0.5;

crates/mesh/src/hex/brep.rs

@@ -372,7 +372,7 @@ fn composite_tags(surface: &BoundarySurface, surfaces: &[Option<Surface>]) -> Ve
             unit(cross(sub(b, a), sub(c, a)))
         })
     };
-    let cos_min = SMOOTH_ANGLE.cos();
+    let cos_min = libm::cos(SMOOTH_ANGLE);
     // Per pair of faces (low tag, high tag): whether every segment between them is smooth.
     let mut pairs: BTreeMap<(u32, u32), bool> = BTreeMap::new();
     for fe in surface.feature_edges().iter().filter(|e| e.kind == FeatureKind::FaceBoundary) {

crates/mesh/src/hex/decompose.rs

@@ -166,7 +166,7 @@ fn sub(a: Q, b: Q) -> Q {
 }
 
 fn len(a: Q) -> f64 {
-    a[0].hypot(a[1])
+    libm::hypot(a[0], a[1])
 }
 
 fn cross(a: Q, b: Q) -> f64 {
@@ -186,14 +186,17 @@ fn unit(a: Q) -> Option<Q> {
 /// region on the left: π on a straight line, less at a convex corner, more at a concave one.
 fn interior(a: Q, o: Q, b: Q) -> f64 {
     let (u, v) = (sub(o, a), sub(b, o));
-    std::f64::consts::PI - cross(u, v).atan2(dot(u, v))
+    std::f64::consts::PI - crate::vec3::atan2(cross(u, v), dot(u, v))
 }
 
 /// The inward bisector at a corner of interior angle `theta` whose outgoing edge runs along unit
 /// `w`: `w` turned left by half the angle.
 fn bisector(w: Q, theta: f64) -> Q {
     let half = 0.5 * theta;
-    [w[0] * half.cos() - w[1] * half.sin(), w[0] * half.sin() + w[1] * half.cos()]
+    {
+        let (s, c) = (libm::sin(half), libm::cos(half));
+        [w[0] * c - w[1] * s, w[0] * s + w[1] * c]
+    }
 }
 
 fn centroid(pts: &[Q]) -> Q {
@@ -1320,7 +1323,7 @@ mod tests {
     /// The half crank's cross-section (y ≤ 0) as the extrusion projects it: the outline, the
     /// bosses' circles where they rise above the arm, and the pocket's outline.
     fn crank_half() -> (Profile, Vec<Vec<Q>>) {
-        let deg = |p: Q, c: Q| (p[1] - c[1]).atan2(p[0] - c[0]).to_degrees();
+        let deg = |p: Q, c: Q| crate::vec3::atan2(p[1] - c[1], p[0] - c[0]).to_degrees();
         let spindle = [0.0, 0.0];
         let pedal = [170.0, 0.0];
         let j1 = [13.455063072993756, -18.65372021077189];
@@ -1344,7 +1347,7 @@ mod tests {
         };
         let outer = vec![
             a(spindle, 23.0, [-23.0, 0.0], j1, true),
-            a(c1, (j1[0] - c1[0]).hypot(j1[1] - c1[1]), j1, r1, false),
+            a(c1, libm::hypot(j1[0] - c1[0], j1[1] - c1[1]), j1, r1, false),
             line(r1, r2),
             a(c2, 6.0, r2, j2, false),
             a(pedal, 15.0, j2, [185.0, 0.0], true),

crates/mesh/src/hex/extrude_to.rs

@@ -410,7 +410,7 @@ fn section_boundary(quads: &[[u32; 4]]) -> Result<SectionBoundary> {
 /// The layer parameters `0 = t_0 < … < t_n = 1`, thicknesses in geometric progression with the
 /// last over the first equal to `bias`.
 pub(crate) fn parameters(n: usize, bias: f64) -> Vec<f64> {
-    let r = if n > 1 { bias.powf(1.0 / (n - 1) as f64) } else { 1.0 };
+    let r = if n > 1 { libm::pow(bias, 1.0 / (n - 1) as f64) } else { 1.0 };
     let mut out = Vec::with_capacity(n + 1);
     let (mut sum, mut h) = (0.0, 1.0);
     out.push(0.0);

crates/mesh/src/hex/geom.rs

@@ -101,7 +101,7 @@ impl Surface {
                 let h = dot(w, axis);
                 let radial = unit(sub(w, scale(axis, h)))?;
                 let rho = norm(sub(w, scale(axis, h)));
-                let (c, s) = (half_angle.cos(), half_angle.sin());
+                let (c, s) = (libm::cos(half_angle), libm::sin(half_angle));
                 let t = (h * c + rho * s).max(0.0);
                 add(apex, add(scale(axis, t * c), scale(radial, t * s)))
             }
@@ -128,7 +128,7 @@ impl Surface {
             Surface::Cone { apex, axis, half_angle } => {
                 let w = sub(p, apex);
                 let radial = unit(sub(w, scale(axis, dot(w, axis))))?;
-                unit(sub(scale(radial, half_angle.cos()), scale(axis, half_angle.sin())))
+                unit(sub(scale(radial, libm::cos(half_angle)), scale(axis, libm::sin(half_angle))))
             }
             Surface::Sphere { center, .. } => unit(sub(p, center)),
             Surface::Torus { center, axis, major, .. } => {
@@ -383,7 +383,7 @@ impl<'a> Inside<'a> {
             let (a, b, c) = (sub(*a, x), sub(*b, x), sub(*c, x));
             let (la, lb, lc) = (norm(a), norm(b), norm(c));
             let den = la * lb * lc + dot(a, b) * lc + dot(a, c) * lb + dot(b, c) * la;
-            sum += 2.0 * dot(a, cross(b, c)).atan2(den);
+            sum += 2.0 * crate::vec3::atan2(dot(a, cross(b, c)), den);
         }
         Ok(sum / (4.0 * std::f64::consts::PI))
     }

crates/mesh/src/hex/layered.rs

@@ -98,8 +98,8 @@ impl Frame {
                 let w = sub(p, o);
                 let z = dot(w, a);
                 let (x, y) = (dot(w, e1), dot(w, e2));
-                let t = y.atan2(x);
-                ([z, x.hypot(y)], if t < 0.0 { t + TAU } else { t })
+                let t = crate::vec3::atan2(y, x);
+                ([z, libm::hypot(x, y)], if t < 0.0 { t + TAU } else { t })
             }
         }
     }
@@ -107,7 +107,7 @@ impl Frame {
     pub(crate) fn point(&self, q: Q, h: f64) -> P {
         match *self {
             Frame::Translate { o, u, v, d } => add(o, add(add(scale(u, q[0]), scale(v, q[1])), scale(d, h))),
-            Frame::Rotate { o, a, e1, e2 } => add(o, add(scale(a, q[0]), scale(add(scale(e1, h.cos()), scale(e2, h.sin())), q[1]))),
+            Frame::Rotate { o, a, e1, e2 } => add(o, add(scale(a, q[0]), scale(add(scale(e1, libm::cos(h)), scale(e2, libm::sin(h))), q[1]))),
         }
     }
 

crates/mesh/src/hex/mod.rs

@@ -479,7 +479,7 @@ impl CellGrid {
             }
             (lo, hi)
         });
-        let per_axis = ((ctxs.len() as f64).cbrt().ceil() as usize * 2).clamp(1, MAX_GRID_DIM);
+        let per_axis = (crate::vec3::cbrt(ctxs.len() as f64).ceil() as usize * 2).clamp(1, MAX_GRID_DIM);
         let extent = (0..3).map(|k| hi[k] - lo[k]).fold(0.0, f64::max);
         let step = if extent > 0.0 && extent.is_finite() { extent / per_axis as f64 } else { 1.0 };
         let dim = |k: usize| (((hi[k] - lo[k]) / step).ceil() as usize).clamp(1, MAX_GRID_DIM);

crates/mesh/src/hex/region.rs

@@ -68,7 +68,7 @@ fn angles(lp: &[Vec<Q>], ccw: bool) -> Option<Vec<f64>> {
         let next = lp.get(k)?;
         let (a, o, b) = (*prev.get(prev.len().checked_sub(2)?)?, *next.first()?, *next.get(1)?);
         let (u, v) = ([o[0] - a[0], o[1] - a[1]], [b[0] - o[0], b[1] - o[1]]);
-        let turn = (u[0] * v[1] - u[1] * v[0]).atan2(u[0] * v[0] + u[1] * v[1]);
+        let turn = crate::vec3::atan2(u[0] * v[1] - u[1] * v[0], u[0] * v[0] + u[1] * v[1]);
         out.push(PI - if ccw { turn } else { -turn });
     }
     Some(out)
@@ -443,8 +443,11 @@ fn disk(loops: &[Vec<Q>], size: &dyn SizeField, opts: &QuadOptions) -> Result<Op
     // ∫y², ∫x² (both × 12) and ∫xy (× 24), signed as the loop: the region's long direction is
     // the principal axis of the largest ∫ along², the angle of 0.5 · atan2(2 ∫xy, ∫x² − ∫y²).
     let sign = a2.signum();
-    let theta = 0.5 * (sign * ixy).atan2(sign * (iyy - ixx));
-    let (u, v) = ([theta.cos(), theta.sin()], [-theta.sin(), theta.cos()]);
+    let theta = 0.5 * crate::vec3::atan2(sign * ixy, sign * (iyy - ixx));
+    let (u, v) = {
+        let (s, c) = (libm::sin(theta), libm::cos(theta));
+        ([c, s], [-s, c])
+    };
     let local = |p: Q| {
         let d = [p[0] - centre[0], p[1] - centre[1]];
         [d[0] * u[0] + d[1] * u[1], d[0] * v[0] + d[1] * v[1]]

crates/mesh/src/octree/mod.rs

@@ -462,7 +462,7 @@ fn to_frame(f: &[[f64; 3]; 3], p: P) -> P {
 /// The grid axes (rows): the coordinate axes, or the normals of the two largest orthogonal planar
 /// patches when more boundary area is aligned with them.
 fn choose_frame(g: &Geometry) -> [[f64; 3]; 3] {
-    let tol = 1f64.to_radians().sin();
+    let tol = libm::sin(1f64.to_radians());
     let normal = |t: u32| g.normals.get(t as usize).copied().unwrap_or([0.0; 3]);
     let mut planar: Vec<(f64, usize, P)> = Vec::new();
     for (i, p) in g.patches.iter().enumerate() {
@@ -501,7 +501,7 @@ fn choose_frame(g: &Geometry) -> [[f64; 3]; 3] {
 /// that only measures how near the edge it started; the faces of a wall face each other.
 fn triangle_thickness(g: &Geometry) -> Vec<f64> {
     let far = 2.0 * norm(sub(g.bounds.hi, g.bounds.lo)) + 1.0;
-    let opposed = -WALL_OPENING.to_radians().cos();
+    let opposed = -libm::cos(WALL_OPENING.to_radians());
     (0..g.surface.triangles().len() as u32)
         .map(|t| {
             let n = g.normals.get(t as usize).copied().unwrap_or([0.0; 3]);
@@ -625,7 +625,7 @@ impl Context<'_> {
             return Ok(true);
         }
         let reach = 0.5 * 3f64.sqrt() * e * 1.05;
-        let cos = self.options.curvature_angle.to_radians().cos();
+        let cos = libm::cos(self.options.curvature_angle.to_radians());
         let mut first: BTreeMap<u32, P> = BTreeMap::new();
         for t in self.g.overlapping(&bx) {
             let [a, b, c] = self.g.triangle(t);

crates/mesh/src/planar.rs

@@ -266,7 +266,7 @@ pub fn face_boundary(positions: &[[f64; 3]], triangles: &[[u32; 3]], tri_face: &
             .map(|i| {
                 let (p, q, r) = (at(i + n - 1), at(i), at(i + 1));
                 let (u, w) = ([q[0] - p[0], q[1] - p[1]], [r[0] - q[0], r[1] - q[1]]);
-                let turn = (u[0] * w[1] - u[1] * w[0]).atan2(u[0] * w[0] + u[1] * w[1]).abs();
+                let turn = crate::vec3::atan2(u[0] * w[1] - u[1] * w[0], u[0] * w[0] + u[1] * w[1]).abs();
                 ids[i] != ids[(i + n - 1) % n] || turn > CORNER_ANGLE
             })
             .collect();

crates/mesh/src/project/mod.rs

@@ -551,7 +551,7 @@ pub fn project_mesh(t: &Topology, surface: &BoundarySurface, g: &dyn Geometry, o
                     *m = m.max(measured);
                 }
                 if let Some(l) = class_log.get_mut(cls[d]).filter(|_| len > tol) {
-                    l.0 += measured.ln();
+                    l.0 += libm::log(measured);
                     l.1 += 1;
                 }
             }
@@ -598,7 +598,7 @@ pub fn project_mesh(t: &Topology, surface: &BoundarySurface, g: &dyn Geometry, o
                             let Some(nx) = g.normal(fa, qx) else { continue };
                             let h = dist(pa, qx);
                             if h > 0.5 * probe {
-                                kappa = kappa.max(crate::vec3::dot(na, nx).clamp(-1.0, 1.0).acos() / h);
+                                kappa = kappa.max(libm::acos(crate::vec3::dot(na, nx).clamp(-1.0, 1.0)) / h);
                             }
                         }
                     }
@@ -616,7 +616,7 @@ pub fn project_mesh(t: &Topology, surface: &BoundarySurface, g: &dyn Geometry, o
         .zip(&class_max)
         .zip(&class_arc)
         .map(|((&(sum, k), &longest), &arc)| {
-            let mean = if k == 0 { 1.0 } else { (sum / k as f64).exp().round() };
+            let mean = if k == 0 { 1.0 } else { libm::exp(sum / k as f64).round() };
             (mean.max((longest / MAX_STRETCH).ceil()) as usize).max(arc).clamp(1, crate::block::MAX_DIVISIONS)
         })
         .collect();

crates/mesh/src/relax/strong.rs

@@ -65,7 +65,7 @@ enum Phase {
 /// The objective's term for one element of quality `q` whose size went from `v0` to `v`.
 fn phi(phase: Phase, q: f64, valid: bool, v: f64, v0: f64) -> f64 {
     if valid && q > 0.0 {
-        let g = if v > 0.0 && v0 > 0.0 { (v / v0).ln() } else { 0.0 };
+        let g = if v > 0.0 && v0 > 0.0 { libm::log(v / v0) } else { 0.0 };
         let shape = match phase {
             Phase::Smooth => (1.0 / q - 1.0).powi(2),
             Phase::MinMax => q.powi(-MIN_MAX_POWER) - 1.0,

crates/mesh/src/surface/topology.rs

@@ -236,7 +236,7 @@ fn solid_angle(s: &Surface, t: usize, x: [f64; 3]) -> f64 {
     let [a, b, c] = s.tri_xyz(t).map(|p| sub(p, x));
     let (la, lb, lc) = (norm(a), norm(b), norm(c));
     let den = la * lb * lc + dot(a, b) * lc + dot(a, c) * lb + dot(b, c) * la;
-    2.0 * det(a, b, c).atan2(den)
+    2.0 * crate::vec3::atan2(det(a, b, c), den)
 }
 
 /// Edges between different B-rep faces, and edges inside a face whose normals turn by more than