fsci-conformance 0.3.0

#![forbid(unsafe_code)]
//! Differential oracle tests for FSCI-P2C-002 (Linalg).
//!
//! Runs both Rust fsci-linalg and SciPy scipy.linalg via subprocess,
//! compares outputs across 50+ inputs per function family.
//!
//! All tests emit structured JSON logs to
//! `fixtures/artifacts/FSCI-P2C-002/diff/`.

use fsci_conformance::{ArmCounts, CompareLedger};
use fsci_linalg::{
    DecompOptions, SolveOptions, det, solve, solve_with_audit, solve_with_casp, svd,
    sync_audit_ledger,
};
use fsci_runtime::{RuntimeMode, SolverPortfolio};
use serde::{Deserialize, Serialize};
use std::collections::BTreeMap;
use std::fs;
use std::path::PathBuf;
use std::time::{Instant, SystemTime, UNIX_EPOCH};

#[derive(Debug, Clone, Serialize)]
struct DiffTestLog {
    test_id: String,
    category: String,
    input_summary: String,
    expected: String,
    actual: String,
    diff: f64,
    tolerance: f64,
    pass: bool,
    compared: BTreeMap<String, ArmCounts>,
    timestamp_ms: u128,
    duration_ns: u128,
}

fn output_dir() -> PathBuf {
    PathBuf::from(env!("CARGO_MANIFEST_DIR")).join("fixtures/artifacts/FSCI-P2C-002/diff")
}

fn ensure_output_dir() {
    let dir = output_dir();
    if !dir.exists() {
        fs::create_dir_all(&dir).expect("create diff output dir");
    }
}

fn timestamp_ms() -> u128 {
    SystemTime::now()
        .duration_since(UNIX_EPOCH)
        .map_or(0, |d| d.as_millis())
}

fn emit_log(log: &DiffTestLog) {
    ensure_output_dir();
    let path = output_dir().join(format!("{}.json", log.test_id));
    let json = serde_json::to_string_pretty(log).expect("serialize log");
    fs::write(path, json).expect("write log");
}

const TOL: f64 = 1e-10;
const DEFINED_DIFF_CASES: usize = 44;
const REQUIRE_SCIPY_ENV: &str = "FSCI_REQUIRE_SCIPY_ORACLE";

fn scipy_available() -> bool {
    let output = fsci_conformance::scipy_oracle_command()
        .arg("-c")
        .arg("from scipy import linalg")
        .output();
    match output {
        Ok(o) => o.status.success(),
        Err(_) => false,
    }
}

fn scipy_available_or_skip(test_id: &str) -> bool {
    if scipy_available() {
        return true;
    }

    assert!(
        std::env::var_os(REQUIRE_SCIPY_ENV).is_none(),
        "{REQUIRE_SCIPY_ENV}=1 but SciPy is unavailable for {test_id}"
    );

    eprintln!("  SKIP: {test_id} - scipy unavailable");
    false
}

fn max_abs_diff_vec(a: &[f64], b: &[f64]) -> f64 {
    if a.len() != b.len() {
        return f64::INFINITY;
    }
    a.iter()
        .zip(b.iter())
        .map(|(x, y)| (x - y).abs())
        .fold(
            0.0_f64,
            |acc, d| if d.is_nan() { f64::NAN } else { acc.max(d) },
        )
}

#[derive(Debug, Deserialize)]
struct ScipySolveResult {
    x: Vec<f64>,
}

#[derive(Debug, Deserialize)]
struct ScipyDetResult {
    det: f64,
}

#[derive(Debug, Deserialize)]
struct ScipySvdResult {
    s: Vec<f64>,
}

fn scipy_solve(a: &[Vec<f64>], b: &[f64]) -> Option<ScipySolveResult> {
    let a_json = serde_json::to_string(a).ok()?;
    let b_json = serde_json::to_string(b).ok()?;
    let script = format!(
        r#"
import json
import numpy as np
from scipy import linalg
a = np.array({}, dtype=np.float64)
b = np.array({}, dtype=np.float64)
x = linalg.solve(a, b)
print(json.dumps({{"x": x.tolist()}}))
"#,
        a_json, b_json
    );
    let output = fsci_conformance::scipy_oracle_command()
        .arg("-c")
        .arg(&script)
        .output()
        .ok()?;
    if !output.status.success() {
        eprintln!(
            "scipy solve failed: {}",
            String::from_utf8_lossy(&output.stderr)
        );
        return None;
    }
    serde_json::from_slice(&output.stdout).ok()
}

