Release SolveCraft 0.1.4

Resolve hover before rendering and suppress it throughout navigation gestures and animated views. Closes #92.

ATTRIBUTION.md

@@ -49,3 +49,4 @@ Generated-in-code assets are original and have no file to list:
 | vendor/truck-meshalgo | truck-meshalgo 0.4.0, https://github.com/ricosjp/truck (RICOS Co. Ltd.) | Apache-2.0 | meshing: parameter searches kept inside the surface domain (vendor/README.md) |
 | vendor/truck-stepio | truck-stepio 0.3.0, https://github.com/ricosjp/truck (RICOS Co. Ltd.) | Apache-2.0 | STEP writer: intersection-curve surfaces written at unique ids (vendor/README.md) |
 | vendor/truck-shapeops | truck-shapeops 0.4.0, https://github.com/ricosjp/truck (RICOS Co. Ltd.) | Apache-2.0 | robust meshing and nearest-parameter fallback in booleans; trace macro (vendor/README.md) |
+| vendor/matext4cgmath | matext4cgmath 0.1.0, https://github.com/ricosjp/matext4cgmath (RICOS Co. Ltd.) | Apache-2.0 | root refinement stops after 100 Newton steps (vendor/README.md) |

Cargo.lock

@@ -1753,8 +1753,6 @@ checksum = "f9f8bd3e56ce4dfc153cf470fffbfa98c7620958b312ca5c3a4b8d5181fd13c6"
 [[package]]
 name = "matext4cgmath"
 version = "0.1.0"
