g_math 0.4.2

//! ULP (Unit in Last Place) measurement for decomposition outputs.
//!
//! Measures the ACTUAL ULP error of every decomposition operation against
//! mpmath 80-digit reference values. This establishes the empirical precision
//! guarantee for the storage tier.
//!
//! ULP is measured as: |gmath_raw - reference_raw| where both are in Q-format
//! raw storage representation. 1 ULP = the smallest representable difference.

use g_math::fixed_point::{FixedPoint, FixedVector, FixedMatrix};
use g_math::fixed_point::imperative::decompose::{
    lu_decompose, qr_decompose, cholesky_decompose,
};

fn fp(s: &str) -> FixedPoint {
    if s.starts_with('-') { -FixedPoint::from_str(&s[1..]) }
    else { FixedPoint::from_str(s) }
}

/// Measure ULP difference between two FixedPoint values.
/// Returns the absolute raw storage difference (number of quanta apart).
fn ulp_diff(a: FixedPoint, b: FixedPoint) -> u64 {
    let diff = a - b;
    let raw = diff.raw();

    #[cfg(any(table_format = "q16_16", table_format = "q32_32", table_format = "q64_64"))]
    { raw.unsigned_abs() as u64 }

    #[cfg(table_format = "q128_128")]
    {
        // I256 → take lower 64 bits of absolute value for ULP count
        let abs_raw = if raw.is_negative() { -raw } else { raw };
        abs_raw.as_i128() as u64
    }

    #[cfg(table_format = "q256_256")]
    {
        let abs_raw = if raw.is_negative() { -raw } else { raw };
        abs_raw.as_i128() as u64
    }
}

fn report_ulp(name: &str, got: FixedPoint, reference: FixedPoint) -> u64 {
    let ulp = ulp_diff(got, reference);
    println!("  {:<40} {:>6} ULP  (got={}, ref={})", name, ulp, got, reference);
    ulp
}

fn report_vec_ulp(name: &str, got: &FixedVector, refs: &[&str]) -> u64 {
    let mut max_ulp = 0u64;
    for i in 0..got.len() {
        let ulp = report_ulp(&format!("{}[{}]", name, i), got[i], fp(refs[i]));
        max_ulp = max_ulp.max(ulp);
    }
    max_ulp
}

// ============================================================================
// ULP measurements
// ============================================================================