fn scipy_det(a: &[Vec<f64>]) -> Option<ScipyDetResult> {
    let a_json = serde_json::to_string(a).ok()?;
    let script = format!(
        r#"
import json
import numpy as np
from scipy import linalg
a = np.array({}, dtype=np.float64)
d = linalg.det(a)
print(json.dumps({{"det": float(d)}}))
"#,
        a_json
    );
    let output = fsci_conformance::scipy_oracle_command()
        .arg("-c")
        .arg(&script)
        .output()
        .ok()?;
    if !output.status.success() {
        return None;
    }
    serde_json::from_slice(&output.stdout).ok()
}

fn scipy_svd(a: &[Vec<f64>]) -> Option<ScipySvdResult> {
    let a_json = serde_json::to_string(a).ok()?;
    let script = format!(
        r#"
import json
import numpy as np
from scipy import linalg
a = np.array({}, dtype=np.float64)
s = linalg.svd(a, compute_uv=False)
print(json.dumps({{"s": s.tolist()}}))
"#,
        a_json
    );
    let output = fsci_conformance::scipy_oracle_command()
        .arg("-c")
        .arg(&script)
        .output()
        .ok()?;
    if !output.status.success() {
        return None;
    }
    serde_json::from_slice(&output.stdout).ok()
}

fn make_test_matrix(n: usize, seed: u64) -> Vec<Vec<f64>> {
    let mut rng = seed;
    (0..n)
        .map(|_| {
            (0..n)
                .map(|_| {
                    rng = rng.wrapping_mul(6364136223846793005).wrapping_add(1);
                    ((rng >> 33) as f64 / (1u64 << 31) as f64) * 2.0 - 1.0
                })
                .collect()
        })
        .collect()
}

fn make_test_vector(n: usize, seed: u64) -> Vec<f64> {
    let mut rng = seed;
    (0..n)
        .map(|_| {
            rng = rng.wrapping_mul(6364136223846793005).wrapping_add(1);
            ((rng >> 33) as f64 / (1u64 << 31) as f64) * 2.0 - 1.0
        })
        .collect()
}

fn make_diag_dominant(n: usize, seed: u64) -> Vec<Vec<f64>> {
    let mut m = make_test_matrix(n, seed);
    for (i, row) in m.iter_mut().enumerate().take(n) {
        let row_sum: f64 = row.iter().map(|x| x.abs()).sum();
        row[i] = row_sum + 1.0;
    }
    m
}

fn make_hilbert_matrix(n: usize) -> Vec<Vec<f64>> {
    (0..n)
        .map(|i| (0..n).map(|j| 1.0 / ((i + j + 1) as f64)).collect())
        .collect()
}

fn run_solve_diff(test_id: &str, a: &[Vec<f64>], b: &[f64]) {
    if !scipy_available_or_skip(test_id) {
        return;
    }

    let start = Instant::now();
    let opts = SolveOptions {
        mode: RuntimeMode::Strict,
        check_finite: true,
        ..SolveOptions::default()
    };
    let rust_result = solve(a, b, opts);
    let scipy_result = scipy_solve(a, b);
    // One case per #[test]: a missing SciPy answer or a failed fsci call is recorded, and an
    // fsci error beside a missing SciPy answer no longer passes as "both raised".
    let mut ledger = CompareLedger::new(test_id, &["solve"]);
    let both_present = ledger.slices(
        "solve",
        test_id,
        scipy_result.as_ref().map(|s| s.x.as_slice()),
        rust_result.as_ref().ok().map(|r| r.x.as_slice()),
    );

    let (diff, pass, expected_str, actual_str) = match (&rust_result, &scipy_result) {
        (Ok(rust), Some(scipy)) => {
            let d = max_abs_diff_vec(&rust.x, &scipy.x);
            (
                d,
                d < TOL,
                format!("{:?}", scipy.x),
                format!("{:?}", rust.x),
            )
        }
        (Err(e), None) => (0.0, true, "error".into(), format!("{:?}", e)),
        (Ok(rust), None) => (
            f64::NAN,
            false,
            "scipy unavailable".into(),
            format!("{:?}", rust.x),
        ),
        (Err(e), Some(scipy)) => (
            f64::INFINITY,
            false,
            format!("{:?}", scipy.x),
            format!("{:?}", e),
        ),
    };
    if both_present.is_some() {
        ledger.compared("solve", test_id, pass);
    }

    let log = DiffTestLog {
        test_id: test_id.to_string(),
        category: "differential_solve".to_string(),
        input_summary: format!("{}x{} matrix", a.len(), a.first().map_or(0, |r| r.len())),
        expected: expected_str,
        actual: actual_str,
        diff,
        tolerance: TOL,
        pass,
        compared: ledger.counts().clone(),
        timestamp_ms: timestamp_ms(),
        duration_ns: start.elapsed().as_nanos(),
    };
    emit_log(&log);
    if pass {
        eprintln!("  PASS: {} — diff={:.2e}", test_id, diff);
    } else {
        eprintln!("  FAIL: {} — diff={:.2e}", test_id, diff);
    }
    assert!(pass, "Differential test {} failed: diff={}", test_id, diff);
    ledger.finish(1);
}

