Blends: an edge is smooth when its faces' normals agree all along it, decided exactly

This replaces the five fixed stations of 0a166bf3d, which compared the faces' normals at the
edge's ends and quarters: a sampled guess that misses a crease or a tangency between stations.

Along each use's pcurve (on the edge's own parameter), the face's natural normal Su x Sv is a
polynomial in t on every piece where the pcurve stays in one span of its knots and one knot cell
of its surface: of degree (2p + 2q - 2) r for a B-spline surface of degree (p, q) and a pcurve
of degree r, (4p + 4q - 2) r for a rational one (its homogeneous form), 0 for a plane. The cross
product of the two outward normals is then a polynomial of degree D, the sum, and it vanishes
on the piece iff it vanishes at D + 1 points of it (the piece's Chebyshev points). Pieces end at
the pcurves' knots and where a pcurve crosses a knot line of its surface (periodic copies too),
each a root of U(t) - k W(t) on a Bézier span, isolated on its Bernstein coefficients. A point
where the normals differ is a witness of a crease; agreement at all of them is a proof of
tangency (crate::tangent_edge). Outward normals opposite (a cusp, a knife edge) is not smooth,
as SMLib's GetBlendEdgeuse has it since fork 81f411d; a seam, whose two uses are on one face,
is smooth.

The C API has no edge classification query; SMLib's own (GetBlendEdgeuse at t = 0.34,
IsTangentSector at five stations) is sampled too. A fork item to make it exact and expose it is
written up for the SMLib batch. The tangent-chain walk decides each edge once.

The bar whose front face runs tangent into a cylinder R8, its front top edge rounded r2: the
round's contact with the front face and the front face's tangent line with the wall are smooth;
the round's edge with the wall (tangent at (40, 0, 8) only) and a box edge are creases; the
knife edge of #851 is a cusp, not smooth (blend.rs test). The blend, fillet and chamfer kernel
tests and the fillets rig keep their outcomes.

Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>

crates/kernel/src/blend.rs