-source = "registry+https://github.com/rust-lang/crates.io-index"
-checksum = "e6541e181de37f70f0aceb25441823a3d0efa9cc1d23f475a0d3926678949178"
 dependencies = [
  "cgmath",
  "katexit",
@@ -3059,7 +3057,7 @@ dependencies = [
 
 [[package]]
 name = "solvecraft"
-version = "0.1.3"
+version = "0.1.4"
 dependencies = [
  "eframe",
  "egui",
@@ -3073,7 +3071,7 @@ dependencies = [
 
 [[package]]
 name = "solvecraft-cli"
-version = "0.1.3"
+version = "0.1.4"
 dependencies = [
  "serde_json",
  "solvecraft-engine",
@@ -3082,7 +3080,7 @@ dependencies = [
 
 [[package]]
 name = "solvecraft-doc"
-version = "0.1.3"
+version = "0.1.4"
 dependencies = [
  "log",
  "serde",
@@ -3096,7 +3094,7 @@ dependencies = [
 
 [[package]]
 name = "solvecraft-engine"
-version = "0.1.3"
+version = "0.1.4"
 dependencies = [
  "log",
  "serde",
@@ -3112,7 +3110,7 @@ dependencies = [
 
 [[package]]
 name = "solvecraft-geom"
-version = "0.1.3"
+version = "0.1.4"
 dependencies = [
  "serde",
  "serde_json",
@@ -3120,7 +3118,7 @@ dependencies = [
 
 [[package]]
 name = "solvecraft-io"
-version = "0.1.3"
+version = "0.1.4"
 dependencies = [
  "quick-xml 0.38.4",
  "solvecraft-doc",
@@ -3132,7 +3130,7 @@ dependencies = [
 
 [[package]]
 name = "solvecraft-kernel"
-version = "0.1.3"
+version = "0.1.4"
 dependencies = [
  "log",
  "serde",
@@ -3147,7 +3145,7 @@ dependencies = [
 
 [[package]]
 name = "solvecraft-mcp"
-version = "0.1.3"
+version = "0.1.4"
 dependencies = [
  "serde_json",
  "solvecraft-engine",
@@ -3155,7 +3153,7 @@ dependencies = [
 
 [[package]]
 name = "solvecraft-render"
-version = "0.1.3"
+version = "0.1.4"
 dependencies = [
  "image",
  "serde",
@@ -3164,7 +3162,7 @@ dependencies = [
 
 [[package]]
 name = "solvecraft-sketch"
-version = "0.1.3"
+version = "0.1.4"
 dependencies = [
  "serde",
  "serde_json",
@@ -3175,7 +3173,7 @@ dependencies = [
 
 [[package]]
 name = "solvecraft-ui-egui"
-version = "0.1.3"
+version = "0.1.4"
 dependencies = [
  "egui",
  "egui-wgpu",
@@ -3189,7 +3187,7 @@ dependencies = [
 
 [[package]]
 name = "solvecraft-web"
-version = "0.1.3"
+version = "0.1.4"
 dependencies = [
  "eframe",
  "egui",
@@ -4488,7 +4486,7 @@ checksum = "e450f9b2ed1dff33c94c12589a87338689467b9c4f5d8a5710bd09a847d2c8a7"
 
 [[package]]
 name = "xtask"
-version = "0.1.3"
+version = "0.1.4"
 dependencies = [
  "serde_json",
 ]

Cargo.toml

@@ -5,7 +5,7 @@ members = ["crates/*", "apps/*", "xtask"]
 exclude = ["vendor"]
 
 [workspace.package]
-version = "0.1.3"
+version = "0.1.4"
 edition = "2024"
 license = "MIT OR Apache-2.0"
 rust-version = "1.90"
@@ -98,3 +98,4 @@ debug = 0
 truck-shapeops = { path = "vendor/truck-shapeops" }
 truck-stepio = { path = "vendor/truck-stepio" }
 truck-meshalgo = { path = "vendor/truck-meshalgo" }
+matext4cgmath = { path = "vendor/matext4cgmath" }

crates/doc/src/document.rs

@@ -125,7 +125,8 @@ pub struct Extent {
     /// Optional second side distance (two-sided extrude), along the negative normal.
     #[serde(default, skip_serializing_if = "Option::is_none")]
     pub distance2: Option<String>,
-    /// Start offset from the sketch plane (expression).
+    /// Start offset from the sketch plane (expression), along the sketch normal whatever the
+    /// direction.
     #[serde(default, skip_serializing_if = "Option::is_none")]
     pub start_offset: Option<String>,
     /// Through all bodies in the direction(s) instead of `distance`.

crates/engine/src/boolean_regression_tests.rs

@@ -0,0 +1,66 @@
+//! Booleans that failed, panicked or hung, replayed from the scripts in
+//! `examples/regressions/` (each script's description says what went wrong).
+
+use serde_json::{Value, json};
+
+use crate::Session;
+
+fn replay(script: &str) -> Session {
+    let mut s = Session::default();
+    let script: Value = serde_json::from_str(script).unwrap();
+    s.run_script(&script).unwrap();
+    s
+}
+
+fn volume(s: &Session, body: &str) -> f64 {
+    s.model.state().body(body).map(|b| solvecraft_kernel::measure(&b.body).unwrap().volume).unwrap_or(0.0)
+}
+
+/// Slab less a 557.68 × 100 × 120 pocket and the loft below its floor (556.68 × 99 down 4 mm to
+/// 80 × 80; prismatoid: h/6 · (A_top + A_bottom + 4 A_mid)).
+fn slab_pocket_loft() -> (f64, f64) {
+    let slab = 900.0 * 480.0 * 160.0;
+    let loft = 4.0 / 6.0 * (556.68 * 99.0 + 80.0 * 80.0 + 4.0 * (556.68 + 80.0) / 2.0 * (99.0 + 80.0) / 2.0);
+    (slab - 557.68 * 100.0 * 120.0, loft)
+}
+
+#[test]
+fn loft_cut_beside_a_round_hole() {
+    let s = replay(include_str!("../../../examples/regressions/loft_cut_beside_round_hole.json"));
+    let (pocketed, loft) = slab_pocket_loft();
+    let hole = std::f64::consts::PI * 17.5 * 17.5 * 25.0;
+    let want = pocketed - loft - hole;
+    let got = volume(&s, "Body1");
+    assert!((got - want).abs() < 1e-6 * want, "{got} vs {want}");
+}
+
+#[test]
+fn shallow_loft_cut_from_the_top_face() {
+    let s = replay(include_str!("../../../examples/regressions/shallow_loft_cut_from_top_face.json"));
+    let (_, loft) = slab_pocket_loft();
+    let want = 900.0 * 480.0 * 160.0 - loft;
+    let got = volume(&s, "Body1");
+    assert!((got - want).abs() < 1e-9 * want, "{got} vs {want}");
+}
+
+#[test]
+fn oblique_cut_through_a_small_arc_matches_a_split() {
+    let scripts = [
+        include_str!("../../../examples/regressions/oblique_cut_of_cylinder.json"),
+        include_str!("../../../examples/regressions/oblique_cut_of_cylinder_hang.json"),
+    ];
+    // The same bar split by the prism's tilted top: the piece above it is the cut's result.
+    let mut split = Session::default();
+    let script: Value = serde_json::from_str(scripts[0]).unwrap();
+    let steps = script.get("commands").and_then(Value::as_array).unwrap();
+    split.run_script(&Value::Array(steps[..5].to_vec())).unwrap();
+    let whole = volume(&split, "Tool");
+    split.execute("solid.split_body", &json!({"body": "Tool", "plane": {"origin": [0, -472, 853.91], "normal": [0, 15.44, 40]}})).unwrap();
+    let above = volume(&split, "Tool");
+    assert!(above < whole - 10.0, "{above} of {whole}");
+    for script in scripts {
+        let s = replay(script);
+        let got = volume(&s, "Tool");
+        assert!((got - above).abs() < 1e-6 * above, "{got} vs {above}");
+    }
+}

crates/engine/src/cmd/features.rs

@@ -13,7 +13,7 @@ pub static COMMANDS: &[CommandSpec] = &[
         .at("SOLID", "CREATE")
         .icon("extrude")
         .key("E")
-        .params("distance: expr (or through_all: true); taper?: angle expr; sketch?: id|name (default: active or last sketch); profiles?: all | [index] | [[curve ids]] | [{point:[x,y]}]; face?: [x,y,z] (extrude a planar body face instead of a sketch profile); direction?: positive|negative|symmetric; distance2?; start_offset?; to?: [x,y,z] (To Object: a face, vertex or point; to_offset?: expr past it); from?: [x,y,z] (From Object: start at that face or point); operation?: new|join|cut|intersect|auto (cut into a body, join out of one, else new); targets?: [body]; name?; body_name?"),
+        .params("distance: expr (or through_all: true); taper?: angle expr; sketch?: id|name (default: active or last sketch); profiles?: all | [index] | [[curve ids]] | [{point:[x,y]}]; face?: [x,y,z] (extrude a planar body face instead of a sketch profile); direction?: positive|negative|symmetric; distance2?; start_offset?: expr (along the sketch normal, also with direction: negative); to?: [x,y,z] (To Object: a face, vertex or point; to_offset?: expr past it); from?: [x,y,z] (From Object: start at that face or point); operation?: new|join|cut|intersect|auto (cut into a body, join out of one, else new); targets?: [body]; name?; body_name?"),
     CommandSpec::new("solid.revolve", "Revolve", revolve)
         .at("SOLID", "CREATE")
         .icon("revolve")

crates/engine/src/lib.rs

@@ -503,6 +503,8 @@ impl Session {
     }
 }
 
+#[cfg(test)]
+mod boolean_regression_tests;
 #[cfg(test)]
 mod construct_tests;
 #[cfg(test)]

crates/kernel/src/build.rs

@@ -795,13 +795,68 @@ pub enum PathSeg {
     },
 }
 
+impl PathSeg {
+    /// Unit directions of the path where the segment starts and ends.
+    fn tangents(&self) -> Option<(Vec3, Vec3)> {
+        match *self {
+            PathSeg::Line { a, b } => {
+                let t = (b - a).normalized()?;
+                Some((t, t))
+            }
+            PathSeg::Arc { a, center, axis, angle } => {
+                let k = axis.normalized()?;
+                let end = rotate_about(a, center, k, angle);
+                Some((k.cross(a - center).normalized()?, k.cross(end - center).normalized()?))
+            }
+        }
+    }
+
+    /// Points along the segment (an arc within `tol` of its chords), without its start.
+    fn points_after_start(&self, tol: f64) -> Vec<Vec3> {
+        match *self {
+            PathSeg::Line { b, .. } => vec![b],
+            PathSeg::Arc { a, center, axis, angle } => {
+                let Some(k) = axis.normalized() else { return vec![a] };
+                let r = (a - center).len();
+                let step = if r > tol { 2.0 * (1.0 - tol / r).acos() } else { angle.abs() };
+                let n = ((angle.abs() / step.max(1e-6)).ceil() as usize).clamp(1, 512);
+                (1..=n).map(|i| rotate_about(a, center, k, angle * i as f64 / n as f64)).collect()
+            }
+        }
+    }
+}
+
+/// `p` turned by `angle` about the unit axis `k` through `center` (right hand).
+fn rotate_about(p: Vec3, center: Vec3, k: Vec3, angle: f64) -> Vec3 {
+    let v = p - center;
+    let (s, c) = angle.sin_cos();
+    center + v * c + k.cross(v) * s + k * (k.dot(v) * (1.0 - c))
+}
+
 /// Sweep a planar region (outer loop) along a chain of path segments, keeping the profile
-/// perpendicular to the path (lines translate it, arcs rotate it).
+/// perpendicular to the path (lines translate it, arcs rotate it). A path with sharp corners
+/// is swept as a polyline, mitred at the corners.
 pub fn sweep(plane: &Plane, region: &Region2, path: &[PathSeg]) -> Result<Body> {
     region_ok(region)?;
     if path.is_empty() || path.len() > 1000 {
         return Err(KernelError::Invalid("the path needs 1…1000 segments".into()));
     }
+    // Segment by segment the profile only follows a path whose direction runs on across each
+    // joint: at a corner the next line would slide the profile within its own plane.
+    let ends: Vec<(Vec3, Vec3)> =
+        path.iter().map(PathSeg::tangents).collect::<Option<_>>().ok_or_else(|| KernelError::Invalid("a path segment has no direction".into()))?;
+    let corner = ends.windows(2).any(|w| w[0].1.dot(w[1].0) < (1e-3f64).cos());
+    if corner {
+        let first = match path.first() {
+            Some(PathSeg::Line { a, .. } | PathSeg::Arc { a, .. }) => *a,
+            None => return Err(KernelError::Invalid("the path needs 1…1000 segments".into())),
+        };
+        let mut pts = vec![first];
+        for seg in path {
+            pts.extend(seg.points_after_start(1e-2));
+        }
+        return crate::sweep_path::sweep_path(plane, region, &pts);
+    }
     // The outer loop and each hole, swept along the path; the ends are capped with holes.
     let mut starts = vec![wire(plane, &region.outer.ccw())?];
     for h in &region.holes {

crates/kernel/src/heal.rs

@@ -404,6 +404,12 @@ pub fn heal(solid: Solid, size: f64) -> Solid {
 
 /// [`heal`], keeping faces apart along lines they were split along on purpose.
 pub(crate) fn heal_keep(solid: Solid, size: f64, keep: &[Vec<Vec3>]) -> Solid {
+    // Edges a boolean made are intersection curves, which panic where they no longer
+    // evaluate; the solid then stays as it is.
+    crate::guard("heal", || Ok(heal_shells(&solid, size, keep))).ok().flatten().unwrap_or(solid)
+}
+
+fn heal_shells(solid: &Solid, size: f64, keep: &[Vec<Vec3>]) -> Option<Solid> {
     let tol = (size * 1e-7).max(1e-9);
     let mut shells = Vec::new();
     let mut changed = false;
@@ -428,9 +434,9 @@ pub(crate) fn heal_keep(solid: Solid, size: f64, keep: &[Vec<Vec3>]) -> Solid {
         shells.push(cur);
     }
     if !changed {
-        return solid;
+        return None;
     }
-    Solid::try_new(shells).unwrap_or(solid)
+    Solid::try_new(shells).ok()
 }
 
 #[cfg(test)]

crates/kernel/src/ops.rs

@@ -189,11 +189,28 @@ pub fn boolean(a: &Body, b: &Body, op: BoolOp) -> Result<Option<Body>> {
     r
 }
 
+/// Smallest side of the box where the bodies' vertex boxes overlap (infinite when they don't).
+fn meeting_size(a: &Body, b: &Body) -> f64 {
+    let bbox = |x: &Body| {
+        let mut bb = solvecraft_geom::Aabb3::EMPTY;
+        for v in x.solid.vertex_iter() {
+            bb.add(from_p3(v.point()));
+        }
+        bb
+    };
+    let (ba, bb) = (bbox(a), bbox(b));
+    let s = solvecraft_geom::Aabb3 { min: ba.min.max(bb.min), max: ba.max.min(bb.max) }.size();
+    if s.x > 0.0 && s.y > 0.0 && s.z > 0.0 { s.x.min(s.y).min(s.z) } else { f64::INFINITY }
+}
+
 fn boolean_whole(a: &Body, b: &Body, op: BoolOp) -> Result<Option<Body>> {
     a.require_brep("a boolean")?;
     b.require_brep("a boolean")?;
     let size = a.size().max(b.size());
-    let tol_m = size * 5e-4;
+    // The checks mesh both bodies; where they meet in a region much thinner than they are (a
+    // thin part trimmed by a long one), a step from their size would facet the small faces
+    // there coarsely enough to fail the membership test on a correct result.
+    let tol_m = (size * 5e-4).min(meeting_size(a, b) * 2e-3).max(size * 1e-5);
     let meshes = (a.tessellate(tol_m).ok(), b.tessellate(tol_m).ok());
     let (va, vb) = match &meshes {
         (Some(ma), Some(mb)) => (ma.measure().volume, mb.measure().volume),
@@ -296,7 +313,7 @@ fn boolean_whole(a: &Body, b: &Body, op: BoolOp) -> Result<Option<Body>> {
                     empty_votes += 1;
                     last = "empty result".into();
                 }
-                Ok(s) => match Body::new(crate::heal::heal_keep(s, size, &keep)) {
+                Ok(s) => match guard("boolean result", || Body::new(crate::heal::heal_keep(s, size, &keep))) {
                     Ok(body) => {
                         let (ok, v) = check(&body);
                         if ok {
@@ -498,6 +515,15 @@ pub fn transform(body: &Body, translate: Vec3, origin: Vec3, axis: Vec3, angle:
 /// Split a body with a plane: the parts on the positive and the negative side of the plane
 /// (only the non-empty ones).
 pub fn split_by_plane(body: &Body, plane: &solvecraft_geom::Plane) -> Result<Vec<Body>> {
+    let mut out = Vec::new();
+    for positive in [true, false] {
+        out.extend(half_space(body, plane, positive)?);
+    }
+    Ok(out)
+}
+
+/// The part of a body on one side of a plane (`None` when nothing is there).
+pub(crate) fn half_space(body: &Body, plane: &solvecraft_geom::Plane, positive: bool) -> Result<Option<Body>> {
     body.require_brep("split")?;
     let size = body.size();
     let mut bb = solvecraft_geom::Aabb3::EMPTY;
@@ -515,12 +541,7 @@ pub fn split_by_plane(body: &Body, plane: &solvecraft_geom::Plane) -> Result<Vec
         ]),
         holes: vec![],
     };
-    let mut out = Vec::new();
-    for (lo, hi) in [(0.0, reach), (-reach, 0.0)] {
-        let half = crate::build::extrude(plane, std::slice::from_ref(&sq), lo, hi)?.pop().ok_or_else(|| KernelError::Failed("half space".into()))?;
-        if let Some(part) = boolean(body, &half, BoolOp::Intersect)? {
-            out.push(part);
-        }
-    }
-    Ok(out)
+    let (lo, hi) = if positive { (0.0, reach) } else { (-reach, 0.0) };
+    let half = crate::build::extrude(plane, std::slice::from_ref(&sq), lo, hi)?.pop().ok_or_else(|| KernelError::Failed("half space".into()))?;
+    boolean(body, &half, BoolOp::Intersect)
 }

crates/kernel/src/polybool.rs

@@ -306,12 +306,118 @@ pub fn planar_boolean(a: &Body, b: &Body, op: BoolOp) -> Result<Option<Body>> {
     a.require_brep("a boolean")?;
     b.require_brep("a boolean")?;
     let (Some(pa), Some(pb)) = (planar_polys(a), planar_polys(b)) else {
+        if let Some(r) = sliced_planar_boolean(a, b, op) {
+            return r;
+        }
         return Err(KernelError::Failed("not supported yet: coincident faces on curved bodies".into()));
     };
     let size = a.size().max(b.size());
     polys_boolean(pa, pb, size, op, size * 1e-9)
 }
 
+/// [`planar_boolean`] of a body with curved faces and an all-planar one, when a slab between
+/// two parallel planes holds `b` and none of `a`'s curved faces (a pocket cut beside a round
+/// hole): `a` is split across the slab, the middle piece takes the exact boolean, and the
+/// pieces outside are joined back on. `None` when no such slab exists.
+fn sliced_planar_boolean(a: &Body, b: &Body, op: BoolOp) -> Option<Result<Option<Body>>> {
+    let pb = planar_polys(b)?;
+    let size = a.size().max(b.size());
+    let tol = (a.size() * 1e-3).max(1e-3);
+    let infos = a.faces(tol).ok()?;
+    let mesh = a.tessellate(tol).ok()?;
+    let mut curved: HashMap<u32, solvecraft_geom::Aabb3> = HashMap::new();
+    for (t, fi) in mesh.triangles.iter().zip(&mesh.tri_face) {
+        if infos.get(*fi as usize)?.plane_normal.is_some() {
+            continue;
+        }
+        let bx = curved.entry(*fi).or_insert(solvecraft_geom::Aabb3::EMPTY);
+        for p in mesh.tri(t)? {
+            bx.add(p);
+        }
+    }
+    let mut tool = solvecraft_geom::Aabb3::EMPTY;
+    for p in pb.iter().flat_map(|p| &p.v) {
+        tool.add(*p);
+    }
+    let at = |v: Vec3, k: usize| [v.x, v.y, v.z].get(k).copied().unwrap_or(0.0);
+    let axes = [Vec3::X, Vec3::Y, Vec3::Z];
+    // Planar faces of `a` square to an axis: a cut must not land in one of their planes.
+    let flat: Vec<(Vec3, Vec3)> = infos.iter().filter_map(|f| Some((f.plane_normal?, f.centroid))).collect();
+    let gap = size * 1e-4;
+    let place = |k: usize, near: f64, far: f64| -> Option<f64> {
+        let axis = *axes.get(k)?;
+        [0.5, 0.37, 0.63, 0.25, 0.75]
+            .iter()
+            .map(|t| near + (far - near) * t)
+            .find(|x| !flat.iter().any(|(n, c)| n.dot(axis).abs() > 1.0 - 1e-9 && (at(*c, k) - x).abs() < gap))
+    };
+    // An axis along which every curved face lies clear of the tool: the cuts either side.
+    let (k, lo, hi) = (0..3).find_map(|k| {
+        let (t0, t1) = (at(tool.min, k), at(tool.max, k));
+        let (mut below, mut above) = (None::<f64>, None::<f64>);
+        for bx in curved.values() {
+            let (c0, c1) = (at(bx.min, k), at(bx.max, k));
+            if c1 < t0 - gap {
+                below = Some(below.map_or(c1, |x| x.max(c1)));
+            } else if c0 > t1 + gap {
+                above = Some(above.map_or(c0, |x| x.min(c0)));
+            } else {
+                return None;
+            }
+        }
+        let lo = match below {
+            Some(c) => Some(place(k, t0, c)?),
+            None => None,
+        };
+        let hi = match above {
+            Some(c) => Some(place(k, t1, c)?),
+            None => None,
+        };
+        Some((k, lo, hi))
+    })?;
+    let (u, w) = match k {
+        0 => (Vec3::Y, Vec3::Z),
+        1 => (Vec3::Z, Vec3::X),
+        _ => (Vec3::X, Vec3::Y),
+    };
+    let axis = *axes.get(k)?;
+    let run = || -> Result<Option<Body>> {
+        let fail = |m: &str| KernelError::Failed(format!("planar boolean beside curved faces: {m}"));
+        let mut middle = a.clone();
+        let mut outside = Vec::new();
+        for (cut, keep_positive) in [(lo, true), (hi, false)] {
+            let Some(x) = cut else { continue };
+            let plane = solvecraft_geom::Plane::new(axis * x, u, w).ok_or_else(|| fail("plane"))?;
+            outside.extend(crate::ops::half_space(&middle, &plane, !keep_positive)?);
+            middle = crate::ops::half_space(&middle, &plane, keep_positive)?.ok_or_else(|| fail("nothing beside the tool"))?;
+        }
+        let pm = planar_polys(&middle).ok_or_else(|| fail("the middle piece is not planar"))?;
+        let r = polys_boolean(pm, pb, size, op, size * 1e-9)?;
+        if op == BoolOp::Intersect {
+            return Ok(r);
+        }
+        // The rebuild leaves the cut faces' edges split where other polygons were cut; healed,
+        // they match the pieces outside edge for edge.
+        let mut pieces: Vec<Body> = match r {
+            Some(m) => vec![Body::new(crate::heal::heal_keep(m.deep_copy(), size, a.split_keep()))?],
+            None => Vec::new(),
+        };
+        // Each piece outside meets the middle along its cut face (glued one shell to another).
+        if pieces.iter().chain(&outside).any(|x| x.solid.boundaries().len() != 1) {
+            return Err(fail("a piece in several parts"));
+        }
+        for o in outside {
+            let joined = pieces.iter().enumerate().find_map(|(i, m)| Some((i, crate::coplanar::glue(m, &o)?.ok()?)));
+            let (i, j) = joined.ok_or_else(|| fail("a piece does not join back"))?;
+            if let Some(slot) = pieces.get_mut(i) {
+                *slot = j;
+            }
+        }
+        crate::ops::join_pieces(pieces)
+    };
+    Some(run())
+}
+
 fn polys_boolean(pa: Vec<Poly>, pb: Vec<Poly>, size: f64, op: BoolOp, eps: f64) -> Result<Option<Body>> {
     let polys = csg(pa, pb, op, eps)?;
     if polys.is_empty() {

crates/kernel/src/tests.rs

@@ -258,6 +258,38 @@ fn loft_and_sweep() {
     assert!(rel(measure(&t).unwrap().volume, expect) < 1e-3, "{} vs {expect}", measure(&t).unwrap().volume);
 }
 
+#[test]
+fn sweep_round_a_sharp_corner_is_mitred() {
+    // A 2 × 2 square centred on an L of two 10 mm lines: lines alone would slide the profile
+    // along the second leg within its own plane (an empty, inside-out body). Mitred, the
+    // volume is the area times the path length.
+    let sq =
+        Region2 { outer: Loop2::polygon(&[Vec2::new(-1.0, -1.0), Vec2::new(1.0, -1.0), Vec2::new(1.0, 1.0), Vec2::new(-1.0, 1.0)]), holes: vec![] };
+    let path = [
+        PathSeg::Line { a: Vec3::ZERO, b: Vec3::new(0.0, 0.0, 10.0) },
+        PathSeg::Line { a: Vec3::new(0.0, 0.0, 10.0), b: Vec3::new(10.0, 0.0, 10.0) },
+    ];
+    let b = sweep(&Plane::XY, &sq, &path).unwrap();
+    let m = measure(&b).unwrap();
+    assert!(rel(m.volume, 80.0) < 1e-9, "{}", m.volume);
+}
+
+#[test]
+fn moving_a_body_with_an_intersection_curve_edge_ends() {
+    // A cylinder cut by a tilted box keeps an elliptical intersection curve as an edge. Moving
+    // such a body scales the curve's tolerance by the matrix norm, whose root refinement cycled
+    // forever for about one translation in eight (this one included): the app hung.
+    let cyl = cylinder(Vec3::ZERO, Vec3::Z, 5.0, 40.0).unwrap();
+    let block = box_solid(Vec3::new(-20.0, -20.0, 20.0), Vec3::new(20.0, 20.0, 60.0)).unwrap();
+    let block = crate::ops::transform(&block, Vec3::ZERO, Vec3::new(0.0, 0.0, 20.0), Vec3::X, 0.4).unwrap();
+    let cut = boolean(&cyl, &block, BoolOp::Cut).unwrap().unwrap();
+    assert!(cut.solid.edge_iter().any(|e| matches!(e.curve(), truck_modeling::Curve::IntersectionCurve(_))));
+    let t = Vec3::new(0.6911, -0.3337, -0.6127) * 2.334;
+    let moved = crate::ops::transform(&cut, t, Vec3::ZERO, Vec3::Z, 0.0).unwrap();
+    let (v0, v1) = (measure(&cut).unwrap().volume, measure(&moved).unwrap().volume);
+    assert!(rel(v0, v1) < 1e-6, "{v0} vs {v1}");
+}
+
 #[test]
 fn planar_booleans_with_coincident_faces() {
     // L-bracket: a wall flush with three sides of its base.

crates/ui-egui/src/documents.rs

@@ -298,13 +298,13 @@ pub fn tabs(app: &mut SolveApp, ui: &mut egui::Ui, r: Rect, x0: f32, x1: f32) {
     let response = ui.interact(home, ui.id().with("home_tab"), Sense::click());
     ui.painter().rect_filled(
         home,
-        5.0,
+        egui::CornerRadius { nw: 5, ne: 5, sw: 0, se: 0 },
         if app.home.open {
             t.toolbar
         } else if response.hovered() {
             t.hover
         } else {
-            Color32::from_white_alpha(10)
+            t.app_bar
         },
     );
     ui.painter().text(home.center(), Align2::CENTER_CENTER, "Home", FontId::proportional(13.0), if app.home.open { t.text } else { t.app_bar_text });
@@ -338,7 +338,7 @@ pub fn tabs(app: &mut SolveApp, ui: &mut egui::Ui, r: Rect, x0: f32, x1: f32) {
         } else if resp.hovered() {
             Color32::from_white_alpha(28)
         } else {
-            Color32::from_white_alpha(10)
+            t.app_bar
         };
         ui.painter().rect_filled(tr, egui::CornerRadius { nw: 5, ne: 5, sw: 0, se: 0 }, fill);
         let ink = if active { t.text } else { t.app_bar_text };

crates/ui-egui/src/theme.rs

@@ -16,7 +16,6 @@ pub fn is_dark() -> bool {
 pub struct Tokens {
     pub app_bar: Color32,
     pub toolbar: Color32,
-    pub tab_active: Color32,
     pub panel: Color32,
     pub panel_header: Color32,
     pub border: Color32,
@@ -132,7 +131,6 @@ impl Tokens {
         Tokens {
             app_bar: Color32::from_rgb(52, 56, 64),
             toolbar: Color32::from_rgb(246, 247, 249),
-            tab_active: Color32::from_rgb(255, 255, 255),
             panel: Color32::from_rgb(250, 251, 252),
             panel_header: Color32::from_rgb(236, 238, 242),
             border: Color32::from_rgb(208, 212, 220),
@@ -212,7 +210,6 @@ impl Tokens {
         let mut t = Tokens::light();
         t.app_bar = Color32::from_rgb(28, 31, 38);
         t.toolbar = Color32::from_rgb(37, 41, 49);
-        t.tab_active = Color32::from_rgb(46, 52, 64);
         t.panel = Color32::from_rgb(39, 43, 52);
         t.panel_header = Color32::from_rgb(47, 52, 62);
         t.border = Color32::from_rgb(63, 69, 82);

crates/ui-egui/src/toolbar.rs

@@ -38,7 +38,7 @@ pub fn promoted(tab: &str, panel: &str) -> usize {
 
 pub fn app_bar(app: &mut SolveApp, ui: &mut egui::Ui) {
     let t = Tokens::get();
-    egui::Panel::top("sc_appbar").exact_size(34.0).frame(egui::Frame::NONE.fill(t.app_bar)).show(ui, |ui| {
+    egui::Panel::top("sc_appbar").exact_size(34.0).show_separator_line(false).frame(egui::Frame::NONE.fill(t.app_bar)).show(ui, |ui| {
         let r = ui.max_rect();
         // The bar is the window's title bar: its empty space moves the window.
         if app.custom_titlebar || app.integrated_titlebar {
@@ -305,7 +305,6 @@ pub fn toolbar(app: &mut SolveApp, ui: &mut egui::Ui) {
             let active = app.ui.tab == tab;
             let resp = ui.interact(tr, ui.id().with(("tab", tab)), Sense::click());
             if active {
-                painter.rect_filled(tr, egui::CornerRadius { nw: 4, ne: 4, sw: 0, se: 0 }, t.tab_active);
                 let c = if tab == "SKETCH" { t.sketch_accent } else { t.accent };
                 painter.line_segment([pos2(tr.left() + 4.0, tr.bottom() - 1.0), pos2(tr.right() - 4.0, tr.bottom() - 1.0)], Stroke::new(2.0, c));
             } else if resp.hovered() {

examples/regressions/loft_cut_beside_round_hole.json

@@ -0,0 +1,24 @@
+{
+ "description": "A loft cut from a pocket floor beside a round hole it never touches (the hole made first). Slab 900 x 480 x 160, a 557.68 x 100 pocket 120 deep, a 35 mm hole 25 deep 66 mm from the pocket, then a loft cut from a rectangle on the pocket floor (0.5 mm inside its walls) down 4 mm to an 80 x 80 square. The loft's top lies in the pocket floor, and the hole's cylinder kept the exact planar boolean from taking it: the cut failed with a kernel panic. Result: the same body as with the hole made after the loft.",
+ "commands": [
+  {"command": "sketch.create", "params": {"plane": "XY", "offset": 690, "name": "Outline"}},
+  {"command": "sketch.rectangle.two_point", "params": {"p0": [0, -480], "p1": [900, 0]}},
+  {"command": "sketch.finish", "params": {}},
+  {"command": "solid.extrude", "params": {"distance": 160, "name": "Slab"}},
+  {"command": "sketch.create", "params": {"plane": "XY", "offset": 850, "name": "BowlTop"}},
+  {"command": "sketch.rectangle.two_point", "params": {"p0": [171.16, -151.16], "p1": [728.84, -51.16]}},
+  {"command": "sketch.finish", "params": {}},
+  {"command": "solid.extrude", "params": {"distance": -120, "operation": "cut", "name": "BowlPocket"}},
+  {"command": "sketch.create", "params": {"plane": "XY", "offset": 850, "name": "TapHole"}},
+  {"command": "sketch.circle.center", "params": {"center": [795, -101.16], "radius": 17.5}},
+  {"command": "sketch.finish", "params": {}},
+  {"command": "solid.extrude", "params": {"distance": -25, "operation": "cut", "name": "TapCut"}},
+  {"command": "sketch.create", "params": {"plane": "XY", "offset": 730, "name": "TroughInset"}},
+  {"command": "sketch.rectangle.two_point", "params": {"p0": [171.66, -150.66], "p1": [728.34, -51.66]}},
+  {"command": "sketch.finish", "params": {}},
+  {"command": "sketch.create", "params": {"plane": "XY", "offset": 726, "name": "DrainSeat"}},
+  {"command": "sketch.rectangle.two_point", "params": {"p0": [410, -141.16], "p1": [490, -61.16]}},
+  {"command": "sketch.finish", "params": {}},
+  {"command": "solid.loft", "params": {"sections": [{"sketch": "TroughInset"}, {"sketch": "DrainSeat"}], "operation": "cut", "name": "TroughFall"}}
+ ]
+}

examples/regressions/oblique_cut_of_cylinder.json

@@ -0,0 +1,15 @@
+{
+ "description": "A thin 437 mm bar whose profile has a 3 mm concave arc, cut by a long prism whose top is tilted 21 degrees to the bar (an elliptical cut through the arc's cylinder). The check on the result meshed both bodies at a step from the 920 mm prism, which faceted the arc so coarsely that a correct cut failed the membership test. Result: the bar with its end trimmed by the tilted plane, as a split by that plane gives.",
+ "commands": [
+  {"command": "sketch.create", "params": {"plane": "XZ", "offset": -462, "name": "Profile"}},
+  {"command": "sketch.line", "params": {"points": [[147.4407, -849.95], [147.4407, -851.0], [151.3952, -851.0], [151.3952, -847.4709], [150.3952, -847.4709]]}},
+  {"command": "sketch.arc.center_point", "params": {"center": [147.4407, -846.95], "start": "l1.start", "end": "l4.end"}},
+  {"command": "sketch.finish", "params": {}},
+  {"command": "solid.extrude", "params": {"sketch": "Profile", "profiles": "all", "direction": "positive", "start_offset": -1, "distance": 437.52, "operation": "new", "body_name": "Tool"}},
+  {"command": "sketch.create", "params": {"plane": "YZ", "name": "Clip"}},
+  {"command": "sketch.line", "params": {"points": [[-472, 853.91], [-432, 838.47], [-432, 820], [-472, 820]], "closed": true}},
+  {"command": "sketch.finish", "params": {}},
+  {"command": "solid.extrude", "params": {"sketch": "Clip", "profiles": "all", "start_offset": -10, "distance": 920, "operation": "new", "body_name": "ClipBody"}},
+  {"command": "solid.combine", "params": {"target": "Tool", "tools": ["ClipBody"], "operation": "cut", "keep_tools": true}}
+ ]
+}

examples/regressions/oblique_cut_of_cylinder_hang.json

@@ -0,0 +1,15 @@
+{
+ "description": "The oblique cut of oblique_cut_of_cylinder.json with the prism running past the bar's end. Moving a boolean result back after a shifted retry computed a matrix norm whose root refinement never ended: the app hung at full CPU. Result: the same trimmed bar.",
+ "commands": [
+  {"command": "sketch.create", "params": {"plane": "XZ", "offset": -462, "name": "Profile"}},
+  {"command": "sketch.line", "params": {"points": [[147.4407, -849.95], [147.4407, -851.0], [151.3952, -851.0], [151.3952, -847.4709], [150.3952, -847.4709]]}},
+  {"command": "sketch.arc.center_point", "params": {"center": [147.4407, -846.95], "start": "l1.start", "end": "l4.end"}},
+  {"command": "sketch.finish", "params": {}},
+  {"command": "solid.extrude", "params": {"sketch": "Profile", "profiles": "all", "direction": "positive", "start_offset": -1, "distance": 437.52, "operation": "new", "body_name": "Tool"}},
+  {"command": "sketch.create", "params": {"plane": "YZ", "name": "Clip"}},
+  {"command": "sketch.line", "params": {"points": [[-480, 856.9979999999999], [-432, 838.47], [-432, 820], [-480, 820]], "closed": true}},
+  {"command": "sketch.finish", "params": {}},
+  {"command": "solid.extrude", "params": {"sketch": "Clip", "profiles": "all", "start_offset": -10, "distance": 920, "operation": "new", "body_name": "ClipBody"}},
+  {"command": "solid.combine", "params": {"target": "Tool", "tools": ["ClipBody"], "operation": "cut", "keep_tools": true}}
+ ]
+}

examples/regressions/shallow_loft_cut_from_top_face.json

@@ -0,0 +1,16 @@
+{
+ "description": "A shallow, wide loft cut starting in a body's top face: slab 900 x 480 x 160, loft cut from a 556.68 x 99 rectangle on the top face down 4 mm to an 80 x 80 square (sides about 0.9 degrees from horizontal). Healing an intermediate boolean result panicked outside the boolean's guard and stopped every later attempt. Result: the slab less the loft, 69 003 015.31 mm³.",
+ "commands": [
+  {"command": "sketch.create", "params": {"plane": "XY", "name": "Base"}},
+  {"command": "sketch.rectangle.two_point", "params": {"p0": [0, 0], "p1": [900, 480]}},
+  {"command": "sketch.finish", "params": {}},
+  {"command": "solid.extrude", "params": {"distance": 160, "name": "Slab"}},
+  {"command": "sketch.create", "params": {"plane": "XY", "offset": 160, "name": "LoftTop"}},
+  {"command": "sketch.rectangle.two_point", "params": {"p0": [171.66, 329.34], "p1": [728.34, 428.34]}},
+  {"command": "sketch.finish", "params": {}},
+  {"command": "sketch.create", "params": {"plane": "XY", "offset": 156, "name": "LoftBottom"}},
+  {"command": "sketch.rectangle.two_point", "params": {"p0": [410, 338.84], "p1": [490, 418.84]}},
+  {"command": "sketch.finish", "params": {}},
+  {"command": "solid.loft", "params": {"sections": [{"sketch": "LoftTop"}, {"sketch": "LoftBottom"}], "operation": "cut", "name": "LoftCut"}}
+ ]
+}

vendor/README.md

@@ -49,3 +49,15 @@ Changes:
   direction (long triangles from a trim loop to the grid's far ends cut through cylinders).
 - `lib.rs`: compiler warnings no longer fail the build. Example, test and bench targets are
   dropped (their sources are not vendored).
+
+## matext4cgmath 0.1.0 (Apache-2.0, © RICOS Co. Ltd., https://github.com/ricosjp/matext4cgmath)
+
+truck-geometry uses it for a matrix's operator norm (an intersection curve scales its tolerance
+by it whenever it is moved).
+
+Changes:
+- `solver.rs`: the Newton refinement of cubic and quartic roots stops after 100 steps. Near a
+  double root it could cycle forever: moving a body whose edges a boolean made (the retries
+  move both bodies and the result back) hung the app.
+- `lib.rs`: compiler warnings no longer fail the build. Test targets are dropped (their sources
+  are not vendored).

vendor/matext4cgmath/Cargo.toml

@@ -0,0 +1,40 @@
+# THIS FILE IS AUTOMATICALLY GENERATED BY CARGO
+#
+# When uploading crates to the registry Cargo will automatically
+# "normalize" Cargo.toml files for maximal compatibility
+# with all versions of Cargo and also rewrite `path` dependencies
+# to registry (e.g., crates.io) dependencies.
+#
+# If you are reading this file be aware that the original Cargo.toml
+# will likely look very different (and much more reasonable).
+# See Cargo.toml.orig for the original contents.
+
+[package]
+edition = "2021"
+name = "matext4cgmath"
+version = "0.1.0"
+authors = ["Yoshinori Tanimura <tanimura@ricos.co.jp>"]
+description = "Unofficial third-party cgmath extensions for calculate eigenvalues, operator norms and Iwasawa decomposition."
+homepage = "https://github.com/ricosjp/matext4cgmath"
+keywords = [
+    "cgmath",
+    "graphics",
+    "linear-algebra",
+]
+categories = ["graphics"]
+license = "Apache-2.0"
+repository = "https://github.com/ricosjp/matext4cgmath"
+resolver = "2"
+autotests = false
+autoexamples = false
+autobenches = false
+
+[dependencies.cgmath]
+version = "0.18.0"
+
+[dependencies.katexit]
+version = "0.1.2"
+
+[dependencies.num-complex]
+version = "0.4.0"
+

vendor/matext4cgmath/LICENSE

@@ -0,0 +1,201 @@
+                                 Apache License
+                           Version 2.0, January 2004
+                        http://www.apache.org/licenses/
+
+   TERMS AND CONDITIONS FOR USE, REPRODUCTION, AND DISTRIBUTION
+
+   1. Definitions.
+
+      "License" shall mean the terms and conditions for use, reproduction,
+      and distribution as defined by Sections 1 through 9 of this document.
+
+      "Licensor" shall mean the copyright owner or entity authorized by
+      the copyright owner that is granting the License.
+
+      "Legal Entity" shall mean the union of the acting entity and all
+      other entities that control, are controlled by, or are under common
+      control with that entity. For the purposes of this definition,
+      "control" means (i) the power, direct or indirect, to cause the
+      direction or management of such entity, whether by contract or
+      otherwise, or (ii) ownership of fifty percent (50%) or more of the
+      outstanding shares, or (iii) beneficial ownership of such entity.
+
+      "You" (or "Your") shall mean an individual or Legal Entity
+      exercising permissions granted by this License.
+
+      "Source" form shall mean the preferred form for making modifications,
+      including but not limited to software source code, documentation
+      source, and configuration files.
+
+      "Object" form shall mean any form resulting from mechanical
+      transformation or translation of a Source form, including but
+      not limited to compiled object code, generated documentation,
+      and conversions to other media types.
+
+      "Work" shall mean the work of authorship, whether in Source or
+      Object form, made available under the License, as indicated by a
+      copyright notice that is included in or attached to the work
+      (an example is provided in the Appendix below).
+
+      "Derivative Works" shall mean any work, whether in Source or Object
+      form, that is based on (or derived from) the Work and for which the
+      editorial revisions, annotations, elaborations, or other modifications
+      represent, as a whole, an original work of authorship. For the purposes
+      of this License, Derivative Works shall not include works that remain
+      separable from, or merely link (or bind by name) to the interfaces of,
+      the Work and Derivative Works thereof.
+
+      "Contribution" shall mean any work of authorship, including
+      the original version of the Work and any modifications or additions
+      to that Work or Derivative Works thereof, that is intentionally
+      submitted to Licensor for inclusion in the Work by the copyright owner
+      or by an individual or Legal Entity authorized to submit on behalf of
+      the copyright owner. For the purposes of this definition, "submitted"
+      means any form of electronic, verbal, or written communication sent
+      to the Licensor or its representatives, including but not limited to
+      communication on electronic mailing lists, source code control systems,
+      and issue tracking systems that are managed by, or on behalf of, the
+      Licensor for the purpose of discussing and improving the Work, but
+      excluding communication that is conspicuously marked or otherwise
+      designated in writing by the copyright owner as "Not a Contribution."
+
+      "Contributor" shall mean Licensor and any individual or Legal Entity
+      on behalf of whom a Contribution has been received by Licensor and
+      subsequently incorporated within the Work.
+
+   2. Grant of Copyright License. Subject to the terms and conditions of
+      this License, each Contributor hereby grants to You a perpetual,
+      worldwide, non-exclusive, no-charge, royalty-free, irrevocable
+      copyright license to reproduce, prepare Derivative Works of,
+      publicly display, publicly perform, sublicense, and distribute the
+      Work and such Derivative Works in Source or Object form.
+
+   3. Grant of Patent License. Subject to the terms and conditions of
+      this License, each Contributor hereby grants to You a perpetual,
+      worldwide, non-exclusive, no-charge, royalty-free, irrevocable
+      (except as stated in this section) patent license to make, have made,
+      use, offer to sell, sell, import, and otherwise transfer the Work,
+      where such license applies only to those patent claims licensable
+      by such Contributor that are necessarily infringed by their
+      Contribution(s) alone or by combination of their Contribution(s)
+      with the Work to which such Contribution(s) was submitted. If You
+      institute patent litigation against any entity (including a
+      cross-claim or counterclaim in a lawsuit) alleging that the Work
+      or a Contribution incorporated within the Work constitutes direct
+      or contributory patent infringement, then any patent licenses
+      granted to You under this License for that Work shall terminate
+      as of the date such litigation is filed.
+
+   4. Redistribution. You may reproduce and distribute copies of the
+      Work or Derivative Works thereof in any medium, with or without
+      modifications, and in Source or Object form, provided that You
+      meet the following conditions:
+
+      (a) You must give any other recipients of the Work or
+          Derivative Works a copy of this License; and
+
+      (b) You must cause any modified files to carry prominent notices
+          stating that You changed the files; and
+
+      (c) You must retain, in the Source form of any Derivative Works
+          that You distribute, all copyright, patent, trademark, and
+          attribution notices from the Source form of the Work,
+          excluding those notices that do not pertain to any part of
+          the Derivative Works; and
+
+      (d) If the Work includes a "NOTICE" text file as part of its
+          distribution, then any Derivative Works that You distribute must
+          include a readable copy of the attribution notices contained
+          within such NOTICE file, excluding those notices that do not
+          pertain to any part of the Derivative Works, in at least one
+          of the following places: within a NOTICE text file distributed
+          as part of the Derivative Works; within the Source form or
+          documentation, if provided along with the Derivative Works; or,
+          within a display generated by the Derivative Works, if and
+          wherever such third-party notices normally appear. The contents
+          of the NOTICE file are for informational purposes only and
+          do not modify the License. You may add Your own attribution
+          notices within Derivative Works that You distribute, alongside
+          or as an addendum to the NOTICE text from the Work, provided
+          that such additional attribution notices cannot be construed
+          as modifying the License.
+
+      You may add Your own copyright statement to Your modifications and
+      may provide additional or different license terms and conditions
+      for use, reproduction, or distribution of Your modifications, or
+      for any such Derivative Works as a whole, provided Your use,
+      reproduction, and distribution of the Work otherwise complies with
+      the conditions stated in this License.
+
+   5. Submission of Contributions. Unless You explicitly state otherwise,
+      any Contribution intentionally submitted for inclusion in the Work
+      by You to the Licensor shall be under the terms and conditions of
+      this License, without any additional terms or conditions.
+      Notwithstanding the above, nothing herein shall supersede or modify
+      the terms of any separate license agreement you may have executed
+      with Licensor regarding such Contributions.
+
+   6. Trademarks. This License does not grant permission to use the trade
+      names, trademarks, service marks, or product names of the Licensor,
+      except as required for reasonable and customary use in describing the
+      origin of the Work and reproducing the content of the NOTICE file.
+
+   7. Disclaimer of Warranty. Unless required by applicable law or
+      agreed to in writing, Licensor provides the Work (and each
+      Contributor provides its Contributions) on an "AS IS" BASIS,
+      WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or
+      implied, including, without limitation, any warranties or conditions
+      of TITLE, NON-INFRINGEMENT, MERCHANTABILITY, or FITNESS FOR A
+      PARTICULAR PURPOSE. You are solely responsible for determining the
+      appropriateness of using or redistributing the Work and assume any
+      risks associated with Your exercise of permissions under this License.
+
+   8. Limitation of Liability. In no event and under no legal theory,
+      whether in tort (including negligence), contract, or otherwise,
+      unless required by applicable law (such as deliberate and grossly
+      negligent acts) or agreed to in writing, shall any Contributor be
+      liable to You for damages, including any direct, indirect, special,
+      incidental, or consequential damages of any character arising as a
+      result of this License or out of the use or inability to use the
+      Work (including but not limited to damages for loss of goodwill,
+      work stoppage, computer failure or malfunction, or any and all
+      other commercial damages or losses), even if such Contributor
+      has been advised of the possibility of such damages.
+
+   9. Accepting Warranty or Additional Liability. While redistributing
+      the Work or Derivative Works thereof, You may choose to offer,
+      and charge a fee for, acceptance of support, warranty, indemnity,
+      or other liability obligations and/or rights consistent with this
+      License. However, in accepting such obligations, You may act only
+      on Your own behalf and on Your sole responsibility, not on behalf
+      of any other Contributor, and only if You agree to indemnify,
+      defend, and hold each Contributor harmless for any liability
+      incurred by, or claims asserted against, such Contributor by reason
+      of your accepting any such warranty or additional liability.
+
+   END OF TERMS AND CONDITIONS
+
+   APPENDIX: How to apply the Apache License to your work.
+
+      To apply the Apache License to your work, attach the following
+      boilerplate notice, with the fields enclosed by brackets "[]"
+      replaced with your own identifying information. (Don't include
+      the brackets!)  The text should be enclosed in the appropriate
+      comment syntax for the file format. We also recommend that a
+      file or class name and description of purpose be included on the
+      same "printed page" as the copyright notice for easier
+      identification within third-party archives.
+
+   Copyright 2022 Research Institute for Computational Science Co. Ltd
+
+   Licensed under the Apache License, Version 2.0 (the "License");
+   you may not use this file except in compliance with the License.
+   You may obtain a copy of the License at
+
+       http://www.apache.org/licenses/LICENSE-2.0
+
+   Unless required by applicable law or agreed to in writing, software
+   distributed under the License is distributed on an "AS IS" BASIS,
+   WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+   See the License for the specific language governing permissions and
+   limitations under the License.

vendor/matext4cgmath/README.md

@@ -0,0 +1,5 @@
+# matext4cgmath
+
+[![Crates.io](https://img.shields.io/crates/v/matext4cgmath.svg)](https://crates.io/crates/matext4cgmath) [![Docs.rs](https://docs.rs/matext4cgmath/badge.svg)](https://docs.rs/matext4cgmath)
+
+Unofficial third-party cgmath extensions for calculate eigenvalues, operator norms and Iwasawa decomposition.

vendor/matext4cgmath/src/eigens.rs

@@ -0,0 +1,153 @@
+use crate::*;
+
+impl<F: BaseFloat> EigenValues for Matrix2<F> {
+    type EigenValues = [Complex<F>; 2];
+    fn eigenvalues(self) -> [Complex<F>; 2] {
+        solver::solve_quadratic(-self.trace(), self.determinant())
+    }
+}
+
+impl<F: BaseFloat> EigenValues for Matrix3<F> {
+    type EigenValues = [Complex<F>; 3];
+    fn eigenvalues(self) -> [Complex<F>; 3] {
+        let poly2 = self[0][0] * self[1][1] + self[0][0] * self[2][2] + self[1][1] * self[2][2]
+            - self[0][1] * self[1][0]
+            - self[1][2] * self[2][1]
+            - self[0][2] * self[2][0];
+        solver::solve_cubic(-self.trace(), poly2, -self.determinant())
+    }
+}
+
+impl<F: BaseFloat> EigenValues for Matrix4<F> {
+    type EigenValues = [Complex<F>; 4];
+    fn eigenvalues(self) -> [Complex<F>; 4] {
+        let poly2 = self[0][0] * self[1][1]
+            + self[0][0] * self[2][2]
+            + self[0][0] * self[3][3]
+            + self[1][1] * self[2][2]
+            + self[1][1] * self[3][3]
+            + self[2][2] * self[3][3]
+            - self[0][1] * self[1][0]
+            - self[0][2] * self[2][0]
+            - self[0][3] * self[3][0]
+            - self[1][2] * self[2][1]
+            - self[1][3] * self[3][1]
+            - self[2][3] * self[3][2];
+        let poly3 = (self[0][0] * self[1][1] + self[0][0] * self[2][2] + self[1][1] * self[2][2]
+            - self[0][1] * self[1][0]
+            - self[0][2] * self[2][0]
+            - self[1][2] * self[2][1])
+            * self[3][3]
+            + (self[2][0] * self[3][2] - self[2][2] * self[3][0] + self[1][0] * self[3][1]
+                - self[1][1] * self[3][0])
+                * self[0][3]
+            + (self[2][1] * self[3][2] - self[2][2] * self[3][1] + self[3][0] * self[0][1]
+                - self[3][1] * self[0][0])
+                * self[1][3]
+            + (self[3][1] * self[1][2] - self[3][2] * self[1][1] + self[3][0] * self[0][2]
+                - self[3][2] * self[0][0])
+                * self[2][3]
+            + Matrix3::from_cols(self[0].truncate(), self[1].truncate(), self[2].truncate())
+                .determinant();
+        solver::solve_quartic(-self.trace(), poly2, -poly3, self.determinant())
+    }
+}
+
+impl<F: BaseFloat> OperatorNorm for Matrix2<F> {
+    #[inline]
+    fn norm_l1(self) -> F {
+        F::max(
+            F::abs(self[0][0]) + F::abs(self[0][1]),
+            F::abs(self[1][0]) + F::abs(self[1][1]),
+        )
+    }
+    #[inline]
+    fn norm_linf(self) -> F {
+        F::max(
+            F::abs(self[0][0]) + F::abs(self[1][0]),
+            F::abs(self[0][1]) + F::abs(self[1][1]),
+        )
+    }
+    #[inline]
+    fn norm_l2(self) -> F {
+        let eigens = (self.transpose() * self).eigenvalues();
+        F::sqrt(F::max(eigens[0].re, eigens[1].re))
+    }
+}
+
+impl<F: BaseFloat> OperatorNorm for Matrix3<F> {
+    #[inline]
+    fn norm_l1(self) -> F {
+        F::max(
+            F::max(
+                F::abs(self[0][0]) + F::abs(self[0][1]) + F::abs(self[0][2]),
+                F::abs(self[1][0]) + F::abs(self[1][1]) + F::abs(self[1][2]),
+            ),
+            F::abs(self[2][0]) + F::abs(self[2][1]) + F::abs(self[2][2]),
+        )
+    }
+    #[inline]
+    fn norm_linf(self) -> F {
+        F::max(
+            F::max(
+                F::abs(self[0][0]) + F::abs(self[1][0]) + F::abs(self[2][0]),
+                F::abs(self[0][1]) + F::abs(self[1][1]) + F::abs(self[2][1]),
+            ),
+            F::abs(self[0][2]) + F::abs(self[1][2]) + F::abs(self[2][2]),
+        )
+    }
+    #[inline]
+    fn norm_l2(self) -> F {
+        let eigens = (self.transpose() * self).eigenvalues();
+        F::sqrt(F::max(F::max(eigens[0].re, eigens[1].re), eigens[2].re))
+    }
+}
+
+impl<F: BaseFloat> OperatorNorm for Matrix4<F> {
+    #[inline]
+    fn norm_l1(self) -> F {
+        F::max(
+            F::max(
+                F::max(
+                    F::abs(self[0][0])
+                        + F::abs(self[0][1])
+                        + F::abs(self[0][2])
+                        + F::abs(self[0][3]),
+                    F::abs(self[1][0])
+                        + F::abs(self[1][1])
+                        + F::abs(self[1][2])
+                        + F::abs(self[1][3]),
+                ),
+                F::abs(self[2][0]) + F::abs(self[2][1]) + F::abs(self[2][2]) + F::abs(self[2][3]),
+            ),
+            F::abs(self[3][0]) + F::abs(self[3][1]) + F::abs(self[3][2]) + F::abs(self[3][3]),
+        )
+    }
+    #[inline]
+    fn norm_linf(self) -> F {
+        F::max(
+            F::max(
+                F::max(
+                    F::abs(self[0][0])
+                        + F::abs(self[1][0])
+                        + F::abs(self[2][0])
+                        + F::abs(self[3][0]),
+                    F::abs(self[0][1])
+                        + F::abs(self[1][1])
+                        + F::abs(self[2][1])
+                        + F::abs(self[3][1]),
+                ),
+                F::abs(self[0][2]) + F::abs(self[1][2]) + F::abs(self[2][2]) + F::abs(self[3][2]),
+            ),
+            F::abs(self[0][3]) + F::abs(self[1][3]) + F::abs(self[2][3]) + F::abs(self[3][3]),
+        )
+    }
+    #[inline]
+    fn norm_l2(self) -> F {
+        let eigens = (self.transpose() * self).eigenvalues();
+        F::sqrt(F::max(
+            F::max(F::max(eigens[0].re, eigens[1].re), eigens[2].re),
+            eigens[3].re,
+        ))
+    }
+}

vendor/matext4cgmath/src/exp_decomp.rs

@@ -0,0 +1,106 @@
+use crate::*;
+
+impl<F: BaseFloat> Exponential for Matrix2<F> {}
+impl<F: BaseFloat> Exponential for Matrix3<F> {}
+impl<F: BaseFloat> Exponential for Matrix4<F> {}
+
+impl<F: BaseFloat> Decomposition for Matrix2<F> {
+    fn iwasawa_decomposition(self) -> Option<(Self, Self, Self)> {
+        let v0 = self[0];
+        let mag0 = dot(v0, v0);
+        if mag0 <= F::zero() {
+            return None;
+        }
+        let n0 = dot(self[1], v0) / mag0;
+        let v1 = self[1] - v0 * n0;
+        let mag1 = dot(v1, v1);
+        if mag1 <= F::zero() {
+            return None;
+        }
+        let (a0, a1) = (F::sqrt(mag0), F::sqrt(mag1));
+        Some((
+            Matrix2::from_cols(v0 / a0, v1 / a1),
+            Matrix2::from_diagonal(Vector2::new(a0, a1)),
+            Matrix2::new(F::one(), F::zero(), n0, F::one()),
+        ))
+    }
+}
+
+impl<F: BaseFloat> Decomposition for Matrix3<F> {
+    fn iwasawa_decomposition(self) -> Option<(Self, Self, Self)> {
+        let v0 = self[0];
+        let mag0 = dot(v0, v0);
+        if mag0 <= F::zero() {
+            return None;
+        }
+        let n0 = dot(self[1], v0) / mag0;
+        let v1 = self[1] - v0 * n0;
+        let mag1 = dot(v1, v1);
+        if mag1 <= F::zero() {
+            return None;
+        }
+        let n1 = dot(self[2], v0) / mag0;
+        let n2 = dot(self[2], v1) / mag1;
+        let v2 = self[2] - v0 * n1 - v1 * n2;
+        let mag2 = dot(v2, v2);
+        if mag2 <= F::zero() {
+            return None;
+        }
+        let (a0, a1, a2) = (F::sqrt(mag0), F::sqrt(mag1), F::sqrt(mag2));
+        #[rustfmt::skip]
+        let res = Some((
+            Matrix3::from_cols(v0 / a0, v1 / a1, v2 / a2),
+            Matrix3::from_diagonal(Vector3::new(a0, a1, a2)),
+            Matrix3::new(
+                F::one(), F::zero(), F::zero(),
+                n0, F::one(), F::zero(),
+                n1, n2, F::one(),
+            ),
+        ));
+        res
+    }
+}
+
+impl<F: BaseFloat> Decomposition for Matrix4<F> {
+    fn iwasawa_decomposition(self) -> Option<(Self, Self, Self)> {
+        let v0 = self[0];
+        let mag0 = dot(v0, v0);
+        if mag0 <= F::zero() {
+            return None;
+        }
+        let n0 = dot(self[1], v0) / mag0;
+        let v1 = self[1] - v0 * n0;
+        let mag1 = dot(v1, v1);
+        if mag1 <= F::zero() {
+            return None;
+        }
+        let n1 = dot(self[2], v0) / mag0;
+        let n2 = dot(self[2], v1) / mag1;
+        let v2 = self[2] - v0 * n1 - v1 * n2;
+        let mag2 = dot(v2, v2);
+        if mag2 <= F::zero() {
+            return None;
+        }
+        let n3 = dot(self[3], v0) / mag0;
+        let n4 = dot(self[3], v1) / mag1;
+        let n5 = dot(self[3], v2) / mag2;
+        let v3 = self[3] - v0 * n3 - v1 * n4 - v2 * n5;
+        let mag3 = dot(v3, v3);
+        if mag3 <= F::zero() {
+            return None;
+        }
+        let (a0, a1, a2, a3) = (F::sqrt(mag0), F::sqrt(mag1), F::sqrt(mag2), F::sqrt(mag3));
+        #[rustfmt::skip]
+        let res = Some((
+            Matrix4::from_cols(v0 / a0, v1 / a1, v2 / a2, v3 / a3),
+            Matrix4::from_diagonal(Vector4::new(a0, a1, a2, a3)),
+            Matrix4::new(
+                F::one(), F::zero(), F::zero(), F::zero(),
+                n0, F::one(), F::zero(), F::zero(),
+                n1, n2, F::one(), F::zero(),
+                n3, n4, n5, F::one(),
+            ),
+        ));
+        res
+    }
+}

vendor/matext4cgmath/src/lib.rs

@@ -0,0 +1,120 @@
+//! Unofficial third-party cgmath extensions for calculate eigenvalues, operator norms and Iwasawa decomposition.
+
+// SolveCraft: compiler warnings no longer fail the build.
+#![deny(clippy::all, rust_2018_idioms)]
+#![warn(
+    missing_docs,
+    missing_debug_implementations,
+    trivial_casts,
+    trivial_numeric_casts,
+    unsafe_code,
+    unstable_features,
+    unused_import_braces,
+    unused_qualifications
+)]
+#![allow(clippy::many_single_char_names)]
+
+pub use cgmath;
+use cgmath::*;
+use num_complex::Complex;
+
+mod eigens;
+mod exp_decomp;
+/// solvers for low dimensional algebraic equations.
+pub mod solver;
+
+/// extension for eigen values
+pub trait EigenValues: VectorSpace {
+    /// the type of the array of eigen values
+    type EigenValues;
+    /// calculate eigen values.
+    ///
+    /// # Examples
+    ///
+    /// ```
+    /// use cgmath::*;
+    /// use matext4cgmath::*;
+    /// use num_complex::Complex;
+    /// const EPS: f64 = 1.0e-10;
+    ///
+    /// let mat = Matrix2::new(4.0, -2.0, 3.0, -1.0);
+    /// let mut eigens = mat.eigenvalues();
+    /// // Even in the case of real solutions, the order in the array is not guaranteed.
+    /// eigens.sort_by(|x, y| x.re.partial_cmp(&y.re).unwrap());
+    /// assert!(Complex::norm(eigens[0] - 1.0) < EPS);
+    /// assert!(Complex::norm(eigens[1] - 2.0) < EPS);
+    /// ```
+    fn eigenvalues(self) -> Self::EigenValues;
+}
+
+#[cfg_attr(doc, katexit::katexit)]
+/// [operator norms](https://en.wikipedia.org/wiki/Matrix_norm): $L^1$, $L^2$, and $L^\infty$.
+/// # Examples
+///
+/// ```
+/// use cgmath::*;
+/// use matext4cgmath::*;
+///
+/// let mat = Matrix2::new(1.0, 3.0, -2.0, 4.0);
+/// // L^1 operator norm is the maximum column absolute sumation.
+/// assert_eq!(mat.norm_l1(), 6.0);
+/// // L^2 operator norm is the maximum singular value.
+/// let ans_norm_l2 = (5.0 + f64::sqrt(5.0)) / f64::sqrt(2.0);
+/// assert!(f64::abs(mat.norm_l2() - ans_norm_l2) < 1.0e-10);
+/// // L^∞ operator norm is the maximum row absolute sumation.
+/// assert_eq!(mat.norm_linf(), 7.0);
+/// ```
+pub trait OperatorNorm: VectorSpace {
+    /// operator norm for absolute value sumation: $L^1$.
+    fn norm_l1(self) -> Self::Scalar;
+    /// operator norm for Euclidean norm: $L^2$.
+    fn norm_l2(self) -> Self::Scalar;
+    /// operator norm for absolute value maximum: $L^{\infty}$.
+    fn norm_linf(self) -> Self::Scalar;
+}
+
+/// calculate exponential value
+pub trait Exponential: OperatorNorm + std::ops::AddAssign<Self> + One
+where
+    Self::Scalar: BaseFloat,
+{
+    /// calculate exponential
+    ///
+    /// # Examples
+    ///
+    /// ```
+    /// use cgmath::*;
+    /// use matext4cgmath::*;
+    ///
+    /// let x = Matrix2::new(0.0, 1.0, -1.0, 0.0);
+    /// let res = x.exp();
+    /// let ans = Matrix2::from_angle(Rad(1.0));
+    /// assert!((res - ans).norm_l1() < 1.0e-10);
+    /// ```
+    fn exp(self) -> Self {
+        use num_traits::{Float, NumCast};
+        let eps = <Self::Scalar as Float>::epsilon();
+        let mut a = self;
+        let mut res = Self::one();
+        for i in 2_u16..=64 {
+            res += a;
+            a = a * self / <Self::Scalar as NumCast>::from(i).unwrap();
+            if a.norm_linf() < eps * self.norm_linf() {
+                break;
+            }
+        }
+        res
+    }
+}
+
+#[cfg_attr(doc, katexit::katexit)]
+/// some decompositions of matrix
+pub trait Decomposition: VectorSpace {
+    #[cfg_attr(doc, katexit::katexit)]
+    /// Returns $(K, A, N)$ of the Iwasawa decomposition: $M = KAN$.
+    ///
+    /// - $K$: orthonormal matrix
+    /// - $A$: diagonal matrix
+    /// - $N$: upper-half unipotent matrix
+    fn iwasawa_decomposition(self) -> Option<(Self, Self, Self)>;
+}

vendor/matext4cgmath/src/solver.rs

@@ -0,0 +1,246 @@
+use crate::*;
+
+#[cfg_attr(doc, katexit::katexit)]
+/// Solve equation: $x^2 + ax + b = 0$.
+///
+/// # Examples
+///
+/// ```
+/// use matext4cgmath::solver;
+/// use num_complex::Complex;
+/// const EPS: f64 = 1.0e-10;
+///
+/// let res = solver::solve_quadratic(3.0, -4.0);
+/// // Real solutions are always sorted in order of solution size.
+/// assert!(Complex::norm(res[0] + 4.0) < EPS);
+/// assert!(Complex::norm(res[1] - 1.0) < EPS);
+/// ```
+pub fn solve_quadratic<F: BaseFloat>(a: F, b: F) -> [Complex<F>; 2] {
+    let two = F::one() + F::one();
+    let four = two + two;
+    let det = a * a - four * b;
+    match det >= F::zero() {
+        true => {
+            let h = F::sqrt(det);
+            [
+                Complex::new((-a - h) / two, F::zero()),
+                Complex::new((-a + h) / two, F::zero()),
+            ]
+        }
+        false => {
+            let h = F::sqrt(-det);
+            [Complex::new(-a, h) / two, Complex::new(-a, -h) / two]
+        }
+    }
+}
+
+#[cfg_attr(doc, katexit::katexit)]
+/// Solve equation: $x^3 + px + q = 0$.
+///
+/// # Examples
+///
+/// ```
+/// use matext4cgmath::solver;
+/// use num_complex::Complex;
+/// const EPS: f64 = 1.0e-10;
+///
+/// let mut res = solver::pre_solve_cubic(-7.0, -6.0);
+/// // Even in the case of real solutions, the order in the array is not guaranteed.
+/// res.sort_by(|x, y| x.re.partial_cmp(&y.re).unwrap());
+/// let ans = [Complex::from(-2.0), Complex::from(-1.0), Complex::from(3.0)];
+/// res.iter().zip(ans).for_each(|(x, y)| {
+///     assert!(Complex::norm(x - y) < EPS);
+/// });
+/// ```
+pub fn pre_solve_cubic<F: BaseFloat>(p: F, q: F) -> [Complex<F>; 3] {
+    let two = F::one() + F::one();
+    let three = two + F::one();
+    let sqrt3_2 = F::sqrt(three) / two;
+    let omega = Complex::new(-F::one() / two, sqrt3_2);
+    let omega2 = Complex::new(-F::one() / two, -sqrt3_2);
+    let eps_2 = F::sqrt(F::epsilon());
+
+    let p_3 = p / three;
+    let q_2 = q / two;
+    let alpha2 = q_2 * q_2 + p_3 * p_3 * p_3;
+    let (x, y) = match alpha2 >= F::zero() {
+        true => {
+            let alpha = F::sqrt(alpha2);
+            let tmpx = -q_2 - alpha;
+            let tmpy = -q_2 + alpha;
+            (
+                Complex::new(
+                    F::signum(tmpx) * F::powf(F::abs(tmpx), F::one() / three),
+                    F::zero(),
+                ),
+                Complex::new(
+                    F::signum(tmpy) * F::powf(F::abs(tmpy), F::one() / three),
+                    F::zero(),
+                ),
+            )
+        }
+        false => {
+            let alphai = F::sqrt(-alpha2);
+            (
+                Complex::powf(Complex::new(-q_2, alphai), F::one() / three),
+                Complex::powf(Complex::new(-q_2, -alphai), F::one() / three),
+            )
+        }
+    };
+    let mut res = [x + y, omega * x + omega2 * y, omega2 * x + omega * y];
+    // precision by Newton method
+    res.iter_mut().for_each(|x| {
+        let mut f = *x * *x * *x + *x * p + q;
+        let mut f_prime = *x * *x * three + p;
+        // SolveCraft: at most 100 steps. Near a double root (the eigenvalues of a
+        // translation matrix) the iteration can cycle forever without meeting the tolerance.
+        let mut steps = 0;
+        while f.norm() > eps_2 * f_prime.norm() && steps < 100 {
+            steps += 1;
+            if f_prime.norm() < eps_2 {
+                return;
+            }
+            *x -= f / f_prime;
+            f = *x * *x * *x + *x * p + q;
+            f_prime = *x * *x * three + p;
+        }
+    });
+    res
+}
+
+#[cfg_attr(doc, katexit::katexit)]
+/// solve equation: $x^3 + ax^2 + bx + c = 0$.
+///
+/// # Examples
+///
+/// ```
+/// use matext4cgmath::solver;
+/// use num_complex::Complex;
+/// const EPS: f64 = 1.0e-10;
+///
+/// let mut res = solver::solve_cubic(-3.0, 0.0, 4.0);
+/// // Even in the case of real solutions, the order in the array is not guaranteed.
+/// res.sort_by(|x, y| x.re.partial_cmp(&y.re).unwrap());
+/// let ans = [Complex::from(-1.0), Complex::from(2.0), Complex::from(2.0)];
+/// res.iter().zip(ans).for_each(|(x, y)| {
+///     assert!(Complex::norm(x - y) < EPS);
+/// });
+/// ```
+pub fn solve_cubic<F: BaseFloat>(a: F, b: F, c: F) -> [Complex<F>; 3] {
+    let two = F::one() + F::one();
+    let three = two + F::one();
+    let twenty_seven = three * three * three;
+    let p = b - a * a / three;
+    let q = c - a * b / three + two * a * a * a / twenty_seven;
+    let mut res = pre_solve_cubic(p, q);
+    res.iter_mut().for_each(|x| {
+        *x -= a / three;
+    });
+    res
+}
+
+#[cfg_attr(doc, katexit::katexit)]
+/// solve equation: $x^4 + px^2 + qx + r = 0$.
+///
+/// # Examples
+///
+/// ```
+/// use matext4cgmath::solver;
+/// use num_complex::Complex;
+/// const EPS: f64 = 1e-7;
+///
+/// let mut res = solver::pre_solve_quartic(-5.0, 0.0, 4.0);
+/// // Even in the case of real solutions, the order in the array is not guaranteed.
+/// res.sort_by(|x, y| x.re.partial_cmp(&y.re).unwrap());
+/// let ans = [Complex::from(-2.0), Complex::from(-1.0), Complex::from(1.0), Complex::from(2.0)];
+/// res.iter().zip(ans).for_each(|(x, y)| {
+///     assert!(Complex::norm(x - y) < EPS);
+/// });
+/// ```
+pub fn pre_solve_quartic<F: BaseFloat>(p: F, q: F, r: F) -> [Complex<F>; 4] {
+    let one = F::one();
+    let two = one + one;
+    let four = two + two;
+    let eps_2 = F::sqrt(F::epsilon());
+
+    let a = two * p;
+    let b = p * p - four * r;
+    let c = -q * q;
+    let f = solve_cubic(a, b, c);
+    let a = f[0].sqrt() / two;
+    let b = f[1].sqrt() / two;
+    let c = f[2].sqrt() / two;
+
+    let mut res = (0..8)
+        .map(|i| {
+            let a = a * F::powi(-F::one(), i % 2);
+            let b = b * F::powi(-F::one(), (i / 2) % 2);
+            let c = c * F::powi(-F::one(), (i / 4) % 2);
+            [-a - b - c, -a + b + c, a - b + c, a + b - c]
+        })
+        .map(|x| {
+            let f = x
+                .iter()
+                .map(|t| (t * t * t * t + t * t * p + t * q + r).norm_sqr())
+                .max_by(|x, y| x.partial_cmp(y).unwrap())
+                .unwrap();
+            (x, f)
+        })
+        .min_by(|x, y| x.1.partial_cmp(&y.1).unwrap())
+        .unwrap()
+        .0;
+    // refinement by Newton method
+    res.iter_mut().for_each(|x| {
+        let mut f = *x * *x * *x * *x + *x * *x * p + *x * q + r;
+        let mut f_prime = *x * *x * *x * four + *x * p * two + q;
+        // SolveCraft: at most 100 steps. Near a double root (the eigenvalues of a
+        // translation matrix) the iteration can cycle forever without meeting the tolerance.
+        let mut steps = 0;
+        while f.norm() > eps_2 * f_prime.norm() && steps < 100 {
+            steps += 1;
+            if f_prime.norm() < eps_2 {
+                return;
+            }
+            *x -= f / f_prime;
+            f_prime = *x * *x * *x * four + *x * p * two + q;
+            f = *x * *x * *x * *x + *x * *x * p + *x * q + r;
+        }
+    });
+    res
+}
+
+#[cfg_attr(doc, katexit::katexit)]
+/// solve equation: $x^4 + ax^3 + bx^2 + cx + d = 0$.
+///
+/// # Examples
+///
+/// ```
+/// use matext4cgmath::solver;
+/// use num_complex::Complex;
+/// const EPS: f64 = 1.0e-10;
+///
+/// let mut res = solver::solve_quartic(1.0, -7.0, -1.0, 6.0);
+/// // Even in the case of real solutions, the order in the array is not guaranteed.
+/// res.sort_by(|x, y| x.re.partial_cmp(&y.re).unwrap());
+/// let ans = [Complex::from(-3.0), Complex::from(-1.0), Complex::from(1.0), Complex::from(2.0)];
+/// res.iter().zip(ans).for_each(|(x, y)| {
+///     assert!(Complex::norm(x * x * x * x + x * x * x - 7.0 * x * x - x + 6.0) < EPS);
+///     assert!(Complex::norm(x - y) < EPS, "{x} {y}");
+/// });
+/// ```
+pub fn solve_quartic<F: BaseFloat>(a: F, b: F, c: F, d: F) -> [Complex<F>; 4] {
+    let one = F::one();
+    let two = one + one;
+    let three = one + two;
+    let four = two + two;
+    let six = two * three;
+    let eight = four + four;
+
+    let a_4 = a / four;
+    let p = b - six * a_4 * a_4;
+    let q = c - two * b * a_4 + eight * a_4 * a_4 * a_4;
+    let r = d - c * a_4 + b * a_4 * a_4 - three * a_4 * a_4 * a_4 * a_4;
+    let mut res = pre_solve_quartic(p, q, r);
+    res.iter_mut().for_each(|x| *x -= a_4);
+    res
+}