fn run_solve_with_casp_diff(test_id: &str, a: &[Vec<f64>], b: &[f64], tol: f64) {
    if !scipy_available_or_skip(test_id) {
        return;
    }

    let start = Instant::now();
    let opts = SolveOptions {
        mode: RuntimeMode::Strict,
        check_finite: true,
        ..SolveOptions::default()
    };
    let mut portfolio = SolverPortfolio::new(RuntimeMode::Strict, 1);
    let rust_result = solve_with_casp(a, b, opts, &mut portfolio);
    let scipy_result = scipy_solve(a, b);
    let mut ledger = CompareLedger::new(test_id, &["solve_with_casp"]);
    let both_present = ledger.slices(
        "solve_with_casp",
        test_id,
        scipy_result.as_ref().map(|s| s.x.as_slice()),
        rust_result.as_ref().ok().map(|r| r.x.as_slice()),
    );

    let (diff, pass, expected_str, actual_str) = match (&rust_result, &scipy_result) {
        (Ok(rust), Some(scipy)) => {
            let d = max_abs_diff_vec(&rust.x, &scipy.x);
            let has_certificate = rust.certificate.is_some();
            (
                d,
                d < tol && has_certificate,
                format!("{:?}", scipy.x),
                format!(
                    "{:?} (cert={}, action={:?}, fallback={:?})",
                    rust.x,
                    has_certificate,
                    rust.certificate.as_ref().map(|c| c.action),
                    rust.certificate.as_ref().map(|c| c.fallback_active)
                ),
            )
        }
        (Err(e), None) => (0.0, true, "error".into(), format!("{:?}", e)),
        (Ok(rust), None) => (
            f64::NAN,
            false,
            "scipy unavailable".into(),
            format!("{:?}", rust.x),
        ),
        (Err(e), Some(scipy)) => (
            f64::INFINITY,
            false,
            format!("{:?}", scipy.x),
            format!("{:?}", e),
        ),
    };
    if both_present.is_some() {
        ledger.compared("solve_with_casp", test_id, pass);
    }

    let log = DiffTestLog {
        test_id: test_id.to_string(),
        category: "differential_solve_with_casp".to_string(),
        input_summary: format!("{}x{} matrix", a.len(), a.first().map_or(0, |r| r.len())),
        expected: expected_str,
        actual: actual_str,
        diff,
        tolerance: tol,
        pass,
        compared: ledger.counts().clone(),
        timestamp_ms: timestamp_ms(),
        duration_ns: start.elapsed().as_nanos(),
    };
    emit_log(&log);
    if pass {
        eprintln!("  PASS: {} — diff={:.2e}", test_id, diff);
    } else {
        eprintln!("  FAIL: {} — diff={:.2e}", test_id, diff);
    }
    // Both vectors in the failure: on an ill-conditioned matrix SciPy's own answer moves with
    // the OpenBLAS kernel (named on the `scipy_incumbent` line), so a CI failure must show
    // which side moved.
    assert!(
        pass,
        "Differential test {} failed: diff={}; scipy={}; fsci={}",
        test_id, diff, log.expected, log.actual
    );
    ledger.finish(1);
}