@@ -657,9 +657,10 @@ pub fn blend_edges(body: &Body, picks: &[Vec3]) -> Result<Vec<usize>> {
 /// tangentially (within 0.5° at a shared vertex), until no more do. Each pick takes its edge as a
 /// blend does (a pick at a vertex where a seam meets a rim takes the rim). An edge only one face
 /// uses ends a chain, and so do a seam and a smooth edge (its two faces meet tangentially, within
-/// 0.5°): they have no corner to round. A chain does not run on into an edge of the other
-/// convexity (a boss's concave base running on tangentially as the convex rim of its overhang):
-/// a round cannot pass from one side of the faces to the other, so it ends at that vertex.
+/// 0.5°, all along it): they have no corner to round. A chain does not run on into an edge of the
+/// other convexity (a boss's concave base running on tangentially as the convex rim of its
+/// overhang): a round cannot pass from one side of the faces to the other, so it ends at that
+/// vertex.
 /// Picks that find no edge are refused.
 pub fn tangent_chain(body: &Body, picks: &[Vec3]) -> Result<Vec<Vec3>> {
     if picks.is_empty() || picks.len() > 1000 || picks.iter().any(|p| !p.is_finite()) {
@@ -673,7 +674,7 @@ pub fn tangent_chain(body: &Body, picks: &[Vec3]) -> Result<Vec<Vec3>> {
             picked.push(e);
         }
     }
-    let chain = chain_of(&rec, &picked, &hidden);
+    let chain = chain_of(&rec, &picked, &hidden)?;
     let mut out = picks.to_vec();
     for e in chain.into_iter().skip(picked.len()) {
         out.push(edge_middle(&rec, e).ok_or_else(|| KernelError::Failed("a point on a chain edge".into()))?);
@@ -696,35 +697,34 @@ fn normal_reach(rec: &Faces) -> f64 {
     rec.rec.tolerance.max(1e-7) * 10.0
 }
 
-/// Stations along an edge, as fractions of its interval, where [`smooth_join`] compares its
-/// faces' normals: its ends, and between them.
-const SMOOTH_STATIONS: [f64; 5] = [0.0, 0.25, 0.5, 0.75, 1.0];
-
-/// Whether edge `e`'s two faces meet tangentially (within [`TANGENT_DEGREES`]) all along it: a
-/// smooth edge, with no corner to blend. An edge whose faces are tangent at one end only (a
-/// round running into a wall it is tangent to where the wall meets its face) is a crease
-/// everywhere else, which is not smooth.
-fn smooth_join(rec: &Faces, e: usize) -> bool {
-    let cos = TANGENT_DEGREES.to_radians().cos();
-    let accept = normal_reach(rec);
-    let faces: Vec<u32> = rec.faces_of(e).into_iter().collect();
-    let ([a, b], Some(edge)) = (faces.as_slice(), rec.rec.edges.get(e)) else { return false };
-    let Ok(curve) = CurveEval::new(edge.parameters()) else { return false };
-    let [t0, t1] = edge.interval;
-    SMOOTH_STATIONS.iter().all(|s| {
-        let Some(p) = curve.point(t0 + (t1 - t0) * s) else { return false };
-        match (rec.normal_at(*a as usize, p, accept), rec.normal_at(*b as usize, p, accept)) {
-            (Some(na), Some(nb)) => na.dot(nb) > cos,
-            _ => false,
-        }
-    })
+/// Whether edge `e`'s two uses meet tangentially (within [`TANGENT_DEGREES`]) all along it: a
+/// smooth edge, with no corner to blend (a seam a face meets itself across is one). An edge
+/// whose faces are tangent at one end only (a round running into a wall it is tangent to where
+/// the wall meets its face) is a crease everywhere else, which is not smooth; nor is a cusp, its
+/// faces folded back on each other (outward normals opposite), which SMLib rounds like any other
+/// edge. Decided exactly from the faces' surfaces and the uses' pcurves
+/// ([`crate::tangent_edge`]).
+fn smooth_join(rec: &Faces, e: usize) -> Result<bool> {
+    crate::tangent_edge::tangent_all_along(&rec.rec, |f| rec.side(f), e, TANGENT_DEGREES.to_radians())
 }
 
 /// `seeds` (edges in the record's order, each once), then every other edge that continues one
 /// of them tangentially, as [`tangent_chain`] follows them; `hidden` marks the edges a pick
 /// passes over.
-fn chain_of(rec: &Faces, seeds: &[usize], hidden: &[bool]) -> Vec<usize> {
+fn chain_of(rec: &Faces, seeds: &[usize], hidden: &[bool]) -> Result<Vec<usize>> {
     let mut chain = seeds.to_vec();
+    // Each edge's smoothness, decided once.
+    let mut smooth: Vec<Option<bool>> = vec![None; rec.rec.edges.len()];
+    let mut smooth_at = |f: usize| -> Result<bool> {
+        if let Some(known) = smooth.get(f).copied().flatten() {
+            return Ok(known);
+        }
+        let s = smooth_join(rec, f)?;
+        if let Some(slot) = smooth.get_mut(f) {
+            *slot = Some(s);
+        }
+        Ok(s)
+    };
     // The unit direction edge `e` leaves its end `v` in.
     let leaving = |e: usize, v: u32| -> Option<Vec3> {
         let [a, b] = rec.rec.edges.get(e)?.vertices;
@@ -761,7 +761,7 @@ fn chain_of(rec: &Faces, seeds: &[usize], hidden: &[bool]) -> Vec<usize> {
             let Some(d) = leaving(e, v) else { continue };
             for (f, other) in rec.rec.edges.iter().enumerate() {
                 let passed_over = hidden.get(f).copied().unwrap_or(false);
-                if passed_over || chain.contains(&f) || other.uses.len() < 2 || !other.vertices.contains(&v) || smooth_join(rec, f) {
+                if passed_over || chain.contains(&f) || other.uses.len() < 2 || !other.vertices.contains(&v) || smooth_at(f)? {
                     continue;
                 }
                 let turns_over = matches!((convex(e, v), convex(f, v)), (Some(a), Some(b)) if a != b);
@@ -771,7 +771,7 @@ fn chain_of(rec: &Faces, seeds: &[usize], hidden: &[bool]) -> Vec<usize> {
             }
         }
     }
-    chain
+    Ok(chain)
 }
 
 /// Before blending exactly `edges` (the record's edges `picks` took, in order): refuses an edge
@@ -783,7 +783,7 @@ pub(crate) fn unfinished_chain(body: &Body, picks: &[Vec3], edges: &[usize], hid
     let rec = Faces::of(body)?;
     let at = |i: usize| picks.get(i).map(|p| format!("({:.3}, {:.3}, {:.3})", p.x, p.y, p.z)).unwrap_or_default();
     for (i, &e) in edges.iter().enumerate() {
-        if smooth_join(&rec, e) {
+        if smooth_join(&rec, e)? {
             return Err(KernelError::Invalid(format!("the edge at {} is smooth: its faces meet tangentially, so there is nothing to round", at(i))));
         }
     }
@@ -793,7 +793,7 @@ pub(crate) fn unfinished_chain(body: &Body, picks: &[Vec3], edges: &[usize], hid
             seeds.push(e);
         }
     }
-    Ok(chain_of(&rec, &seeds, hidden).get(seeds.len()).and_then(|&next| edge_middle(&rec, next)))
+    Ok(chain_of(&rec, &seeds, hidden)?.get(seeds.len()).and_then(|&next| edge_middle(&rec, next)))
 }
 
 /// A blend's failure where a picked edge continues tangentially into `next`, an edge left out
@@ -857,6 +857,33 @@ mod tests {
         assert_eq!(tangent_chain(&rounded, &[side]).unwrap(), vec![side], "the chain ran into the round's seam");
     }
 
+    /// Smoothness is decided all along an edge. A bar 40 × 16 × 10 whose front face runs tangent
+    /// into a cylinder R8 at its far end, its front top edge rounded r2: the round's contact with
+    /// the front face and the front face's tangent line with the wall are smooth (plane against
+    /// cylinder, tangent all along); the round's edge with the wall is tangent at its lower end
+    /// (40, 0, 8) only, a crease everywhere else; a box edge is a crease. A knife edge (a wall of
+    /// radius 21 running out into a plane, #851: its faces tangent all along but folded back on
+    /// each other) is a cusp, which is not smooth.
+    #[test]
+    fn an_edge_is_smooth_only_where_its_faces_are_tangent_all_along_it() {
+        let bar = crate::build::box_solid(Vec3::ZERO, Vec3::new(40.0, 16.0, 10.0)).unwrap();
+        let post = crate::build::cylinder(Vec3::new(40.0, 8.0, 0.0), Vec3::Z, 8.0, 20.0).unwrap();
+        let peg = crate::boolean(&bar, &post, crate::BoolOp::Union).unwrap().unwrap();
+        let rounded = fillet(&peg, &[Vec3::new(20.0, 0.0, 10.0)], 2.0).unwrap();
+        let rec = Faces::of(&rounded).unwrap();
+        let smooth = |p: Vec3| smooth_join(&rec, blend_edges(&rounded, &[p]).unwrap()[0]).unwrap();
+        assert!(smooth(Vec3::new(20.0, 0.0, 8.0)), "the round's contact with the front face");
+        assert!(smooth(Vec3::new(40.0, 0.0, 4.0)), "the front face's tangent line with the wall");
+        assert!(!smooth(Vec3::new(37.2641, 0.4827, 9.3025)), "the round's edge with the wall");
+        assert!(!smooth(Vec3::new(20.0, 0.0, 0.0)), "the bar's front bottom edge");
+        let tip = 62.0 + 1.0 / 6.0;
+        let slab = crate::build::box_solid(Vec3::new(0.0, -21.0, 0.0), Vec3::new(tip, 0.0, 10.0)).unwrap();
+        let disc = crate::build::cylinder(Vec3::new(tip, -21.0, 0.0), Vec3::Z, 21.0, 10.0).unwrap();
+        let knife = crate::boolean(&slab, &disc, crate::BoolOp::Cut).unwrap().unwrap();
+        let rec = Faces::of(&knife).unwrap();
+        assert!(!smooth_join(&rec, blend_edges(&knife, &[Vec3::new(tip, 0.0, 5.0)]).unwrap()[0]).unwrap(), "the knife edge");
+    }
+
     #[test]
     fn g2_fillet_has_the_quintics_area() {
         let body = crate::build::box_solid(Vec3::ZERO, Vec3::new(20.0, 20.0, 20.0)).unwrap();

crates/kernel/src/lib.rs

@@ -139,6 +139,7 @@ mod surfaces;
 mod sweep_path;
 #[cfg(test)]
 mod sweep_rails_tests;
+mod tangent_edge;
 mod tess;
 #[cfg(test)]
 mod thin_face_tests;

crates/kernel/src/tangent_edge.rs

@@ -0,0 +1,349 @@
+//! Whether the two faces on either side of an edge meet tangentially all along it, decided
+//! exactly from the record.
+//!
+//! Each use of the edge carries its pcurve `(u(t), v(t))` in its face's parameters
+//! ([`Face::parameters`]: a plane, or the B-spline surface SMLib holds the face as), on the edge's
+//! own parameter `t`, so the two uses' points at one `t` are one point of the edge. The face's
+//! natural normal there is `S_u × S_v` at the pcurve's point. On a piece of `t` where the pcurve
+//! stays within one span of its own knots and one knot cell of its surface, both are
+//! polynomials: `S` of degree `(p, q)` in `(u, v)` (rational: its homogeneous form, of the same
+//! degree), the pcurve of degree `r` in `t`. So `S_u × S_v`, scaled by a positive power of the
+//! weights, is a polynomial in `t` of degree at most `(2p + 2q - 2) r` (rational: `(4p + 4q - 2) r`,
+//! from `(A_u w - A w_u) × (A_v w - A w_v)`), and a plane's is constant. The cross product of the
+//! two faces' normals is then a polynomial vector of degree at most the sum `D` of theirs, and a
+//! polynomial of degree `D` that vanishes at `D + 1` points vanishes everywhere: the faces are
+//! parallel all along the piece exactly when they are at `D + 1` of its points, its Chebyshev
+//! points, and on one side (smooth, not folded back on each other) when their outward normals
+//! point the same way there. No point is a guess at where the faces part: one where they do is
+//! a witness of a crease (or a cusp, the faces folded back), and agreement at all of them is
+//! a proof of tangency, within the angle the caller allows at each (between them, the cross
+//! product is bounded by its largest value at them times the Lebesgue constant of `D + 1`
+//! Chebyshev points, under `1 + (2/π) ln(D + 1)`).
+//!
+//! The pieces' ends are the pcurves' knots and the parameters where a pcurve crosses a knot line
+//! of its surface (periodic copies included), each a root of `U(t) - k W(t)` on one Bézier span
+//! of the pcurve, isolated by subdividing its Bernstein coefficients (no sign change: no root;
+//! one: exactly one) and found to rounding by bisection.
+
+use solidcraft_geom::Vec3;
+use solidcraft_geom::brep::eval::SurfaceEval;
+use solidcraft_geom::brep::{Brep, Curve2, Face, Surface};
+use solidcraft_geom::bspline::{HPoint, segment};
+
+use crate::{KernelError, Result};
+
+/// Most points one edge's decision evaluates.
+const MAX_POINTS: usize = 1_000_000;
+/// Largest degree of a piece's cross product the decision takes.
+const MAX_DEGREE: usize = 4096;
+/// Subdivisions of one Bernstein polynomial while isolating its roots: its pieces are then
+/// 2⁻⁴⁸ of a pcurve span, where a root left unisolated (a pcurve touching a knot line) ends a
+/// piece there.
+const MAX_DEPTH: u32 = 48;
+/// Subdivisions one root isolation may make: a root that touches its knot line (a double root,
+/// the only kind that keeps two sign changes as the pieces shrink) costs one a level.
+const MAX_SUBDIVISIONS: usize = 4096;
+/// Below this sine of the angle between `S_u` and `S_v` a face's normal is undefined at that
+/// point (a pole): the normal, and so the cross product, is zero there.
+const DEGENERATE: f64 = 1e-9;
+
+/// One use of the edge: its face's surface, the sign turning its natural normal outward, its
+/// pcurve and the degree of its normal along a piece.
+struct Use<'a> {
+    surface: SurfaceEval,
+    sign: f64,
+    pcurve: &'a Curve2,
+    degree: usize,
+    /// The knot lines of the surface in u and v (its distinct knots over its domain) and its
+    /// periods: a polynomial cell's sides.
+    lines: [Vec<f64>; 2],
+    period: [Option<f64>; 2],
+}
+
+fn failed(msg: String) -> KernelError {
+    KernelError::Failed(msg)
+}
+
+/// A B-spline's degree from its knot and control point counts.
+fn degree_of(knots: usize, count: usize) -> Option<usize> {
+    knots.checked_sub(count.checked_add(1)?)
+}
+
+/// The distinct knots of a knot vector over its domain `[knots[p], knots[n]]`.
+fn lines(knots: &[f64], count: usize) -> Option<Vec<f64>> {
+    let p = degree_of(knots.len(), count)?;
+    let mut out: Vec<f64> = knots.get(p..=count)?.to_vec();
+    out.dedup();
+    Some(out)
+}
+
+impl<'a> Use<'a> {
+    fn new(face: &Face, sign: f64, pcurve: &'a Curve2) -> Result<Use<'a>> {
+        let surface = face.parameters();
+        let r = match pcurve {
+            Curve2::Spline { knots, points } => degree_of(knots.len(), points.len()),
+            Curve2::RationalSpline { knots, homogeneous } => degree_of(knots.len(), homogeneous.len()),
+        }
+        .ok_or_else(|| failed("a pcurve's degree".into()))?;
+        // The normal's degree in (u, v) and the knot lines.
+        let (uv_degree, lines) = match surface {
+            Surface::Plane { .. } | Surface::PlaneThrough { .. } => (0, [Vec::new(), Vec::new()]),
+            Surface::Spline { u_knots, v_knots, points } => {
+                let (nu, nv) = (points.len(), points.first().map_or(0, Vec::len));
+                let (p, q) = (degree_of(u_knots.len(), nu), degree_of(v_knots.len(), nv));
+                let (Some(p), Some(q), Some(lu), Some(lv)) = (p, q, lines(u_knots, nu), lines(v_knots, nv)) else {
+                    return Err(failed("a spline face's degree".into()));
+                };
+                ((2 * p + 2 * q).saturating_sub(2), [lu, lv])
+            }
+            Surface::RationalSpline { u_knots, v_knots, homogeneous } => {
+                let (nu, nv) = (homogeneous.len(), homogeneous.first().map_or(0, Vec::len));
+                let (p, q) = (degree_of(u_knots.len(), nu), degree_of(v_knots.len(), nv));
+                let (Some(p), Some(q), Some(lu), Some(lv)) = (p, q, lines(u_knots, nu), lines(v_knots, nv)) else {
+                    return Err(failed("a spline face's degree".into()));
+                };
+                ((4 * p + 4 * q).saturating_sub(2), [lu, lv])
+            }
+            // SMLib gives every curved analytic face its NURBS form, which its pcurves are on.
+            _ => return Err(crate::blend::unsupported("deciding tangency along an edge of a face that is no plane and has no NURBS form")),
+        };
+        let eval = SurfaceEval::new(surface).map_err(|e| failed(format!("a face's surface: {e}")))?;
+        let period = eval.domain().period;
+        Ok(Use { surface: eval, sign, pcurve, degree: uv_degree.saturating_mul(r), lines, period })
+    }
+
+    /// The outward unit normal at the pcurve's point at `t`; `Ok(None)` where it is undefined (a
+    /// pole: a zero of the normal).
+    fn normal(&self, t: f64) -> Result<Option<Vec3>> {
+        let uv = self.pcurve.eval(t).ok_or_else(|| failed(format!("a pcurve at {t}")))?;
+        let (_, su, sv) = self.surface.derivs(uv[0], uv[1]).ok_or_else(|| failed(format!("a face at ({}, {})", uv[0], uv[1])))?;
+        let n = su.cross(sv);
+        if !(n.len() > DEGENERATE * su.len() * sv.len()) {
+            return Ok(None);
+        }
+        Ok(n.normalized().map(|n| n * self.sign))
+    }
+
+    /// The parameters in `(t0, t1)` where the normal stops being one polynomial: the pcurve's
+    /// own knots, and where it crosses a knot line of the surface.
+    fn breaks(&self, [t0, t1]: [f64; 2], out: &mut Vec<f64>) -> Result<()> {
+        let (knots, ctrl): (&[f64], Vec<HPoint>) = match self.pcurve {
+            Curve2::Spline { knots, points } => (knots, points.iter().map(|p| HPoint { x: p[0], y: p[1], z: 0.0, w: 1.0 }).collect()),
+            Curve2::RationalSpline { knots, homogeneous } => {
+                (knots, homogeneous.iter().map(|h| HPoint { x: h[0], y: h[1], z: 0.0, w: h[2] }).collect())
+            }
+        };
+        let r = degree_of(knots.len(), ctrl.len()).ok_or_else(|| failed("a pcurve's degree".into()))?;
+        let mut spans: Vec<f64> = vec![t0];
+        spans.extend(knots.iter().copied().filter(|k| *k > t0 && *k < t1));
+        spans.push(t1);
+        spans.dedup();
+        out.extend(spans.iter().copied());
+        if self.lines.iter().all(Vec::is_empty) {
+            return Ok(());
+        }
+        for w in spans.windows(2) {
+            let &[a, b] = w else { continue };
+            let (_, bez) = segment(r, knots, &ctrl, a, b).map_err(|e| failed(format!("a pcurve span: {e}")))?;
+            if bez.len() != r + 1 {
+                return Err(failed("a pcurve span is no Bézier piece".into()));
+            }
+            for d in 0..2 {
+                let coord = |h: &HPoint| if d == 0 { h.x } else { h.y };
+                // The span's range in this direction: its control points' hull (weights > 0).
+                let along: Vec<f64> = bez.iter().map(|h| coord(h) / h.w).collect();
+                let lo = along.iter().copied().fold(f64::INFINITY, f64::min);
+                let hi = along.iter().copied().fold(f64::NEG_INFINITY, f64::max);
+                if !(lo.is_finite() && hi.is_finite()) {
+                    return Err(failed("a pcurve's control points".into()));
+                }
+                for k in cell_sides(&self.lines[d], self.period[d], lo, hi)? {
+                    // Coefficients within rounding of zero are zero: a pcurve running along
+                    // the knot line (a seam, a boundary) crosses none there.
+                    let scale = bez.iter().map(|h| coord(h).abs() + (k * h.w).abs()).fold(0.0, f64::max);
+                    let floor = 64.0 * f64::EPSILON * scale;
+                    let c: Vec<f64> = bez.iter().map(|h| coord(h) - k * h.w).map(|x| if x.abs() <= floor { 0.0 } else { x }).collect();
+                    if c.iter().all(|x| *x == 0.0) {
+                        continue;
+                    }
+                    let mut roots = Vec::new();
+                    let mut work = MAX_SUBDIVISIONS;
+                    bernstein_roots(&c, 0.0, 1.0, MAX_DEPTH, &mut work, &mut roots);
+                    out.extend(roots.into_iter().map(|s| a + (b - a) * s));
+                }
+            }
+        }
+        Ok(())
+    }
+}
+
+/// The knot lines `lines` (and, with a period, their copies a whole number of periods away) that
+/// lie in `[lo, hi]`.
+fn cell_sides(lines: &[f64], period: Option<f64>, lo: f64, hi: f64) -> Result<Vec<f64>> {
+    let Some(p) = period.filter(|p| *p > 0.0) else {
+        return Ok(lines.iter().copied().filter(|k| *k >= lo && *k <= hi).collect());
+    };
+    let (Some(&first), Some(&last)) = (lines.first(), lines.last()) else { return Ok(Vec::new()) };
+    let (n0, n1) = (((lo - last) / p).floor(), ((hi - first) / p).ceil());
+    if !(n1 - n0 <= 64.0) {
+        return Err(failed("a pcurve running round its surface too many times".into()));
+    }
+    let mut out = Vec::new();
+    let mut n = n0;
+    while n <= n1 {
+        out.extend(lines.iter().map(|k| k + n * p).filter(|k| *k >= lo && *k <= hi));
+        n += 1.0;
+    }
+    Ok(out)
+}
+
+/// The value at `s` of the polynomial with Bernstein coefficients `c` on `[0, 1]` (de Casteljau).
+fn bernstein_value(c: &[f64], s: f64) -> f64 {
+    let mut b = c.to_vec();
+    while b.len() > 1 {
+        b = b.windows(2).map(|w| if let [x, y] = w { x * (1.0 - s) + y * s } else { 0.0 }).collect();
+    }
+    b.first().copied().unwrap_or(0.0)
+}
+
+/// The Bernstein coefficients of the polynomial `c` on `[0, ½]` and `[½, 1]`.
+fn bernstein_halves(c: &[f64]) -> (Vec<f64>, Vec<f64>) {
+    let mut left = Vec::with_capacity(c.len());
+    let mut right = Vec::with_capacity(c.len());
+    let mut b = c.to_vec();
+    while let (Some(&f), Some(&l)) = (b.first(), b.last()) {
+        left.push(f);
+        right.push(l);
+        b = b.windows(2).map(|w| if let [x, y] = w { 0.5 * (x + y) } else { 0.0 }).collect();
+    }
+    right.reverse();
+    (left, right)
+}
+
+/// The roots in `[lo, hi]` of the polynomial with Bernstein coefficients `c` there, mapped
+/// into `[lo, hi]`'s place in `[0, 1]`: none where the coefficients keep one sign, exactly one
+/// where they change sign once (found by bisection to rounding), otherwise those of each half.
+fn bernstein_roots(c: &[f64], lo: f64, hi: f64, depth: u32, work: &mut usize, out: &mut Vec<f64>) {
+    let (Some(&first), Some(&last)) = (c.first(), c.last()) else { return };
+    if first == 0.0 {
+        out.push(lo);
+    }
+    if last == 0.0 {
+        out.push(hi);
+    }
+    let signs: Vec<bool> = c.iter().filter(|x| **x != 0.0).map(|x| *x > 0.0).collect();
+    let changes = signs.windows(2).filter(|w| matches!(w, [x, y] if x != y)).count();
+    if changes == 0 {
+        return;
+    }
+    if changes == 1 && first != 0.0 && last != 0.0 {
+        let (mut a, mut b) = (0.0, 1.0);
+        for _ in 0..64 {
+            let m = 0.5 * (a + b);
+            if !(m > a && m < b) {
+                break;
+            }
+            if (bernstein_value(c, m) > 0.0) == (first > 0.0) {
+                a = m;
+            } else {
+                b = m;
+            }
+        }
+        out.push(lo + (hi - lo) * 0.5 * (a + b));
+        return;
+    }
+    let mid = 0.5 * (lo + hi);
+    // Out of depth or work, the roots are within this piece: it ends a piece of the edge at
+    // both its ends.
+    if depth == 0 || *work == 0 {
+        out.extend([lo, hi]);
+        return;
+    }
+    *work -= 1;
+    let (l, r) = bernstein_halves(c);
+    bernstein_roots(&l, lo, mid, depth - 1, work, out);
+    bernstein_roots(&r, mid, hi, depth - 1, work, out);
+}
+
+/// Whether edge `e`'s two uses' faces meet with outward normals within `max_angle` (radians) of
+/// each other all along it, `sign(f)` turning face `f`'s natural normal outward. An edge with
+/// other than two uses has no such pair: false. See the module documentation.
+pub(crate) fn tangent_all_along(rec: &Brep, sign: impl Fn(usize) -> Option<f64>, e: usize, max_angle: f64) -> Result<bool> {
+    let edge = rec.edges.get(e).ok_or_else(|| failed(format!("edge {e}")))?;
+    let [ua, ub] = edge.uses.as_slice() else { return Ok(false) };
+    let use_of = |u: &solidcraft_geom::brep::EdgeUseRef| -> Result<Use<'_>> {
+        let face = rec.faces.get(u.face as usize).ok_or_else(|| failed(format!("face {}", u.face)))?;
+        let oriented = match face.loops.get(u.r#loop as usize) {
+            Some(solidcraft_geom::brep::Loop::Edges { edges }) => edges.get(u.position as usize),
+            _ => None,
+        }
+        .ok_or_else(|| failed(format!("a use of edge {e}")))?;
+        let s = sign(u.face as usize).ok_or_else(|| failed(format!("face {} is in no shell", u.face)))?;
+        Use::new(face, s, &oriented.pcurve)
+    };
+    let (a, b) = (use_of(ua)?, use_of(ub)?);
+    let degree = a.degree.saturating_add(b.degree);
+    if degree > MAX_DEGREE {
+        return Err(failed(format!("deciding tangency along edge {e}: its faces' normals are of degree {degree}")));
+    }
+    let interval = edge.interval;
+    let mut breaks = Vec::new();
+    a.breaks(interval, &mut breaks)?;
+    b.breaks(interval, &mut breaks)?;
+    breaks.retain(|t| *t >= interval[0] && *t <= interval[1]);
+    breaks.sort_by(f64::total_cmp);
+    breaks.dedup();
+    let nodes = degree + 1;
+    if breaks.len().saturating_mul(nodes) > MAX_POINTS {
+        return Err(failed(format!("deciding tangency along edge {e}: {} pieces of degree {degree}", breaks.len())));
+    }
+    let cos = max_angle.cos();
+    let mut defined = false;
+    for w in breaks.windows(2) {
+        let &[t0, t1] = w else { continue };
+        let (m, h) = (0.5 * (t0 + t1), 0.5 * (t1 - t0));
+        if !(h > 0.0) {
+            continue;
+        }
+        for j in 0..nodes {
+            let t = m + h * ((2 * j + 1) as f64 * std::f64::consts::PI / (2 * nodes) as f64).cos();
+            if let (Some(na), Some(nb)) = (a.normal(t)?, b.normal(t)?) {
+                if !(na.dot(nb) >= cos) {
+                    return Ok(false);
+                }
+                defined = true;
+            }
+        }
+    }
+    if !defined {
+        return Err(failed(format!("edge {e}: its faces' normals are undefined all along it")));
+    }
+    Ok(true)
+}
+
+#[cfg(test)]
+mod tests {
+    use super::*;
+
+    #[test]
+    fn bernstein_roots_are_isolated_and_exact() {
+        // (s - 0.25)(s - 0.75) = s² - s + 3/16 in Bernstein form on [0, 1]: 3/16, 3/16 - 1/2, 3/16.
+        let c = [0.1875, 0.1875 - 0.5, 0.1875];
+        let mut roots = Vec::new();
+        let mut work = MAX_SUBDIVISIONS;
+        bernstein_roots(&c, 0.0, 1.0, MAX_DEPTH, &mut work, &mut roots);
+        roots.sort_by(f64::total_cmp);
+        assert_eq!(roots.len(), 2, "{roots:?}");
+        assert!((roots[0] - 0.25).abs() < 1e-15 && (roots[1] - 0.75).abs() < 1e-15, "{roots:?}");
+        // No sign change: no root.
+        let mut none = Vec::new();
+        bernstein_roots(&[1.0, 0.5, 2.0], 0.0, 1.0, MAX_DEPTH, &mut work, &mut none);
+        assert!(none.is_empty());
+    }
+
+    #[test]
+    fn periodic_knot_lines_repeat_by_whole_periods() {
+        let sides = cell_sides(&[0.0, 1.0, 2.0], Some(2.0), 1.5, 4.5).unwrap();
+        assert_eq!(sides, vec![2.0, 2.0, 3.0, 4.0, 4.0]);
+    }
+}