#[test]
fn test_ulp_measurement_report() {
    println!("\n╔══════════════════════════════════════════════════════════════╗");
    println!("║        DECOMPOSITION ULP MEASUREMENT REPORT                ║");
    println!("╚══════════════════════════════════════════════════════════════╝\n");

    let mut max_ulp_overall = 0u64;

    // ── 1. LU Solve: well-conditioned 3x3 (exact rational solution) ──
    println!("── LU Solve: 3x3 well-conditioned (exact rational answers) ──");
    {
        let a = FixedMatrix::from_slice(3, 3, &[
            fp("3"), fp("1"), fp("2"),
            fp("1"), fp("4"), fp("1"),
            fp("2"), fp("1"), fp("5"),
        ]);
        let b = FixedVector::from_slice(&[fp("1"), fp("2"), fp("3")]);
        let lu = lu_decompose(&a).unwrap();
        let x = lu.solve(&b).unwrap();
        // mpmath: x = [-0.2, 0.4, 0.6]
        let ulp = report_vec_ulp("LU solve 3x3", &x, &["-0.2", "0.4", "0.6"]);
        max_ulp_overall = max_ulp_overall.max(ulp);
    }
    println!();

    // ── 2. QR Solve: same system ──
    println!("── QR Solve: 3x3 well-conditioned ──");
    {
        let a = FixedMatrix::from_slice(3, 3, &[
            fp("3"), fp("1"), fp("2"),
            fp("1"), fp("4"), fp("1"),
            fp("2"), fp("1"), fp("5"),
        ]);
        let b = FixedVector::from_slice(&[fp("1"), fp("2"), fp("3")]);
        let qr = qr_decompose(&a).unwrap();
        let x = qr.solve(&b).unwrap();
        let ulp = report_vec_ulp("QR solve 3x3", &x, &["-0.2", "0.4", "0.6"]);
        max_ulp_overall = max_ulp_overall.max(ulp);
    }
    println!();

    // ── 3. Cholesky Solve: same system (it's SPD) ──
    println!("── Cholesky Solve: 3x3 SPD ──");
    {
        let a = FixedMatrix::from_slice(3, 3, &[
            fp("3"), fp("1"), fp("2"),
            fp("1"), fp("4"), fp("1"),
            fp("2"), fp("1"), fp("5"),
        ]);
        let b = FixedVector::from_slice(&[fp("1"), fp("2"), fp("3")]);
        let chol = cholesky_decompose(&a).unwrap();
        let x = chol.solve(&b).unwrap();
        let ulp = report_vec_ulp("Cholesky solve 3x3", &x, &["-0.2", "0.4", "0.6"]);
        max_ulp_overall = max_ulp_overall.max(ulp);
    }
    println!();

    // ── 4. LU Determinant ──
    println!("── Determinants ──");
    {
        let a = FixedMatrix::from_slice(3, 3, &[
            fp("3"), fp("1"), fp("2"),
            fp("1"), fp("4"), fp("1"),
            fp("2"), fp("1"), fp("5"),
        ]);
        let lu = lu_decompose(&a).unwrap();
        let ulp = report_ulp("det(3x3 SPD) = 40", lu.determinant(), fp("40"));
        max_ulp_overall = max_ulp_overall.max(ulp);

        let a2 = FixedMatrix::from_slice(2, 2, &[fp("1"), fp("2"), fp("3"), fp("4")]);
        let lu2 = lu_decompose(&a2).unwrap();
        let ulp = report_ulp("det(2x2) = -2", lu2.determinant(), fp("-2"));
        max_ulp_overall = max_ulp_overall.max(ulp);
    }
    println!();

    // ── 5. Cholesky factors with irrational values ──
    println!("── Cholesky factors (irrational entries — sqrt involved) ──");
    {
        // [[4,2],[2,3]] → L = [[2,0],[1,sqrt(2)]]
        let a = FixedMatrix::from_slice(2, 2, &[fp("4"), fp("2"), fp("2"), fp("3")]);
        let chol = cholesky_decompose(&a).unwrap();
        let ulp = report_ulp("L[0][0] = 2", chol.l.get(0, 0), fp("2"));
        max_ulp_overall = max_ulp_overall.max(ulp);
        let ulp = report_ulp("L[1][0] = 1", chol.l.get(1, 0), fp("1"));
        max_ulp_overall = max_ulp_overall.max(ulp);
        // L[1][1] = sqrt(2) — mpmath: 1.41421356237309504880168872420969807856967...
        #[cfg(any(table_format = "q16_16", table_format = "q32_32", table_format = "q64_64"))]
        let sqrt2_first = "1.4142135623730950488";
        #[cfg(table_format = "q128_128")]
        let sqrt2_first = "1.41421356237309504880168872420969807856";
        #[cfg(table_format = "q256_256")]
        let sqrt2_first = "1.4142135623730950488016887242096980785696718753769480731766797379907324784621";
        let ulp = report_ulp("L[1][1] = sqrt(2)",
            chol.l.get(1, 1), fp(sqrt2_first));
        max_ulp_overall = max_ulp_overall.max(ulp);
    }
    println!();

    {
        // [[2,1],[1,2]] → L[0][0] = sqrt(2), L[1][0] = 1/sqrt(2), L[1][1] = sqrt(3/2)
        let a = FixedMatrix::from_slice(2, 2, &[fp("2"), fp("1"), fp("1"), fp("2")]);
        let chol = cholesky_decompose(&a).unwrap();
        // mpmath 50-digit references for Cholesky factors
        // sqrt(2)   = 1.41421356237309504880168872420969807856967187537694...
        // 1/sqrt(2)  = 0.70710678118654752440084436210484903928483593768474...
        // sqrt(3/2)  = 1.22474487139158904909864203735294569598297374032833...
        #[cfg(any(table_format = "q16_16", table_format = "q32_32", table_format = "q64_64"))]
        let (sqrt2, inv_sqrt2, sqrt_3_2) = (
            "1.4142135623730950488",
            "0.70710678118654752440",
            "1.2247448713915890490",
        );
        #[cfg(table_format = "q128_128")]
        let (sqrt2, inv_sqrt2, sqrt_3_2) = (
            "1.41421356237309504880168872420969807856",
            "0.70710678118654752440084436210484903928",
            "1.22474487139158904909864203735294569598",
        );
        #[cfg(table_format = "q256_256")]
        let (sqrt2, inv_sqrt2, sqrt_3_2) = (
            "1.4142135623730950488016887242096980785696718753769480731766797379907324784621",
            "0.7071067811865475244008443621048490392848359376884740365883398689953662392310",
            "1.2247448713915890490986420373529456959829737403283350642163462836254801887286",
        );
        let ulp = report_ulp("L[0][0] = sqrt(2)",
            chol.l.get(0, 0), fp(sqrt2));
        max_ulp_overall = max_ulp_overall.max(ulp);
        let ulp = report_ulp("L[1][0] = 1/sqrt(2)",
            chol.l.get(1, 0), fp(inv_sqrt2));
        max_ulp_overall = max_ulp_overall.max(ulp);
        let ulp = report_ulp("L[1][1] = sqrt(3/2)",
            chol.l.get(1, 1), fp(sqrt_3_2));
        max_ulp_overall = max_ulp_overall.max(ulp);
    }
    println!();

    // ── 6. QR orthogonality: measure Q^T Q - I ──
    println!("── QR Orthogonality: max |Q^T Q - I| in ULP ──");
    {
        let a = FixedMatrix::from_slice(3, 3, &[
            fp("12"), fp("-51"), fp("4"),
            fp("6"), fp("167"), fp("-68"),
            fp("-4"), fp("24"), fp("-41"),
        ]);
        let qr = qr_decompose(&a).unwrap();
        let qtq = &qr.q.transpose() * &qr.q;
        let mut max_orth_ulp = 0u64;
        for i in 0..3 {
            for j in 0..3 {
                let expected = if i == j { fp("1") } else { fp("0") };
                let ulp = ulp_diff(qtq.get(i, j), expected);
                max_orth_ulp = max_orth_ulp.max(ulp);
            }
        }
        println!("  {:<40} {:>6} ULP", "max |Q^TQ - I| element", max_orth_ulp);
        max_ulp_overall = max_ulp_overall.max(max_orth_ulp);
    }
    println!();

    // ── 7. Inverse: A * A^-1 - I ──
    println!("── Inverse accuracy: max |A * A^-1 - I| in ULP ──");
    {
        let a = FixedMatrix::from_slice(3, 3, &[
            fp("1"), fp("2"), fp("3"),
            fp("0"), fp("1"), fp("4"),
            fp("5"), fp("6"), fp("0"),
        ]);
        let lu = lu_decompose(&a).unwrap();
        let inv = lu.inverse().unwrap();
        let product = &a * &inv;
        let mut max_inv_ulp = 0u64;
        for i in 0..3 {
            for j in 0..3 {
                let expected = if i == j { fp("1") } else { fp("0") };
                let ulp = ulp_diff(product.get(i, j), expected);
                max_inv_ulp = max_inv_ulp.max(ulp);
            }
        }
        println!("  {:<40} {:>6} ULP", "max |A*A^-1 - I| element", max_inv_ulp);
        max_ulp_overall = max_ulp_overall.max(max_inv_ulp);
    }
    println!();

    // ── 8. Hilbert 4x4 (ill-conditioned, cond ≈ 28375) ──
    println!("── Hilbert 4x4 solve (cond ≈ 28375) ──");
    {
        let h = FixedMatrix::from_fn(4, 4, |i, j| {
            FixedPoint::one() / FixedPoint::from_int((i + j + 1) as i32)
        });
        let b = FixedVector::from_slice(&[fp("1"), fp("1"), fp("1"), fp("1")]);
        let lu = lu_decompose(&h).unwrap();
        let x = lu.solve(&b).unwrap();
        let _ulp = report_vec_ulp("Hilbert LU solve (raw)", &x, &["-4", "60", "-180", "140"]);
        // Don't include raw Hilbert in max — it's condition-number dominated

        // Iterative refinement: one step should recover nearly all precision
        let x_refined = lu.refine(&h, &b, &x).unwrap();
        let ulp_refined = report_vec_ulp("Hilbert REFINED (1 step)", &x_refined, &["-4", "60", "-180", "140"]);
        max_ulp_overall = max_ulp_overall.max(ulp_refined);

        // Second refinement step
        let x_refined2 = lu.refine(&h, &b, &x_refined).unwrap();
        let ulp_refined2 = report_vec_ulp("Hilbert REFINED (2 steps)", &x_refined2, &["-4", "60", "-180", "140"]);
        max_ulp_overall = max_ulp_overall.max(ulp_refined2);

        // Round-trip after refinement
        let ax = h.mul_vector(&x_refined2);
        let mut max_rt_ulp = 0u64;
        for i in 0..4 {
            let ulp = ulp_diff(ax[i], b[i]);
            max_rt_ulp = max_rt_ulp.max(ulp);
        }
        println!("  {:<40} {:>6} ULP", "Hilbert refined round-trip |Ax-b|", max_rt_ulp);
        max_ulp_overall = max_ulp_overall.max(max_rt_ulp);
    }
    println!();

    // ── SUMMARY ──
    println!("╔══════════════════════════════════════════════════════════════╗");
    println!("║  MAX ULP ACROSS ALL MEASUREMENTS: {:>6}                    ║", max_ulp_overall);
    println!("╚══════════════════════════════════════════════════════════════╝");

    // The max ULP is dominated by the Hilbert matrix (cond ≈ 28375).
    // This is a mathematical inevitability: error ≈ κ(A) × arithmetic_ULP.
    // With κ ≈ 28375 and arithmetic at 1-2 ULP, the solution error is ~50K-2M ULP.
    // The RESIDUAL (Ax-b) stays small (~136 ULP) — the amplification is in the
    // solution space, not the residual space.
    //
    // Well-conditioned bound: single-digit ULP (verified above in the output).
    // Q16.16 (16 fractional bits): QR orthogonality degrades significantly — the
    // Householder reflections accumulate rounding at ~52M ULP. This is the
    // mathematical consequence of only 16 bits of fractional precision, not a bug.
    #[cfg(table_format = "q16_16")]
    let ulp_limit = 100_000_000; // Q16.16: 16-bit fractions → large QR/Hilbert ULP
    #[cfg(not(table_format = "q16_16"))]
    let ulp_limit = 5_000_000;
    assert!(max_ulp_overall < ulp_limit,
        "ULP errors are unreasonably large: {} ULP max (limit {})", max_ulp_overall, ulp_limit);
}