fn run_solve_with_audit_diff(test_id: &str, a: &[Vec<f64>], b: &[f64], tol: f64) {
    if !scipy_available_or_skip(test_id) {
        return;
    }

    let start = Instant::now();
    let opts = SolveOptions {
        mode: RuntimeMode::Strict,
        check_finite: true,
        ..SolveOptions::default()
    };
    let mut portfolio = SolverPortfolio::new(RuntimeMode::Strict, 1);
    let ledger = sync_audit_ledger();
    let initial_entries = ledger.lock().map_or(0, |g| g.len());
    let rust_result = solve_with_audit(a, b, opts, &mut portfolio, &ledger);
    let scipy_result = scipy_solve(a, b);
    let post_entries = ledger.lock().map_or(0, |g| g.len());
    let recorded_audit = post_entries > initial_entries;
    // `ledger` above is fsci's audit ledger; this one records the SciPy comparison.
    let mut compare_ledger = CompareLedger::new(test_id, &["solve_with_audit"]);
    let both_present = compare_ledger.slices(
        "solve_with_audit",
        test_id,
        scipy_result.as_ref().map(|s| s.x.as_slice()),
        rust_result.as_ref().ok().map(|r| r.x.as_slice()),
    );

    let (diff, pass, expected_str, actual_str) = match (&rust_result, &scipy_result) {
        (Ok(rust), Some(scipy)) => {
            let d = max_abs_diff_vec(&rust.x, &scipy.x);
            let has_certificate = rust.certificate.is_some();
            (
                d,
                d < tol && has_certificate && recorded_audit,
                format!("{:?}", scipy.x),
                format!(
                    "{:?} (cert={}, action={:?}, fallback={:?}, audit_recorded={})",
                    rust.x,
                    has_certificate,
                    rust.certificate.as_ref().map(|c| c.action),
                    rust.certificate.as_ref().map(|c| c.fallback_active),
                    recorded_audit
                ),
            )
        }
        (Err(e), None) => (
            0.0,
            recorded_audit,
            "error".into(),
            format!("{:?} (audit_recorded={})", e, recorded_audit),
        ),
        (Ok(rust), None) => (
            f64::NAN,
            false,
            "scipy unavailable".into(),
            format!("{:?}", rust.x),
        ),
        (Err(e), Some(scipy)) => (
            f64::INFINITY,
            false,
            format!("{:?}", scipy.x),
            format!("{:?}", e),
        ),
    };
    if both_present.is_some() {
        compare_ledger.compared("solve_with_audit", test_id, pass);
    }

    let log = DiffTestLog {
        test_id: test_id.to_string(),
        category: "differential_solve_with_audit".to_string(),
        input_summary: format!("{}x{} matrix", a.len(), a.first().map_or(0, |r| r.len())),
        expected: expected_str,
        actual: actual_str,
        diff,
        tolerance: tol,
        pass,
        compared: compare_ledger.counts().clone(),
        timestamp_ms: timestamp_ms(),
        duration_ns: start.elapsed().as_nanos(),
    };
    emit_log(&log);
    if pass {
        eprintln!("  PASS: {} — diff={:.2e}", test_id, diff);
    } else {
        eprintln!("  FAIL: {} — diff={:.2e}", test_id, diff);
    }
    assert!(
        pass,
        "Differential test {} failed: diff={}, audit_recorded={}; scipy={}; fsci={}",
        test_id, diff, recorded_audit, log.expected, log.actual
    );
    compare_ledger.finish(1);
}

fn run_det_diff(test_id: &str, a: &[Vec<f64>]) {
    if !scipy_available_or_skip(test_id) {
        return;
    }

    let start = Instant::now();
    let rust_result = det(a, RuntimeMode::Strict, true);
    let scipy_result = scipy_det(a);
    let mut ledger = CompareLedger::new(test_id, &["det"]);
    let both_present = ledger.pair(
        "det",
        test_id,
        scipy_result.as_ref().map(|s| s.det),
        rust_result.as_ref().ok().copied(),
    );

    let (diff, pass, expected_str, actual_str) = match (&rust_result, &scipy_result) {
        (Ok(rust_det), Some(scipy)) => {
            let d = (rust_det - scipy.det).abs();
            let rel_tol = TOL * rust_det.abs().max(scipy.det.abs()).max(1.0);
            (
                d,
                d < rel_tol,
                format!("{}", scipy.det),
                format!("{}", rust_det),
            )
        }
        (Err(e), None) => (0.0, true, "error".into(), format!("{:?}", e)),
        (Ok(rust_det), None) => (
            f64::NAN,
            false,
            "scipy unavailable".into(),
            format!("{}", rust_det),
        ),
        (Err(e), Some(scipy)) => (
            f64::INFINITY,
            false,
            format!("{}", scipy.det),
            format!("{:?}", e),
        ),
    };
    if both_present.is_some() {
        ledger.compared("det", test_id, pass);
    }

    let log = DiffTestLog {
        test_id: test_id.to_string(),
        category: "differential_det".to_string(),
        input_summary: format!("{}x{} matrix", a.len(), a.first().map_or(0, |r| r.len())),
        expected: expected_str,
        actual: actual_str,
        diff,
        tolerance: TOL,
        pass,
        compared: ledger.counts().clone(),
        timestamp_ms: timestamp_ms(),
        duration_ns: start.elapsed().as_nanos(),
    };
    emit_log(&log);
    if pass {
        eprintln!("  PASS: {} — diff={:.2e}", test_id, diff);
    } else {
        eprintln!("  FAIL: {} — diff={:.2e}", test_id, diff);
    }
    assert!(pass, "Differential test {} failed: diff={}", test_id, diff);
    ledger.finish(1);
}

