geom::krawczyk: the one Krawczyk operator, the certificate's (sm_GnEncloseOver) order, extracted (#15140)

The fillet certificate's feet enclosure forms Krawczyk's contraction C = I - Y F' with its q, and K = mid - Y G + C (X - mid), inline at its two sites (the step and the tightening). #15140 wants one operator in the crate, each caller keeping its own strategy (centre, first box, preconditioner, widening). This moves the two computations into geom::krawczyk as contraction() and k_operator(), in the certificate's own order (the fork's): C_ij from delta_ij, subtracting y_ik F'_kj in k's order, q a running rounded-up row sum; K_i from mid_i, subtracting y_ij g_j, then adding c_ij (x_j - mid_j), in j's order. The certificate calls them; q < 1 stays in its own test. Same expressions in the same order, so the certificate's boxes are bit for bit; the dead f_at helper goes.

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

solidcore/crates/solidcore/src/blend/fillet_section/certificate.rs

@@ -161,6 +161,7 @@ use super::tay::{
 };
 use super::taylor_model::Tm;
 use crate::error::{Result, SolidError};
+use crate::geom::krawczyk::{contraction, k_operator};
 use crate::kernel::Kernel;
 use crate::nurbs::b_spline_curve::BSplineCurve;
 use crate::surface::surface::Surface;
@@ -808,10 +809,6 @@ fn ivl_at(x: &[Ivl; 4], i: usize) -> Ivl {
     x.get(i).copied().unwrap_or(IZ)
 }
 
-fn f_at(x: &[f64; 4], i: usize) -> f64 {
-    x.get(i).copied().unwrap_or(f64::NAN)
-}
-
 /// The system over the partials `p` (`sSys` with `apP` at `asP`).
 fn system<'a>(
     p: &'a [Partials<Ivl>; 2],
@@ -1961,30 +1958,8 @@ fn enclose_over_cells(
         );
         let jm = jac.map(|row| row.map(Ivl::mid));
         y = gn_try!(inverse4(&jm), Inner::Unproved);
-        q = 0.0;
-        for (i, (yi, ci)) in y.iter().zip(cm.iter_mut()).enumerate() {
-            let mut row = 0.0_f64;
-            for (j, cij) in ci.iter_mut().enumerate() {
-                let mut c = Ivl::from(if i == j { 1.0 } else { 0.0 });
-                for (&yik, jk) in yi.iter().zip(&jac) {
-                    c = c - Ivl::from(yik) * jk.get(j).copied().unwrap_or(IZ);
-                }
-                *cij = c;
-                row = round_up(row + c.mag());
-            }
-            q = cmax(q, row);
-        }
-        let mut kk = [IZ; 4];
-        for (i, ((yi, ci), ki)) in y.iter().zip(&cm).zip(kk.iter_mut()).enumerate() {
-            let mut k = Ivl::from(f_at(&xm, i));
-            for (&yij, gj) in yi.iter().zip(&sg) {
-                k = k - Ivl::from(yij) * gj.at(0, 0);
-            }
-            for (j, &cij) in ci.iter().enumerate() {
-                k = k + cij * (ivl_at(&bx_, j) - Ivl::from(f_at(&xm, j)));
-            }
-            *ki = k;
-        }
+        (cm, q) = contraction(&y, &jac);
+        let kk = k_operator(&xm, &y, &sg.each_ref().map(|g| g.at(0, 0)), &cm, &bx_);
         let verdict = classify_k(&kk, &bx_);
         let inside = q < 1.0 && verdict == KVerdict::Inside;
         if inside {
@@ -2081,34 +2056,15 @@ fn enclose_over_cells(
         let Ok(yc) = inverse4(&jac.map(|row| row.map(Ivl::mid))) else {
             break;
         };
-        let mut cc = [[IZ; 4]; 4];
-        let mut qc = 0.0_f64;
-        for (i, (yi, ci)) in yc.iter().zip(cc.iter_mut()).enumerate() {
-            let mut row = 0.0_f64;
-            for (j, cij) in ci.iter_mut().enumerate() {
-                let mut c = Ivl::from(if i == j { 1.0 } else { 0.0 });
-                for (&yik, jk) in yi.iter().zip(&jac) {
-                    c = c - Ivl::from(yik) * jk.get(j).copied().unwrap_or(IZ);
-                }
-                *cij = c;
-                row = round_up(row + c.mag());
-            }
-            qc = cmax(qc, row);
-        }
+        let (cc, qc) = contraction(&yc, &jac);
         let contracts = qc < 1.0;
         if !contracts {
             break;
         }
         let mut next = bx_;
         let mut empty = false;
-        for (i, ((yi, ci), xi)) in yc.iter().zip(&cc).zip(next.iter_mut()).enumerate() {
-            let mut kv = Ivl::from(f_at(&xc, i));
-            for (&yij, gj) in yi.iter().zip(&sgc) {
-                kv = kv - Ivl::from(yij) * gj.at(0, 0);
-            }
-            for (j, &cij) in ci.iter().enumerate() {
-                kv = kv + cij * (ivl_at(&bx_, j) - Ivl::from(f_at(&xc, j)));
-            }
+        let kc = k_operator(&xc, &yc, &sgc.each_ref().map(|g| g.at(0, 0)), &cc, &bx_);
+        for (xi, kv) in next.iter_mut().zip(kc) {
             let lo = cmax(xi.lo, kv.lo);
             let hi = cmin(xi.hi, kv.hi);
             let holds = lo <= hi;

solidcore/crates/solidcore/src/geom/krawczyk.rs

@@ -17,14 +17,18 @@
 //! |------|----------|
 //! | [`krawczyk`], [`newton`], [`KRAWCZYK_STEPS`], [`NEWTON_STEPS`] | (not upstream: #15009's, moved) |
 //! | [`inverse3`], [`inverse_of_columns2`], [`inverse_of_columns3`], [`inverse_of_columns4`] | (not upstream: #15009's, moved) |
+//! | [`contraction`], [`k_operator`] | `sm_GnEncloseOver`'s `C`, `q` and `K` (the fork's order), extracted (#15140) |
 //!
 //! Deviations: none; the bodies are #15009's as they were in `reach/net.rs`, bit for bit.
 
 use crate::blend::fillet_section::certificate::inverse4;
-use crate::blend::fillet_section::interval::{round_up, Ivl};
+use crate::blend::fillet_section::interval::{cmax, round_up, Ivl};
 use crate::error::Result;
 use crate::geom::vector3d::Vec3;
 
+/// The point interval 0 (a missing entry's value, as `IZ` in the certificate).
+const ZERO: Ivl = Ivl::new(0.0, 0.0);
+
 /// The most Krawczyk steps one root gets ([`krawczyk`]).
 pub(crate) const KRAWCZYK_STEPS: usize = 8;
 
@@ -109,6 +113,60 @@ pub(crate) fn krawczyk<const N: usize>(
     Ok(None)
 }
 
+/// The contraction of Krawczyk's operator for the preconditioner `y` (rows) and `F′` over a
+/// box `jac` (rows: `jac[k][j]` the `k`th equation's partial in the `j`th unknown):
+/// `C = I − Y F′`, each `C_ij` formed from `δ_ij` by subtracting `y_ik F′_kj` in `k`'s
+/// order, and `q = max_i Σ_j |C_ij|`, each row's sum rounded up as it grows. The fork's
+/// order (`sm_GnEncloseOver`), shared with [`k_operator`] (tangency-labs/solidcraft#15140).
+pub(crate) fn contraction<const N: usize>(
+    y: &[[f64; N]; N],
+    jac: &[[Ivl; N]; N],
+) -> ([[Ivl; N]; N], f64) {
+    let mut cm = [[Ivl::from(0.0); N]; N];
+    let mut q = 0.0_f64;
+    for (i, (yi, ci)) in y.iter().zip(cm.iter_mut()).enumerate() {
+        let mut row = 0.0_f64;
+        for (j, cij) in ci.iter_mut().enumerate() {
+            let mut c = Ivl::from(if i == j { 1.0 } else { 0.0 });
+            for (&yik, jk) in yi.iter().zip(jac) {
+                c = c - Ivl::from(yik) * jk.get(j).copied().unwrap_or(ZERO);
+            }
+            *cij = c;
+            row = round_up(row + c.mag());
+        }
+        q = cmax(q, row);
+    }
+    (cm, q)
+}
+
+/// Krawczyk's operator over the box `x` about `mid`: `K_i = mid_i − Σ_j y_ij g_j + Σ_j
+/// c_ij (x_j − mid_j)`, from `mid_i`, each term subtracted then added in `j`'s order (the
+/// fork's `sm_GnEncloseOver`), `g` the system at the point `mid`, `c` the [`contraction`] for
+/// `y`. A box `x` with `K` strictly inside it holds a root of the system (module docs); each
+/// caller keeps its own centre, first box, preconditioner and widening, and a caller whose
+/// proof needs `q < 1` tests it itself (tangency-labs/solidcraft#15140).
+pub(crate) fn k_operator<const N: usize>(
+    mid: &[f64; N],
+    y: &[[f64; N]; N],
+    g: &[Ivl; N],
+    c: &[[Ivl; N]; N],
+    x: &[Ivl; N],
+) -> [Ivl; N] {
+    let mut kk = [Ivl::from(0.0); N];
+    for (i, ((yi, ci), ki)) in y.iter().zip(c).zip(kk.iter_mut()).enumerate() {
+        let mut k = Ivl::from(mid.get(i).copied().unwrap_or(f64::NAN));
+        for (&yij, gj) in yi.iter().zip(g) {
+            k = k - Ivl::from(yij) * *gj;
+        }
+        for (j, &cij) in ci.iter().enumerate() {
+            let xj = x.get(j).copied().unwrap_or(ZERO);
+            k = k + cij * (xj - Ivl::from(mid.get(j).copied().unwrap_or(f64::NAN)));
+        }
+        *ki = k;
+    }
+    kk
+}
+
 /// `F` at a point and its columns there, as numbers ([`newton`]).
 pub(crate) type Linearised<const N: usize> = Option<([f64; N], [[f64; N]; N])>;