diff --git a/Cargo.lock b/Cargo.lock
index e91aaca24..c3770a76d 100644
--- a/Cargo.lock
+++ b/Cargo.lock
@@ -4417,6 +4417,7 @@ name = "solidcraft-smlib"
 version = "0.18.1"
 dependencies = [
  "cc",
+ "flate2",
  "log",
  "serde_json",
  "solidcraft-geom",
diff --git a/crates/smlib/Cargo.toml b/crates/smlib/Cargo.toml
index 5981899a3..c7f117395 100644
--- a/crates/smlib/Cargo.toml
+++ b/crates/smlib/Cargo.toml
@@ -32,4 +32,5 @@ solidcraft-geom = { workspace = true }
 cc = { version = "1.2", features = ["parallel"] }
 
 [dev-dependencies]
+flate2 = "1.1"
 serde_json = { workspace = true }
diff --git a/crates/smlib/tests/ssi_corpus.rs b/crates/smlib/tests/ssi_corpus.rs
new file mode 100644
index 000000000..9b387a72d
--- /dev/null
+++ b/crates/smlib/tests/ssi_corpus.rs
@@ -0,0 +1,360 @@
+//! The SSI corpus (`tests/ssi-corpus`) and its scorer (`ssi_score`, a port of the corpus's
+//! `scoring/score.py`).
+//!
+//! The parity check against the Python scorer is ignored by default; it needs the corpus
+//! runner's raw outputs and score.py's scores for the same run:
+//!
+//! ```text
+//! SSI_PARITY_RAW=<corpus>/results/raw/gsi_1e-6 \
+//! SSI_PARITY_SCORES=<corpus>/results/score_gsi_1e-6.json \
+//! SSI_PARITY_TOL=1e-6 SSI_SM_MESSAGES=<smlib>/source/SMLib/inc/SmMessages.h \
+//!   cargo test -p solidcraft-smlib --test ssi_corpus -- --ignored --nocapture
+//! ```
+
+mod ssi_score;
+
+use serde_json::Value;
+use ssi_score::geom::{Nurbs, V3, linspace, np_max};
+use ssi_score::{Case, Coverage, EngineCurve, EngineResult, Outcome, Score, Tolerances, Verdict, case_ids, fmt_e, load_case, score_case};
+
+fn all_cases() -> Vec<Case> {
+    let ids = case_ids();
+    assert_eq!(ids.len(), 119, "corpus case count");
+    ids.iter().map(|id| load_case(id).unwrap_or_else(|e| panic!("{id}: {e}"))).collect()
+}
+
+/// `f` over `items` on a few threads, results in order.
+fn par_map<T: Sync, R: Send>(items: &[T], f: impl Fn(&T) -> R + Sync) -> Vec<R> {
+    let n = std::thread::available_parallelism().map_or(1, |n| n.get()).clamp(1, 8);
+    let next = std::sync::atomic::AtomicUsize::new(0);
+    let mut out: Vec<(usize, R)> = std::thread::scope(|s| {
+        let workers: Vec<_> = (0..n)
+            .map(|_| {
+                s.spawn(|| {
+                    let mut mine = Vec::new();
+                    loop {
+                        let i = next.fetch_add(1, std::sync::atomic::Ordering::Relaxed);
+                        let Some(item) = items.get(i) else { break };
+                        mine.push((i, f(item)));
+                    }
+                    mine
+                })
+            })
+            .collect();
+        workers.into_iter().flat_map(|w| w.join().unwrap()).collect()
+    });
+    out.sort_by_key(|(i, _)| *i);
+    out.into_iter().map(|(_, r)| r).collect()
+}
+
+/// The expected intersection itself, as an engine answer.
+fn truth_as_output(case: &Case) -> EngineResult {
+    let curve = |p: &Vec<V3>| EngineCurve { curve_type: "CROSSING".into(), xyz: p.clone(), ..EngineCurve::default() };
+    EngineResult {
+        status: ssi_score::SM_SUCCESS,
+        status_name: None,
+        time_ms: 0.0,
+        curves: case.expected.branches.iter().chain(&case.expected.coincident_regions).map(curve).collect(),
+        points: case.expected.isolated_points.clone(),
+    }
+}
+
+#[test]
+fn every_case_parses_and_cites_its_source() {
+    let cases = all_cases();
+    let y23: Vec<&Case> = cases.iter().filter(|c| c.id.starts_with("y23_")).collect();
+    assert_eq!(y23.len(), 25);
+    for c in &y23 {
+        let src = c.json.get("source").and_then(Value::as_str).unwrap_or_default();
+        assert!(src.contains("Yang, Jia & Yan 2023"), "{}: source {src:?}", c.id);
+    }
+    for c in &cases {
+        assert_eq!(c.surfaces.len(), 2, "{}", c.id);
+        assert!(c.cov_tol > 0.0, "{}", c.id);
+    }
+}
+
+/// Every NURBS surface of the corpus lies on its exact truth (the corpus's own self-check): this
+/// checks the NURBS evaluation and the distance functions against each other.
+#[test]
+fn nurbs_surfaces_lie_on_their_truths() {
+    let mut checked = 0;
+    for case in all_cases() {
+        for (k, spec) in case.surfaces.iter().enumerate() {
+            if spec.get("type").and_then(Value::as_str) != Some("nurbs") {
+                continue;
+            }
+            let s = Nurbs::from_spec(spec).unwrap_or_else(|e| panic!("{}: {e}", case.id));
+            let d = s.domain();
+            let pts: Vec<V3> = linspace(d[0], d[1], 13)
+                .into_iter()
+                .flat_map(|u| linspace(d[2], d[3], 13).into_iter().map(move |v| (u, v)))
+                .map(|(u, v)| s.eval(u, v))
+                .collect();
+            let err = np_max(&case.truth[k].distance(&pts));
+            assert!(err < 1e-9, "{} surface {k}: {err:e} off its truth", case.id);
+            checked += 1;
+        }
+    }
+    assert!(checked >= 99, "{checked} NURBS surfaces checked");
+}
+
+/// Handing the scorer the expected intersection itself passes every case, except that the
+/// expected branches run through singular points where an answer must split its curves.
+#[test]
+fn the_truth_itself_passes() {
+    let mut fails = Vec::new();
+    let mut unsplit = 0;
+    let cases = all_cases();
+    let scores = par_map(&cases, |case| score_case(case, &Outcome::Ran(truth_as_output(case)), Tolerances::at(1e-6)));
+    for (case, s) in cases.iter().zip(scores) {
+        let only_unsplit = !case.expected.singular_points.is_empty()
+            && !s.flags.is_empty()
+            && s.flags.iter().all(|f| f.starts_with("singular point ") && f.contains(" not identified "));
+        if s.verdict == Verdict::WrongTopology && only_unsplit {
+            unsplit += 1;
+        } else if s.verdict != Verdict::Pass {
+            fails.push(format!("{} {}: {}", case.id, s.verdict, s.flags.join("; ")));
+        }
+    }
+    assert!(fails.is_empty(), "{}", fails.join("\n"));
+    assert_eq!(unsplit, 10, "cases whose expected branches run through a singular point");
+}
+
+#[test]
+fn verdicts_follow_severity() {
+    let case = load_case("c1_cone_plane_ellipse").unwrap();
+    let tol = Tolerances::at(1e-6);
+    let truth = truth_as_output(&case);
+
+    let crash = score_case(&case, &Outcome::Crashed("killed by signal 11".into()), tol);
+    assert_eq!(crash.verdict, Verdict::Crash);
+    assert_eq!(crash.detail.as_deref(), Some("killed by signal 11"));
+
+    let mut failed = truth.clone();
+    failed.status = 1001;
+    let s = score_case(&case, &Outcome::Ran(failed), tol);
+    assert_eq!(s.verdict, Verdict::Error);
+    assert_eq!(s.flags[0], "ERROR status 1001");
+
+    let empty = EngineResult { curves: Vec::new(), ..truth.clone() };
+    let s = score_case(&case, &Outcome::Ran(empty), tol);
+    assert_eq!(s.verdict, Verdict::Incomplete);
+    assert!(s.flags[0].starts_with("branch 0 covered 0% (max gap inf"), "{:?}", s.flags);
+
+    // half of the ellipse: covered in part, two free ends
+    let mut half = truth.clone();
+    let n = half.curves[0].xyz.len();
+    half.curves[0].xyz.truncate(n / 2);
+    assert_eq!(score_case(&case, &Outcome::Ran(half), tol).verdict, Verdict::Incomplete);
+
+    // the ellipse twice: duplicated output
+    let mut twice = truth.clone();
+    twice.curves.push(twice.curves[0].clone());
+    let s = score_case(&case, &Outcome::Ran(twice), tol);
+    assert_eq!(s.verdict, Verdict::WrongTopology);
+    assert!(s.flags.iter().any(|f| f == "output curves 0 and 1 overlap (100% of one lies on the other)"), "{:?}", s.flags);
+
+    // a stray point far away
+    let mut stray = truth.clone();
+    stray.points.push([50.0, 50.0, 50.0]);
+    assert_eq!(score_case(&case, &Outcome::Ran(stray), tol).verdict, Verdict::WrongTopology);
+
+    // the ellipse lifted by 1e-4: still covered (cov_tol 1e-3) but past acc_tol = 1e-5
+    let mut lifted = truth.clone();
+    for p in &mut lifted.curves[0].xyz {
+        p[2] += 1e-4;
+    }
+    let s = score_case(&case, &Outcome::Ran(lifted), tol);
+    assert_eq!(s.verdict, Verdict::Inaccurate, "{:?}", s.flags);
+    assert!(s.max_dist_3d > 5e-5 && s.max_dist_3d < 1e-4 + 1e-12, "{}", s.max_dist_3d);
+
+    // a TOUCHING curve counts as a point
+    let mut touching = truth.clone();
+    touching.curves[0].curve_type = "TOUCHING".into();
+    let s = score_case(&case, &Outcome::Ran(touching), tol);
+    assert_eq!(s.n_out_points, 1);
+    assert_eq!(s.verdict, Verdict::Incomplete);
+
+    assert_eq!(score_case(&case, &Outcome::Ran(truth), tol).verdict, Verdict::Pass);
+}
+
+#[test]
+fn python_number_formatting() {
+    assert_eq!(fmt_e(1.2345e-5, 2), "1.23e-05");
+    assert_eq!(fmt_e(1e-5, 0), "1e-05");
+    assert_eq!(fmt_e(123.0, 1), "1.2e+02");
+    assert_eq!(fmt_e(f64::INFINITY, 2), "inf");
+    assert_eq!(fmt_e(3.0e-300, 2), "3.00e-300");
+}
+
+// ------------------------------------------------------------------------------------------
+// parity with the Python scorer
+
+/// Python's `json.dump` writes `Infinity`, `-Infinity` and `NaN` bare; quote them (outside
+/// strings) so serde_json reads them, and read numbers through [`pnum`].
+fn python_json(s: &str) -> Value {
+    let mut out = String::with_capacity(s.len());
+    let (mut in_str, mut esc) = (false, false);
+    let mut rest = s;
+    while let Some(c) = rest.chars().next() {
+        if !in_str && let Some(tok) = ["-Infinity", "Infinity", "NaN"].into_iter().find(|t| rest.starts_with(t)) {
+            out.push_str(&format!("\"{tok}\""));
+            rest = &rest[tok.len()..];
+            continue;
+        }
+        if in_str {
+            if esc {
+                esc = false;
+            } else if c == '\\' {
+                esc = true;
+            } else if c == '"' {
+                in_str = false;
+            }
+        } else if c == '"' {
+            in_str = true;
+        }
+        out.push(c);
+        rest = &rest[c.len_utf8()..];
+    }
+    serde_json::from_str(&out).expect("score JSON")
+}
+
+fn pnum(v: Option<&Value>) -> f64 {
+    match v {
+        Some(Value::String(s)) => s.parse().unwrap_or(f64::NAN),
+        Some(v) => v.as_f64().unwrap_or(f64::NAN),
+        None => f64::NAN,
+    }
+}
+
+fn close(a: f64, b: f64) -> bool {
+    // Closest-point distances that converge to round-off (~1e-12) differ in their last digits.
+    (a.is_nan() && b.is_nan()) || a == b || (a - b).abs() <= 1e-9 * a.abs().max(b.abs()) + 1e-11
+}
+
+/// Compares the numbers of one score with score.py's record; returns the mismatches.
+fn compare(rs: &Score, py: &Value) -> Vec<String> {
+    let mut bad = Vec::new();
+    let f = |k: &str| pnum(py.get(k));
+    let mut num = |k: &str, v: f64| {
+        let p = f(k);
+        if !close(v, p) {
+            bad.push(format!("{k}: rust {v:e} python {p:e}"));
+        }
+    };
+    num("cov_tol", rs.cov_tol);
+    if rs.verdict != Verdict::Crash {
+        num("max_dist_3d", rs.max_dist_3d);
+        num("max_dist_uv", rs.max_dist_uv);
+        num("hausdorff_truth_to_out", rs.hausdorff_truth_to_out);
+        num("hausdorff_out_to_truth", rs.hausdorff_out_to_truth);
+        let cov = py.get("coverage").and_then(Value::as_array).cloned().unwrap_or_default();
+        if cov.len() != rs.coverage.len() {
+            bad.push(format!("coverage: {} entries, python {}", rs.coverage.len(), cov.len()));
+        }
+        for (c, p) in rs.coverage.iter().zip(&cov) {
+            let g = |k: &str| pnum(p.get(k));
+            let (a, b) = match c {
+                Coverage::Polyline { max, covered_frac } => ((*max, g("max")), (*covered_frac, g("covered_frac"))),
+                Coverage::Point { dist, .. } => ((*dist, g("dist")), (0.0, 0.0)),
+            };
+            if !close(a.0, a.1) || !close(b.0, b.1) {
+                bad.push(format!("coverage {c:?} python {p}"));
+            }
+        }
+        let strs = |k: &str| -> Vec<String> {
+            py.get(k).and_then(Value::as_array).map(|a| a.iter().filter_map(|s| s.as_str().map(str::to_string)).collect()).unwrap_or_default()
+        };
+        if strs("expected_topology") != rs.expected_topology || strs("output_topology") != rs.output_topology {
+            bad.push(format!(
+                "topology {:?}/{:?} python {:?}/{:?}",
+                rs.output_topology,
+                rs.expected_topology,
+                strs("output_topology"),
+                strs("expected_topology")
+            ));
+        }
+        if strs("types") != rs.types {
+            bad.push("types".into());
+        }
+        let n = |k: &str| py.get(k).and_then(Value::as_u64).unwrap_or(u64::MAX) as usize;
+        if n("n_curves") != rs.n_curves || n("n_out_points") != rs.n_out_points {
+            bad.push("curve/point counts".into());
+        }
+    }
+    let flags: Vec<String> =
+        py.get("flags").and_then(Value::as_array).map(|a| a.iter().filter_map(|s| s.as_str().map(str::to_string)).collect()).unwrap_or_default();
+    if flags != rs.flags {
+        bad.push(format!("flags: rust {:?} python {:?}", rs.flags, flags));
+    }
+    bad
+}
+
+#[test]
+#[ignore = "needs the corpus runner's raw outputs and score.py's scores (SSI_PARITY_RAW, SSI_PARITY_SCORES)"]
+fn parity_with_python_scorer() {
+    let raw = std::env::var("SSI_PARITY_RAW").expect("SSI_PARITY_RAW: the runner's raw output directory");
+    let scores = std::env::var("SSI_PARITY_SCORES").expect("SSI_PARITY_SCORES: score.py's score JSON");
+    let tol: f64 = std::env::var("SSI_PARITY_TOL").map(|s| s.parse().expect("SSI_PARITY_TOL")).unwrap_or(1e-6);
+    let names = std::env::var("SSI_SM_MESSAGES")
+        .ok()
+        .map(|p| ssi_score::status_names(&std::fs::read_to_string(p).expect("SSI_SM_MESSAGES")))
+        .unwrap_or_default();
+    let py = python_json(&std::fs::read_to_string(&scores).expect("scores"));
+    let py = py.as_array().expect("score list");
+    assert_eq!(py.len(), 119, "score.py scored {} cases", py.len());
+
+    let scored = par_map(py, |p| {
+        let id = p["id"].as_str().expect("id");
+        let case = load_case(id).expect("case");
+        let outcome = if p["verdict"] == "CRASH" {
+            Outcome::Crashed(p["detail"].as_str().unwrap_or_default().to_string())
+        } else {
+            let path = std::path::Path::new(&raw).join(format!("{id}.json"));
+            let o: Value = serde_json::from_str(&std::fs::read_to_string(&path).expect("raw output")).expect("raw JSON");
+            let mut r = EngineResult::from_runner_json(&o).expect("runner JSON");
+            if r.status_name.is_none() {
+                r.status_name = names.get(&r.status).cloned();
+            }
+            Outcome::Ran(r)
+        };
+        let t0 = std::time::Instant::now();
+        let rs = score_case(&case, &outcome, Tolerances::at(tol));
+        (rs, t0.elapsed().as_millis())
+    });
+
+    let mut same = 0;
+    let mut problems = Vec::new();
+    let delta = |a: f64, b: f64| if a == b || (a.is_nan() && b.is_nan()) { 0.0 } else { (a - b).abs() };
+    println!("{:<40} {:<15} {:<15} |rust - python|: max_dist_3d, max_dist_uv, hausdorff both ways", "case", "python", "rust");
+    for (p, (rs, ms)) in py.iter().zip(&scored) {
+        let id = rs.id.as_str();
+        let pv = p["verdict"].as_str().unwrap_or_default();
+        let g = |k: &str| pnum(p.get(k));
+        let deltas = [
+            delta(rs.max_dist_3d, g("max_dist_3d")),
+            delta(rs.max_dist_uv, g("max_dist_uv")),
+            delta(rs.hausdorff_truth_to_out, g("hausdorff_truth_to_out")),
+            delta(rs.hausdorff_out_to_truth, g("hausdorff_out_to_truth")),
+        ];
+        println!(
+            "{id:<40} {pv:<15} {:<15} {} {} {} {}  ({ms} ms)",
+            rs.verdict.as_str(),
+            fmt_e(deltas[0], 1),
+            fmt_e(deltas[1], 1),
+            fmt_e(deltas[2], 1),
+            fmt_e(deltas[3], 1)
+        );
+        if pv == rs.verdict.as_str() {
+            same += 1;
+        } else {
+            problems.push(format!("{id}: verdict rust {} python {pv}", rs.verdict));
+        }
+        for b in compare(rs, p) {
+            problems.push(format!("{id}: {b}"));
+        }
+    }
+    println!("{same}/119 identical verdicts");
+    assert!(problems.is_empty(), "{}", problems.join("\n"));
+}
diff --git a/crates/smlib/tests/ssi_score/mod.rs b/crates/smlib/tests/ssi_score/mod.rs
new file mode 100644
index 000000000..10c81201b
--- /dev/null
+++ b/crates/smlib/tests/ssi_score/mod.rs
@@ -0,0 +1,574 @@
+//! The SSI corpus scorer (`tests/ssi-corpus`), ported from the corpus's `scoring/score.py`
+//! (jyviko/ssi-corpus b04c8a6, MIT OR Apache-2.0) as test-support code.
+//!
+//! [`score_case`] grades one engine result against a case's exact truth and returns the
+//! verdict, most severe first: CRASH, ERROR, INCOMPLETE, WRONG TOPOLOGY, INACCURATE, PASS, with
+//! every flag and number score.py reports. The flag texts are score.py's, formatted the same.
+//!
+//! Fill an [`EngineResult`] from sampled curves (per sample: `xyz`, `uv1`, `uv2`, both surfaces
+//! evaluated at the UV samples, and the curve type), isolated points and the status, or read
+//! the corpus runner's JSON with [`EngineResult::from_runner_json`].
+#![allow(dead_code)]
+
+pub mod geom;
+pub mod polyset;
+
+use std::path::{Path, PathBuf};
+
+use serde_json::Value;
+
+use geom::{Truth, V3, np_max, py_max, v3, v3_list};
+use polyset::{PolySet, diameter, signature, stride};
+
+/// SMLib's `SM_SUCCESS`: the only status score.py treats as success.
+pub const SM_SUCCESS: i64 = 1000;
+
+/// The corpus directory: `tests/ssi-corpus/cases`.
+pub fn cases_dir() -> PathBuf {
+    Path::new(env!("CARGO_MANIFEST_DIR")).join("../../tests/ssi-corpus/cases")
+}
+
+/// Every case id, sorted.
+pub fn case_ids() -> Vec<String> {
+    let mut ids: Vec<String> = std::fs::read_dir(cases_dir())
+        .map(|rd| {
+            rd.filter_map(|e| e.ok()).filter_map(|e| e.file_name().to_str().and_then(|n| n.strip_suffix(".json.gz")).map(str::to_string)).collect()
+        })
+        .unwrap_or_default();
+    ids.sort();
+    ids
+}
+
+/// Reads `<id>.json.gz` from the corpus.
+pub fn load_case_json(id: &str) -> Result<Value, String> {
+    use std::io::Read;
+    let path = cases_dir().join(format!("{id}.json.gz"));
+    let f = std::fs::File::open(&path).map_err(|e| format!("{}: {e}", path.display()))?;
+    let mut s = String::new();
+    flate2::read::GzDecoder::new(f).read_to_string(&mut s).map_err(|e| format!("{}: {e}", path.display()))?;
+    serde_json::from_str(&s).map_err(|e| format!("{}: {e}", path.display()))
+}
+
+pub fn load_case(id: &str) -> Result<Case, String> {
+    Case::from_json(&load_case_json(id)?)
+}
+
+/// What a case expects: the true intersection, sampled.
+#[derive(Clone, Debug, Default)]
+pub struct Expected {
+    pub branches: Vec<Vec<V3>>,
+    pub isolated_points: Vec<V3>,
+    pub singular_points: Vec<V3>,
+    /// Boundary polylines of coincident regions.
+    pub coincident_regions: Vec<Vec<V3>>,
+}
+
+/// One corpus case (the fields the scorer reads, plus the raw JSON).
+#[derive(Clone, Debug)]
+pub struct Case {
+    pub id: String,
+    pub class: String,
+    /// The two surface specs as the corpus stores them (`type`: plane, cylinder, cone, sphere,
+    /// torus, nurbs); [`geom::Nurbs::from_spec`] evaluates the NURBS ones.
+    pub surfaces: Vec<Value>,
+    pub truth: Vec<Truth>,
+    pub expected: Expected,
+    pub cov_tol: f64,
+    pub ill_cond: Option<f64>,
+    pub feature_size: Option<f64>,
+    pub json: Value,
+}
+
+fn polylines(v: Option<&Value>, key: &str) -> Result<Vec<Vec<V3>>, String> {
+    match v {
+        None => Ok(Vec::new()),
+        Some(Value::Array(a)) => a.iter().map(|b| v3_list(b.get(key))).collect(),
+        Some(other) => Err(format!("not a list: {other}")),
+    }
+}
+
+impl Case {
+    pub fn from_json(c: &Value) -> Result<Case, String> {
+        let s = |k: &str| c.get(k).and_then(Value::as_str).unwrap_or_default().to_string();
+        let truth = c.get("truth").and_then(Value::as_array).ok_or("truth: not a list")?;
+        let truth: Vec<Truth> = truth.iter().map(Truth::from_json).collect::<Result<_, _>>()?;
+        if truth.len() != 2 {
+            return Err("need two truths".into());
+        }
+        let exp = c.get("expected").ok_or("missing expected")?;
+        let expected = Expected {
+            branches: polylines(exp.get("branches"), "points")?,
+            isolated_points: v3_list(exp.get("isolated_points"))?,
+            singular_points: v3_list(exp.get("singular_points"))?,
+            coincident_regions: polylines(exp.get("coincident_regions"), "boundary")?,
+        };
+        let opt = |k: &str| c.get(k).and_then(Value::as_f64);
+        Ok(Case {
+            id: s("id"),
+            class: s("class"),
+            surfaces: c.get("surfaces").and_then(Value::as_array).cloned().unwrap_or_default(),
+            truth,
+            expected,
+            cov_tol: opt("cov_tol").unwrap_or(1e-3),
+            ill_cond: opt("ill_cond"),
+            feature_size: opt("feature_size"),
+            json: c.clone(),
+        })
+    }
+}
+
+/// One output curve, sampled. Missing or non-finite coordinates are NaN.
+#[derive(Clone, Debug, Default)]
+pub struct EngineCurve {
+    /// The curve type as the corpus runner names it: `TOUCHING`, `CROSSING`, `TANGENT`,
+    /// `COINCIDENT`, `NEAR_TANGENT`, `REGION_BOUNDARY`, `UNKNOWN` (`n/a` without types). Only
+    /// `TOUCHING` changes the scoring: such a curve counts as a point.
+    pub curve_type: String,
+    /// 3D samples.
+    pub xyz: Vec<V3>,
+    /// UV samples on surface 1 and 2 (not scored themselves).
+    pub uv1: Vec<[f64; 2]>,
+    pub uv2: Vec<[f64; 2]>,
+    /// Surface 1 evaluated at `uv1` (must lie on surface 2), surface 2 at `uv2` (on surface 1).
+    pub s1_at_uv: Vec<V3>,
+    pub s2_at_uv: Vec<V3>,
+}
+
+/// An engine's answer for one case.
+#[derive(Clone, Debug, Default)]
+pub struct EngineResult {
+    /// The engine status; [`SM_SUCCESS`] is success.
+    pub status: i64,
+    /// The status' name for the flags (score.py: the runner's `status_name`, else the SMLib
+    /// name, else the number).
+    pub status_name: Option<String>,
+    pub time_ms: f64,
+    pub curves: Vec<EngineCurve>,
+    /// Isolated points.
+    pub points: Vec<V3>,
+}
+
+fn uv_list(v: Option<&Value>) -> Vec<[f64; 2]> {
+    v.and_then(Value::as_array).map(|a| a.iter().map(|p| [geom::num(&p[0]), geom::num(&p[1])]).collect()).unwrap_or_default()
+}
+
+impl EngineResult {
+    /// Reads the corpus runner's output JSON (`runner/ssi_runner.cpp`).
+    pub fn from_runner_json(o: &Value) -> Result<EngineResult, String> {
+        let curves = match o.get("curves") {
+            Some(Value::Array(a)) => a
+                .iter()
+                .map(|c| {
+                    Ok(EngineCurve {
+                        curve_type: c.get("type").and_then(Value::as_str).unwrap_or_default().to_string(),
+                        xyz: v3_list(c.get("xyz"))?,
+                        uv1: uv_list(c.get("uv1")),
+                        uv2: uv_list(c.get("uv2")),
+                        s1_at_uv: v3_list(c.get("s1_at_uv"))?,
+                        s2_at_uv: v3_list(c.get("s2_at_uv"))?,
+                    })
+                })
+                .collect::<Result<Vec<_>, String>>()?,
+            _ => return Err("curves: not a list".into()),
+        };
+        Ok(EngineResult {
+            status: o.get("status").and_then(Value::as_i64).ok_or("status: not an integer")?,
+            status_name: o.get("status_name").and_then(Value::as_str).filter(|s| !s.is_empty()).map(str::to_string),
+            time_ms: o.get("time_ms").map(geom::num).unwrap_or(f64::NAN),
+            curves,
+            points: v3_list(o.get("points"))?,
+        })
+    }
+}
+
+/// The engine ran and answered, or it did not (died, timed out, unreadable output).
+#[derive(Clone, Debug)]
+pub enum Outcome {
+    Ran(EngineResult),
+    Crashed(String),
+}
+
+#[derive(Clone, Copy, Debug, PartialEq, Eq, PartialOrd, Ord)]
+pub enum Verdict {
+    Crash,
+    Error,
+    Incomplete,
+    WrongTopology,
+    Inaccurate,
+    Pass,
+    /// Default-tolerance runs only: the case's feature size is under the judging floor.
+    NotApplicable,
+}
+
+impl Verdict {
+    pub fn as_str(self) -> &'static str {
+        match self {
+            Verdict::Crash => "CRASH",
+            Verdict::Error => "ERROR",
+            Verdict::Incomplete => "INCOMPLETE",
+            Verdict::WrongTopology => "WRONG TOPOLOGY",
+            Verdict::Inaccurate => "INACCURATE",
+            Verdict::Pass => "PASS",
+            Verdict::NotApplicable => "N/A",
+        }
+    }
+}
+
+impl std::fmt::Display for Verdict {
+    fn fmt(&self, f: &mut std::fmt::Formatter<'_>) -> std::fmt::Result {
+        f.write_str(self.as_str())
+    }
+}
+
+/// How well one expected branch / region boundary, or one isolated point, is covered.
+#[derive(Clone, Debug, PartialEq)]
+pub enum Coverage {
+    Polyline { max: f64, covered_frac: f64 },
+    Point { point: usize, dist: f64 },
+}
+
+/// The tolerances of a run.
+#[derive(Clone, Copy, Debug)]
+pub struct Tolerances {
+    /// `10 * approx_tol`; `None` for a run at the engine's default tolerance.
+    pub acc_tol: Option<f64>,
+    /// Coverage and join floor: 5e-3 for default-tolerance runs, else 0.
+    pub cov_floor: f64,
+}
+
+impl Tolerances {
+    /// A run at approximation tolerance `tol` (score.py: `acc_tol = 10 * tol`, no floor).
+    pub fn at(tol: f64) -> Tolerances {
+        Tolerances { acc_tol: Some(10.0 * tol), cov_floor: 0.0 }
+    }
+
+    /// A run at the engine's default tolerance.
+    pub fn default_tol() -> Tolerances {
+        Tolerances { acc_tol: None, cov_floor: 5e-3 }
+    }
+}
+
+/// The score of one case (score.py's result record).
+#[derive(Clone, Debug)]
+pub struct Score {
+    pub id: String,
+    pub class: String,
+    pub verdict: Verdict,
+    pub flags: Vec<String>,
+    /// The effective coverage tolerance.
+    pub cov_tol: f64,
+    /// The crash description (CRASH only).
+    pub detail: Option<String>,
+    pub status: Option<i64>,
+    pub status_name: Option<String>,
+    pub time_ms: Option<f64>,
+    pub n_curves: usize,
+    pub n_out_points: usize,
+    pub types: Vec<String>,
+    pub max_dist_3d: f64,
+    pub max_dist_uv: f64,
+    pub coverage: Vec<Coverage>,
+    pub hausdorff_truth_to_out: f64,
+    pub hausdorff_out_to_truth: f64,
+    pub expected_topology: Vec<String>,
+    pub output_topology: Vec<String>,
+}
+
+/// Python's `"%.<prec>e" % x`.
+pub fn fmt_e(x: f64, prec: usize) -> String {
+    if x.is_nan() {
+        return "nan".into();
+    }
+    if x.is_infinite() {
+        return if x > 0.0 { "inf".into() } else { "-inf".into() };
+    }
+    let s = format!("{x:.prec$e}");
+    let Some((m, e)) = s.split_once('e') else { return s };
+    let (sign, digits) = match e.strip_prefix('-') {
+        Some(d) => ('-', d),
+        None => ('+', e),
+    };
+    format!("{m}e{sign}{digits:0>2}")
+}
+
+/// Python's `"%.0f" % x`.
+fn fmt_f0(x: f64) -> String {
+    format!("{x:.0}")
+}
+
+/// Python's `str(list_of_str)`.
+fn py_list(v: &[String]) -> String {
+    format!("[{}]", v.iter().map(|s| format!("'{s}'")).collect::<Vec<_>>().join(", "))
+}
+
+fn finite(p: &[V3]) -> Vec<V3> {
+    p.iter().copied().filter(|q| q.iter().all(|c| c.is_finite())).collect()
+}
+
+/// Scores one engine outcome against a case (score.py's `score_case`).
+pub fn score_case(case: &Case, outcome: &Outcome, tol: Tolerances) -> Score {
+    let exp = &case.expected;
+    let cov_tol = case.cov_tol;
+    let mut cov_eff = cov_tol;
+    if let (Some(ill), Some(acc)) = (case.ill_cond, tol.acc_tol)
+        && ill != 0.0
+        && acc != 0.0
+    {
+        cov_eff = py_max(cov_tol, geom::py_min(0.5, 2.0 * acc * ill));
+    }
+    cov_eff = py_max(cov_eff, tol.cov_floor);
+    let join_tol = py_max(cov_tol, tol.cov_floor);
+    let mut r = Score {
+        id: case.id.clone(),
+        class: case.class.clone(),
+        verdict: Verdict::Crash,
+        flags: Vec::new(),
+        cov_tol: cov_eff,
+        detail: None,
+        status: None,
+        status_name: None,
+        time_ms: None,
+        n_curves: 0,
+        n_out_points: 0,
+        types: Vec::new(),
+        max_dist_3d: 0.0,
+        max_dist_uv: 0.0,
+        coverage: Vec::new(),
+        hausdorff_truth_to_out: 0.0,
+        hausdorff_out_to_truth: 0.0,
+        expected_topology: Vec::new(),
+        output_topology: Vec::new(),
+    };
+    let out = match outcome {
+        Outcome::Crashed(why) => {
+            r.detail = Some(why.clone());
+            return r;
+        }
+        Outcome::Ran(out) => out,
+    };
+    let status_name = out.status_name.clone().unwrap_or_else(|| out.status.to_string());
+    r.status = Some(out.status);
+    r.status_name = Some(status_name.clone());
+    r.time_ms = Some(out.time_ms);
+
+    let mut out_curves: Vec<Vec<V3>> = Vec::new();
+    let mut out_points: Vec<V3> = Vec::new();
+    let mut nonfinite = 0;
+    for c in &out.curves {
+        let p = finite(&c.xyz);
+        nonfinite += c.xyz.len() - p.len();
+        if p.is_empty() {
+            continue;
+        }
+        if c.curve_type == "TOUCHING" || diameter(&p) <= cov_tol {
+            let n = p.len() as f64;
+            let s = p.iter().fold([0.0; 3], |a, q| [a[0] + q[0], a[1] + q[1], a[2] + q[2]]);
+            out_points.push([s[0] / n, s[1] / n, s[2] / n]);
+        } else {
+            out_curves.push(p);
+        }
+        r.types.push(c.curve_type.clone());
+    }
+    out_points.extend(out.points.iter().copied());
+    r.n_curves = out.curves.len();
+    r.n_out_points = out_points.len();
+
+    // accuracy against the exact surfaces
+    let worst = |pts: &[V3]| -> f64 {
+        let mut m = np_max(&case.truth[0].distance(pts));
+        for t in &case.truth[1..] {
+            m = py_max(m, np_max(&t.distance(pts)));
+        }
+        m
+    };
+    let (mut acc3d, mut accuv) = (0.0_f64, 0.0_f64);
+    for c in &out.curves {
+        let p = finite(&c.xyz);
+        if !p.is_empty() {
+            acc3d = py_max(acc3d, worst(&p));
+        }
+        // S1(uv1) must lie on surface 2, S2(uv2) on surface 1.
+        for (k, q) in [(1, &c.s1_at_uv), (0, &c.s2_at_uv)] {
+            let q = finite(q);
+            if !q.is_empty() {
+                accuv = py_max(accuv, np_max(&case.truth[k].distance(&q)));
+            }
+        }
+    }
+    if !out_points.is_empty() {
+        acc3d = py_max(acc3d, worst(&out_points));
+    }
+    r.max_dist_3d = acc3d;
+    r.max_dist_uv = accuv;
+
+    let truth_polys: Vec<&[V3]> = exp.branches.iter().map(Vec::as_slice).chain(exp.coincident_regions.iter().map(Vec::as_slice)).collect();
+    let truth_pts = &exp.isolated_points;
+    let expected_any = !truth_polys.is_empty() || !truth_pts.is_empty();
+
+    let ok_status = out.status == SM_SUCCESS;
+    if !ok_status {
+        if !expected_any && out.curves.is_empty() && out_points.is_empty() {
+            r.flags.push(format!("status {status_name} with empty output (expected empty)"));
+        } else {
+            r.flags.push(format!("ERROR status {status_name}"));
+        }
+    }
+
+    // coverage truth -> output
+    let curve_refs: Vec<&[V3]> = out_curves.iter().map(Vec::as_slice).collect();
+    let outset = PolySet::new(&curve_refs, &out_points);
+    let mut missed = Vec::new();
+    for (k, p) in truth_polys.iter().enumerate() {
+        let d = outset.query(&stride(p, p.len() / 1000));
+        let frac = d.iter().filter(|&&x| x <= cov_eff).count() as f64 / d.len() as f64;
+        let dmax = np_max(&d);
+        r.coverage.push(Coverage::Polyline { max: dmax, covered_frac: frac });
+        if dmax > cov_eff {
+            let kind = if k >= exp.branches.len() { "region boundary".to_string() } else { format!("branch {k}") };
+            missed.push(format!("{kind} covered {}% (max gap {})", fmt_f0(100.0 * frac), fmt_e(dmax, 2)));
+        }
+    }
+    for (k, &p) in truth_pts.iter().enumerate() {
+        let d = outset.query1(p);
+        r.coverage.push(Coverage::Point { point: k, dist: d });
+        if d > cov_eff {
+            missed.push(format!("isolated point {k} not found (nearest output {})", fmt_e(d, 2)));
+        }
+    }
+    let mut h = 0.0;
+    for (i, c) in r.coverage.iter().enumerate() {
+        let x = match c {
+            Coverage::Polyline { max, .. } => *max,
+            Coverage::Point { dist, .. } => *dist,
+        };
+        h = if i == 0 { x } else { py_max(h, x) };
+    }
+    r.hausdorff_truth_to_out = h;
+
+    // output -> truth (spurious output)
+    let tset = PolySet::new(&truth_polys, truth_pts);
+    let mut spurious = Vec::new();
+    let mut h_out = 0.0;
+    for (k, p) in out_curves.iter().enumerate() {
+        let d = tset.query(&stride(p, p.len() / 1000));
+        let dmax = np_max(&d);
+        h_out = py_max(h_out, dmax);
+        if dmax > cov_eff {
+            let off = d.iter().filter(|&&x| x > cov_eff).count() as f64 / d.len() as f64;
+            spurious.push(format!("output curve {k} off the true intersection by up to {} ({}% of it)", fmt_e(dmax, 2), fmt_f0(100.0 * off)));
+        }
+    }
+    for (k, &p) in out_points.iter().enumerate() {
+        let d = tset.query1(p);
+        h_out = py_max(h_out, d);
+        if d > cov_eff {
+            spurious.push(format!("output point {k} off the true intersection by {}", fmt_e(d, 2)));
+        }
+    }
+    r.hausdorff_out_to_truth = h_out;
+
+    // topology signature
+    let truth_owned: Vec<Vec<V3>> = truth_polys.iter().map(|p| p.to_vec()).collect();
+    let sig_exp = signature(&truth_owned, truth_pts, join_tol);
+    let sig_out = signature(&out_curves, &out_points, join_tol);
+    let mut topo = Vec::new();
+    if nonfinite > 0 {
+        topo.push(format!("{nonfinite} non-finite output sample(s)"));
+    }
+    // duplicated / overlapping output curves
+    let osets: Vec<PolySet> = out_curves.iter().map(|c| PolySet::new(&[c], &[])).collect();
+    for i in 0..out_curves.len() {
+        for j in i + 1..out_curves.len() {
+            let (ci, cj) = (&out_curves[i], &out_curves[j]);
+            let di = osets[j].query(&stride(ci, ci.len() / 500));
+            let dj = osets[i].query(&stride(cj, cj.len() / 500));
+            let fi = di.iter().filter(|&&x| x <= join_tol).count() as f64 / di.len() as f64;
+            let fj = dj.iter().filter(|&&x| x <= join_tol).count() as f64 / dj.len() as f64;
+            let frac = py_max(fi, fj);
+            if frac > 0.2 {
+                topo.push(format!("output curves {i} and {j} overlap ({}% of one lies on the other)", fmt_f0(100.0 * frac)));
+            }
+        }
+    }
+    if sig_exp != sig_out {
+        topo.push(format!("topology {}, expected {}", py_list(&sig_out), py_list(&sig_exp)));
+    }
+    for (k, &s) in exp.singular_points.iter().enumerate() {
+        let ends = out_curves.iter().filter_map(|c| c.first()).chain(out_curves.iter().filter_map(|c| c.last())).chain(out_points.iter());
+        let mut d = f64::INFINITY;
+        for (n, &e) in ends.enumerate() {
+            let x = geom::norm(geom::sub(e, s));
+            d = if n == 0 { x } else { geom::py_min(d, x) };
+        }
+        if d > cov_eff {
+            topo.push(format!("singular point {k} not identified (no curve end/point within {}; nearest {})", fmt_e(cov_eff, 1), fmt_e(d, 2)));
+        }
+    }
+    let mut inacc = Vec::new();
+    if let Some(acc) = tol.acc_tol {
+        if acc3d > acc {
+            inacc.push(format!("3D output up to {} from a surface (acc_tol {})", fmt_e(acc3d, 2), fmt_e(acc, 0)));
+        }
+        if accuv > acc {
+            inacc.push(format!("surface at output UV up to {} from the other surface", fmt_e(accuv, 2)));
+        }
+    }
+
+    let has_output = !out.curves.is_empty() || !out_points.is_empty();
+    let verdict = if !ok_status && (expected_any || has_output) {
+        Verdict::Error
+    } else if !missed.is_empty() {
+        Verdict::Incomplete
+    } else if !spurious.is_empty() || !topo.is_empty() {
+        Verdict::WrongTopology
+    } else if !inacc.is_empty() {
+        Verdict::Inaccurate
+    } else {
+        Verdict::Pass
+    };
+    r.flags.extend(missed);
+    r.flags.extend(spurious);
+    r.flags.extend(topo);
+    r.flags.extend(inacc);
+    r.expected_topology = sig_exp;
+    r.output_topology = sig_out;
+    r.verdict = verdict;
+    if tol.cov_floor != 0.0
+        && let Some(fs) = case.feature_size
+        && fs < tol.cov_floor
+    {
+        r.flags.insert(0, format!("feature size {} below the {} judging floor of default-tolerance runs", fmt_e(fs, 1), fmt_e(tol.cov_floor, 0)));
+        r.verdict = Verdict::NotApplicable;
+    }
+    r
+}
+
+/// Parses `#define SM_NAME 1234L` lines of SMLib's `SmMessages.h` (first name per number), for
+/// status names in flags as score.py prints them.
+pub fn status_names(sm_messages_h: &str) -> std::collections::HashMap<i64, String> {
+    let mut names = std::collections::HashMap::new();
+    for line in sm_messages_h.lines() {
+        let Some(rest) = line.strip_prefix("#define") else { continue };
+        if !rest.starts_with(char::is_whitespace) {
+            continue;
+        }
+        let mut it = rest.split_whitespace();
+        let (Some(name), Some(val)) = (it.next(), it.next()) else { continue };
+        let Some(digits) = val.strip_suffix('L') else { continue };
+        if !name.starts_with("SM_") || !name[3..].chars().all(|c| c.is_ascii_alphanumeric() || c == '_') || name.len() == 3 {
+            continue;
+        }
+        if digits.is_empty() || !digits.chars().all(|c| c.is_ascii_digit()) {
+            continue;
+        }
+        if let Ok(n) = digits.parse::<i64>() {
+            names.entry(n).or_insert_with(|| name.to_string());
+        }
+    }
+    names
+}
+
+/// `v3` re-exported for callers building points from JSON.
+pub fn point(v: &Value) -> Result<V3, String> {
+    v3(v)
+}