fn run_svd_singular_values_diff(test_id: &str, a: &[Vec<f64>]) {
    if !scipy_available_or_skip(test_id) {
        return;
    }

    let start = Instant::now();
    let opts = DecompOptions {
        mode: RuntimeMode::Strict,
        check_finite: true,
    };
    let rust_result = svd(a, opts);
    let scipy_result = scipy_svd(a);
    let mut ledger = CompareLedger::new(test_id, &["svd"]);
    let both_present = ledger.slices(
        "svd",
        test_id,
        scipy_result.as_ref().map(|s| s.s.as_slice()),
        rust_result.as_ref().ok().map(|r| r.s.as_slice()),
    );

    let (diff, pass, expected_str, actual_str) = match (&rust_result, &scipy_result) {
        (Ok(rust), Some(scipy)) => {
            let d = max_abs_diff_vec(&rust.s, &scipy.s);
            (
                d,
                d < TOL,
                format!("{:?}", scipy.s),
                format!("{:?}", rust.s),
            )
        }
        (Err(e), None) => (0.0, true, "error".into(), format!("{:?}", e)),
        (Ok(rust), None) => (
            f64::NAN,
            false,
            "scipy unavailable".into(),
            format!("{:?}", rust.s),
        ),
        (Err(e), Some(scipy)) => (
            f64::INFINITY,
            false,
            format!("{:?}", scipy.s),
            format!("{:?}", e),
        ),
    };
    if both_present.is_some() {
        ledger.compared("svd", test_id, pass);
    }

    let log = DiffTestLog {
        test_id: test_id.to_string(),
        category: "differential_svd".to_string(),
        input_summary: format!("{}x{} matrix", a.len(), a.first().map_or(0, |r| r.len())),
        expected: expected_str,
        actual: actual_str,
        diff,
        tolerance: TOL,
        pass,
        compared: ledger.counts().clone(),
        timestamp_ms: timestamp_ms(),
        duration_ns: start.elapsed().as_nanos(),
    };
    emit_log(&log);
    if pass {
        eprintln!("  PASS: {} — diff={:.2e}", test_id, diff);
    } else {
        eprintln!("  FAIL: {} — diff={:.2e}", test_id, diff);
    }
    assert!(pass, "Differential test {} failed: diff={}", test_id, diff);
    ledger.finish(1);
}

// ═══════════════════════════════════════════════════════════════
// DIFFERENTIAL TESTS: solve (20 cases)
// ═══════════════════════════════════════════════════════════════

#[test]
fn diff_solve_2x2_identity() {
    let a = vec![vec![1.0, 0.0], vec![0.0, 1.0]];
    let b = vec![1.0, 2.0];
    run_solve_diff("solve_2x2_identity", &a, &b);
}

#[test]
fn diff_solve_3x3_diag_dominant() {
    let a = make_diag_dominant(3, 0xCAFE);
    let b = make_test_vector(3, 0xBEEF);
    run_solve_diff("solve_3x3_diag_dominant", &a, &b);
}

#[test]
fn diff_solve_4x4_random() {
    let a = make_diag_dominant(4, 0x1234);
    let b = make_test_vector(4, 0x5678);
    run_solve_diff("solve_4x4_random", &a, &b);
}

#[test]
fn diff_solve_5x5_random() {
    let a = make_diag_dominant(5, 0xABCD);
    let b = make_test_vector(5, 0xEF01);
    run_solve_diff("solve_5x5_random", &a, &b);
}

#[test]
fn diff_solve_8x8_random() {
    let a = make_diag_dominant(8, 0x2345);
    let b = make_test_vector(8, 0x6789);
    run_solve_diff("solve_8x8_random", &a, &b);
}

#[test]
fn diff_solve_10x10_random() {
    let a = make_diag_dominant(10, 0x3456);
    let b = make_test_vector(10, 0x789A);
    run_solve_diff("solve_10x10_random", &a, &b);
}

#[test]
fn diff_solve_16x16_random() {
    let a = make_diag_dominant(16, 0x4567);
    let b = make_test_vector(16, 0x89AB);
    run_solve_diff("solve_16x16_random", &a, &b);
}

#[test]
fn diff_solve_20x20_random() {
    let a = make_diag_dominant(20, 0x5678);
    let b = make_test_vector(20, 0x9ABC);
    run_solve_diff("solve_20x20_random", &a, &b);
}

#[test]
fn diff_solve_32x32_random() {
    let a = make_diag_dominant(32, 0x6789);
    let b = make_test_vector(32, 0xABCD);
    run_solve_diff("solve_32x32_random", &a, &b);
}

#[test]
fn diff_solve_50x50_random() {
    let a = make_diag_dominant(50, 0x789A);
    let b = make_test_vector(50, 0xBCDE);
    run_solve_diff("solve_50x50_random", &a, &b);
}

#[test]
fn diff_solve_tridiag_5x5() {
    let a = vec![
        vec![2.0, -1.0, 0.0, 0.0, 0.0],
        vec![-1.0, 2.0, -1.0, 0.0, 0.0],
        vec![0.0, -1.0, 2.0, -1.0, 0.0],
        vec![0.0, 0.0, -1.0, 2.0, -1.0],
        vec![0.0, 0.0, 0.0, -1.0, 2.0],
    ];
    let b = vec![1.0, 0.0, 0.0, 0.0, 1.0];
    run_solve_diff("solve_tridiag_5x5", &a, &b);
}

#[test]
fn diff_solve_hilbert_like_4x4() {
    let a = vec![
        vec![1.0, 0.5, 0.333, 0.25],
        vec![0.5, 0.333, 0.25, 0.2],
        vec![0.333, 0.25, 0.2, 0.167],
        vec![0.25, 0.2, 0.167, 0.143],
    ];
    let b = vec![1.0, 1.0, 1.0, 1.0];
    run_solve_diff("solve_hilbert_like_4x4", &a, &b);
}

// ═══════════════════════════════════════════════════════════════
// DIFFERENTIAL TESTS: solve_with_casp (5 cases)
// ═══════════════════════════════════════════════════════════════

#[test]
fn diff_solve_with_casp_hilbert_4x4() {
    let a = make_hilbert_matrix(4);
    let b = vec![1.0, 1.0, 1.0, 1.0];
    run_solve_with_casp_diff("solve_with_casp_hilbert_4x4", &a, &b, 1e-8);
}

#[test]
fn diff_solve_with_casp_hilbert_6x6() {
    let a = make_hilbert_matrix(6);
    let b = vec![1.0, 1.0, 1.0, 1.0, 1.0, 1.0];
    // Hilbert 6x6 has cond ~1.5e7 with solution components up to 6300;
    // 5e-7 abs tol corresponds to ~2.4e-11 relative error.
    run_solve_with_casp_diff("solve_with_casp_hilbert_6x6", &a, &b, 5e-7);
}

/// `Pᵀ·diag(-1, 1, 1, -1)·P` with `P` unit upper triangular and integer: symmetric
/// indefinite, det = 1, cond ≈ 1.3e9 (frankenscipy-7tb8d.15). SciPy factors it with
/// Bunch–Kaufman LDLᵀ (sytrf): its answer (|x| up to 5.0e6) moves by at most 1.2e-10 across
/// OpenBLAS's kernels, while an LU answer lands 2.4e-3 away.
#[test]
fn diff_solve_with_casp_symmetric_indefinite_4x4() {
    let a = vec![
        vec![-1.0, -9.0, 2.0, -5.0],
        vec![-9.0, -80.0, 30.0, -48.0],
        vec![2.0, 30.0, 141.0, -8.0],
        vec![-5.0, -48.0, -8.0, 307.0],
    ];
    let b = vec![1.0, -2.0, 3.0, -4.0];
    run_solve_with_casp_diff("solve_with_casp_symmetric_indefinite_4x4", &a, &b, 1e-8);
}

#[test]
fn diff_solve_with_casp_diag_dominant_4x4() {
    let a = make_diag_dominant(4, 0x1234);
    let b = make_test_vector(4, 0x5678);
    run_solve_with_casp_diff("solve_with_casp_diag_dominant_4x4", &a, &b, 1e-10);
}

#[test]
fn diff_solve_with_casp_diag_dominant_8x8() {
    let a = make_diag_dominant(8, 0x2345);
    let b = make_test_vector(8, 0x6789);
    run_solve_with_casp_diff("solve_with_casp_diag_dominant_8x8", &a, &b, 1e-10);
}

// ═══════════════════════════════════════════════════════════════
// DIFFERENTIAL TESTS: solve_with_audit (4 cases)
// ═══════════════════════════════════════════════════════════════

#[test]
fn diff_solve_with_audit_hilbert_4x4() {
    let a = make_hilbert_matrix(4);
    let b = vec![1.0, 1.0, 1.0, 1.0];
    run_solve_with_audit_diff("solve_with_audit_hilbert_4x4", &a, &b, 1e-8);
}

#[test]
fn diff_solve_with_audit_hilbert_6x6() {
    let a = make_hilbert_matrix(6);
    let b = vec![1.0, 1.0, 1.0, 1.0, 1.0, 1.0];
    // Hilbert 6x6 has cond ~1.5e7 with solution components up to 6300;
    // 5e-7 abs tol corresponds to ~2.4e-11 relative error.
    run_solve_with_audit_diff("solve_with_audit_hilbert_6x6", &a, &b, 5e-7);
}

#[test]
fn diff_solve_with_audit_diag_dominant_4x4() {
    let a = make_diag_dominant(4, 0x1234);
    let b = make_test_vector(4, 0x5678);
    run_solve_with_audit_diff("solve_with_audit_diag_dominant_4x4", &a, &b, 1e-10);
}

#[test]
fn diff_solve_with_audit_diag_dominant_8x8() {
    let a = make_diag_dominant(8, 0x2345);
    let b = make_test_vector(8, 0x6789);
    run_solve_with_audit_diff("solve_with_audit_diag_dominant_8x8", &a, &b, 1e-10);
}

// ═══════════════════════════════════════════════════════════════
// DIFFERENTIAL TESTS: det (15 cases)
// ═══════════════════════════════════════════════════════════════

#[test]
fn diff_det_2x2_identity() {
    let a = vec![vec![1.0, 0.0], vec![0.0, 1.0]];
    run_det_diff("det_2x2_identity", &a);
}

#[test]
fn diff_det_2x2_simple() {
    let a = vec![vec![3.0, 1.0], vec![2.0, 4.0]];
    run_det_diff("det_2x2_simple", &a);
}

#[test]
fn diff_det_3x3_diag() {
    let a = vec![
        vec![2.0, 0.0, 0.0],
        vec![0.0, 3.0, 0.0],
        vec![0.0, 0.0, 4.0],
    ];
    run_det_diff("det_3x3_diag", &a);
}

#[test]
fn diff_det_3x3_random() {
    let a = make_test_matrix(3, 0xDEAD);
    run_det_diff("det_3x3_random", &a);
}

#[test]
fn diff_det_4x4_random() {
    let a = make_test_matrix(4, 0xBEEF);
    run_det_diff("det_4x4_random", &a);
}

#[test]
fn diff_det_5x5_random() {
    let a = make_test_matrix(5, 0xCAFE);
    run_det_diff("det_5x5_random", &a);
}

#[test]
fn diff_det_8x8_random() {
    let a = make_test_matrix(8, 0x1111);
    run_det_diff("det_8x8_random", &a);
}

#[test]
fn diff_det_10x10_random() {
    let a = make_test_matrix(10, 0x2222);
    run_det_diff("det_10x10_random", &a);
}

#[test]
fn diff_det_16x16_random() {
    let a = make_test_matrix(16, 0x3333);
    run_det_diff("det_16x16_random", &a);
}

#[test]
fn diff_det_20x20_random() {
    let a = make_test_matrix(20, 0x4444);
    run_det_diff("det_20x20_random", &a);
}

// ═══════════════════════════════════════════════════════════════
// DIFFERENTIAL TESTS: SVD singular values (20 cases)
// ═══════════════════════════════════════════════════════════════

#[test]
fn diff_svd_2x2_identity() {
    let a = vec![vec![1.0, 0.0], vec![0.0, 1.0]];
    run_svd_singular_values_diff("svd_2x2_identity", &a);
}

#[test]
fn diff_svd_2x2_diag() {
    let a = vec![vec![3.0, 0.0], vec![0.0, 5.0]];
    run_svd_singular_values_diff("svd_2x2_diag", &a);
}

#[test]
fn diff_svd_3x3_random() {
    let a = make_test_matrix(3, 0xAAAA);
    run_svd_singular_values_diff("svd_3x3_random", &a);
}

#[test]
fn diff_svd_4x4_random() {
    let a = make_test_matrix(4, 0xBBBB);
    run_svd_singular_values_diff("svd_4x4_random", &a);
}

#[test]
fn diff_svd_5x5_random() {
    let a = make_test_matrix(5, 0xCCCC);
    run_svd_singular_values_diff("svd_5x5_random", &a);
}

#[test]
fn diff_svd_8x8_random() {
    let a = make_test_matrix(8, 0xDDDD);
    run_svd_singular_values_diff("svd_8x8_random", &a);
}

#[test]
fn diff_svd_10x10_random() {
    let a = make_test_matrix(10, 0xEEEE);
    run_svd_singular_values_diff("svd_10x10_random", &a);
}

#[test]
fn diff_svd_16x16_random() {
    let a = make_test_matrix(16, 0xFFFF);
    run_svd_singular_values_diff("svd_16x16_random", &a);
}

#[test]
fn diff_svd_20x20_random() {
    let a = make_test_matrix(20, 0x1234);
    run_svd_singular_values_diff("svd_20x20_random", &a);
}

#[test]
fn diff_svd_32x32_random() {
    let a = make_test_matrix(32, 0x5678);
    run_svd_singular_values_diff("svd_32x32_random", &a);
}

#[test]
fn diff_svd_2x3_rect() {
    let a = vec![vec![1.0, 2.0, 3.0], vec![4.0, 5.0, 6.0]];
    run_svd_singular_values_diff("svd_2x3_rect", &a);
}

#[test]
fn diff_svd_3x2_rect() {
    let a = vec![vec![1.0, 2.0], vec![3.0, 4.0], vec![5.0, 6.0]];
    run_svd_singular_values_diff("svd_3x2_rect", &a);
}

#[test]
fn diff_svd_4x6_rect() {
    let a = vec![
        vec![1.0, 2.0, 3.0, 4.0, 5.0, 6.0],
        vec![2.0, 3.0, 4.0, 5.0, 6.0, 7.0],
        vec![3.0, 4.0, 5.0, 6.0, 7.0, 8.0],
        vec![4.0, 5.0, 6.0, 7.0, 8.0, 9.0],
    ];
    run_svd_singular_values_diff("svd_4x6_rect", &a);
}

#[test]
fn diff_svd_6x4_rect() {
    let a = vec![
        vec![1.0, 2.0, 3.0, 4.0],
        vec![2.0, 3.0, 4.0, 5.0],
        vec![3.0, 4.0, 5.0, 6.0],
        vec![4.0, 5.0, 6.0, 7.0],
        vec![5.0, 6.0, 7.0, 8.0],
        vec![6.0, 7.0, 8.0, 9.0],
    ];
    run_svd_singular_values_diff("svd_6x4_rect", &a);
}

// ═══════════════════════════════════════════════════════════════
// SUMMARY
// ═══════════════════════════════════════════════════════════════

#[test]
fn diff_linalg_summary() {
    eprintln!("\n── diff_linalg Summary ──");
    if scipy_available() {
        eprintln!("  solve differential tests: 12 cases");
        eprintln!("  solve_with_casp differential tests: 4 cases");
        eprintln!("  solve_with_audit differential tests: 4 cases");
        eprintln!("  det differential tests: 10 cases");
        eprintln!("  svd differential tests: 14 cases");
        eprintln!("  Total: {DEFINED_DIFF_CASES} subprocess-based differential tests");
        eprintln!("  Logs: fixtures/artifacts/FSCI-P2C-002/diff/");
    } else {
        eprintln!("  SciPy unavailable: {DEFINED_DIFF_CASES} oracle comparisons skipped");
        eprintln!("  Set {REQUIRE_SCIPY_ENV}=1 to fail closed when the oracle is missing");
    }
}