g_math 0.4.2

//! Validation Benchmark Suite — mpmath-validated accuracy + throughput.
//!
//! Combines:
//! - 50K+ mpmath reference points across 4 domains × 4 arithmetic operations
//! - 18 transcendental functions × 1000+ reference points with ULP validation
//! - Ternary UGOD-tier validation with hardcoded test cases
//! - Iteration-based throughput benchmarking with percentile statistics
//!
//! Reuses patterns from existing test suites (UlpStats, ValidationResult,
//! gmath_safe, rational_exact_equal, build_arith_expr) without modifying them.
//!
//! Run per profile:
//!   GMATH_PROFILE=embedded   cargo test --test validation_benchmark -- --nocapture --test-threads=1
//!   GMATH_PROFILE=balanced   cargo test --test validation_benchmark -- --nocapture --test-threads=1
//!   GMATH_PROFILE=scientific cargo test --test validation_benchmark -- --nocapture --test-threads=1

use g_math::fixed_point::canonical::{gmath, evaluate, LazyExpr, set_gmath_mode, reset_gmath_mode};
use g_math::fixed_point::domains::binary_fixed::i256::I256;
#[cfg(table_format = "q256_256")]
use g_math::fixed_point::domains::binary_fixed::i512::I512;
use g_math::fixed_point::domains::symbolic::rational::{RationalNumber, RationalParts};
use std::time::Instant;

// ════════════════════════════════════════════════════════════════════
// Profile detection
// ════════════════════════════════════════════════════════════════════

#[cfg(table_format = "q16_16")]
const ACTIVE_PROFILE: &str = "Q16.16 (realtime)";
#[cfg(table_format = "q32_32")]
const ACTIVE_PROFILE: &str = "Q32.32 (compact)";
#[cfg(table_format = "q64_64")]
const ACTIVE_PROFILE: &str = "Q64.64 (embedded)";
#[cfg(table_format = "q128_128")]
const ACTIVE_PROFILE: &str = "Q128.128 (balanced)";
#[cfg(table_format = "q256_256")]
const ACTIVE_PROFILE: &str = "Q256.256 (scientific)";

#[cfg(table_format = "q16_16")]
const FRAC_BITS: u32 = 16;
#[cfg(table_format = "q32_32")]
const FRAC_BITS: u32 = 32;
#[cfg(table_format = "q64_64")]
const FRAC_BITS: u32 = 64;
#[cfg(table_format = "q128_128")]
const FRAC_BITS: u32 = 128;
#[cfg(table_format = "q256_256")]
const FRAC_BITS: u32 = 256;

// ════════════════════════════════════════════════════════════════════
// Shared infrastructure
// ════════════════════════════════════════════════════════════════════

/// Build a LazyExpr from an input string, correctly handling negative values.
fn gmath_safe(input: &'static str) -> LazyExpr {
    if input.starts_with('-') {
        let positive: &'static str = unsafe {
            std::str::from_utf8_unchecked(
                std::slice::from_raw_parts(input.as_ptr().add(1), input.len() - 1)
            )
        };
        -gmath(positive)
    } else {
        gmath(input)
    }
}

/// Build expression for a binary arithmetic operation.
fn build_arith_expr(a: &'static str, b: &'static str, op: &str) -> LazyExpr {
    match op {
        "add" => gmath_safe(a) + gmath_safe(b),
        "sub" => gmath_safe(a) - gmath_safe(b),
        "mul" => gmath_safe(a) * gmath_safe(b),
        "div" => gmath_safe(a) / gmath_safe(b),
        _ => panic!("unknown op: {}", op),
    }
}

///// Extract (numerator, denominator) as I256 from a RationalNumber.
/// Works for all tiers up to Massive (I256). Returns None for Ultra/Infinite.
fn rational_to_i256_pair(rational: &RationalNumber) -> Option<(I256, I256)> {
    match rational.extract_native() {
        RationalParts::Tiny(n, d) => Some((I256::from_i128(n as i128), I256::from_i128(d as i128))),
        RationalParts::Small(n, d) => Some((I256::from_i128(n as i128), I256::from_i128(d as i128))),
        RationalParts::Medium(n, d) => Some((I256::from_i128(n as i128), I256::from_i128(d as i128))),
        RationalParts::Large(n, d) => Some((I256::from_i128(n as i128), I256::from_i128(d as i128))),
        RationalParts::Huge(n, d) => Some((I256::from_i128(n), I256::from_u128(d))),
        #[cfg(not(feature = "embedded"))]
        RationalParts::Massive(n, d) => Some((n, d)),
        _ => None,
    }
}

/// Compare actual rational == expected (num/den) via cross-multiplication.
/// a_num/a_den == e_num/e_den  ⟺  a_num * e_den == a_den * e_num
/// Uses I256 for tiers ≤ Massive, falls back to I512 for Ultra (Q256.256).
fn rational_exact_equal(rational: &RationalNumber, expected_num: i128, expected_den: i128) -> bool {
    if expected_num == 0 {
        return rational.numerator_i128() == Some(0);
    }
    // Try I256 cross-multiplication first (works for tiers ≤ Massive)
    if let Some((a_num, a_den)) = rational_to_i256_pair(rational) {
        let e_num = I256::from_i128(expected_num);
        let e_den = I256::from_i128(expected_den);
        return a_num * e_den == a_den * e_num;
    }
    // Q256.256 fallback: I512 cross-multiplication for Ultra-tier rationals
    #[cfg(table_format = "q256_256")]
    {
        if let Some((a_num, a_den)) = rational.extract_native().try_as_i512_pair() {
            let e_num = I512::from_i128(expected_num);
            let e_den = I512::from_i128(expected_den);
            return a_num * e_den == a_den * e_num;
        }
    }
    false
}

/// Check if actual rational is approximately equal to expected within relative tolerance.
/// Returns Some(true) if within tolerance, Some(false) if outside, None if extraction failed.
/// Uses pure integer arithmetic — no floats.
///
/// Relative error = |actual - expected| / |expected| = |a_n*e_d - a_d*e_n| / |a_d * e_n|
/// We check: |a_n*e_d - a_d*e_n| * tolerance_den <= |a_d * e_n| * tolerance_num
fn rational_approx_equal(
    rational: &RationalNumber, expected_num: i128, expected_den: i128,
    tol_num: i128, tol_den: i128,
) -> Option<bool> {
    // Try I256 path first (works for tiers ≤ Massive)
    if let Some((a_num, a_den)) = rational_to_i256_pair(rational) {
        let e_num = I256::from_i128(expected_num);
        let e_den = I256::from_i128(expected_den);
        let lhs = a_num * e_den;
        let rhs = a_den * e_num;
        let diff = if lhs > rhs { lhs - rhs } else { rhs - lhs };
        let reference = a_den * e_num;
        let ref_abs = if reference > I256::from_i128(0) { reference } else { I256::from_i128(0) - reference };
        let diff_scaled = diff * I256::from_i128(tol_den);
        let ref_scaled = ref_abs * I256::from_i128(tol_num);
        return Some(diff_scaled <= ref_scaled);
    }

    // Q256.256 fallback: I512 cross-multiplication for Ultra-tier rationals
    #[cfg(table_format = "q256_256")]
    {
        if let Some((a_num, a_den)) = rational.extract_native().try_as_i512_pair() {
            let e_num = I512::from_i128(expected_num);
            let e_den = I512::from_i128(expected_den);
            let lhs = a_num * e_den;
            let rhs = a_den * e_num;
            let diff = if lhs > rhs { lhs - rhs } else { rhs - lhs };
            let reference = a_den * e_num;
            let zero = I512::from_i128(0);
            let ref_abs = if reference > zero { reference } else { zero - reference };
            let diff_scaled = diff * I512::from_i128(tol_den);
            let ref_scaled = ref_abs * I512::from_i128(tol_num);
            return Some(diff_scaled <= ref_scaled);
        }
    }

    None
}

// ════════════════════════════════════════════════════════════════════
// ULP statistics tracker
// ════════════════════════════════════════════════════════════════════

struct UlpStats {
    pub name: String,
    pub max_ulp: u128,
    pub sum_ulp: u128,
    pub count: usize,
    pub worst_label: String,
    pub errors: usize,
    pub ulps: Vec<u128>,
}

impl UlpStats {
    fn new(name: &str) -> Self {
        Self {
            name: name.to_string(),
            max_ulp: 0,
            sum_ulp: 0,
            count: 0,
            worst_label: String::new(),
            errors: 0,
            ulps: Vec::new(),
        }
    }

    fn record(&mut self, ulp: u128, label: &str) {
        self.ulps.push(ulp);
        self.sum_ulp = self.sum_ulp.saturating_add(ulp);
        self.count += 1;
        if ulp > self.max_ulp {
            self.max_ulp = ulp;
            self.worst_label = label.to_string();
        }
    }

    fn record_error(&mut self, _label: &str) {
        self.errors += 1;
    }

    fn p99(&self) -> u128 {
        if self.ulps.is_empty() { return 0; }
        let mut sorted = self.ulps.clone();
        sorted.sort();
        let idx = (sorted.len() as f64 * 0.99).ceil() as usize;
        sorted[idx.min(sorted.len() - 1)]
    }

    fn guaranteed_decimals(&self) -> u32 {
        if self.max_ulp == 0 {
            return (FRAC_BITS as f64 * std::f64::consts::LOG10_2) as u32;
        }
        let total = FRAC_BITS as f64 * std::f64::consts::LOG10_2;
        let loss = (self.max_ulp as f64).log10();
        let g = total - loss;
        if g < 0.0 { 0 } else { g.floor() as u32 }
    }
}

// ════════════════════════════════════════════════════════════════════
// Validation result tracker (correctness + performance)
// ════════════════════════════════════════════════════════════════════

struct ValidationResult {
    pub name: String,
    pub total: usize,
    pub passed: usize,
    pub max_ulp: u128,
    pub failures: Vec<String>,
    pub errors: Vec<String>,
    pub min_ns: u64,
    pub max_ns: u64,
    pub sum_ns: u64,
    pub times_ns: Vec<u64>,
}

impl ValidationResult {
    fn new(name: &str) -> Self {
        Self {
            name: name.to_string(),
            total: 0,
            passed: 0,
            max_ulp: 0,
            failures: Vec::new(),
            errors: Vec::new(),
            min_ns: u64::MAX,
            max_ns: 0,
            sum_ns: 0,
            times_ns: Vec::new(),
        }
    }

    fn record_timing(&mut self, ns: u64) {
        self.times_ns.push(ns);
        self.sum_ns = self.sum_ns.saturating_add(ns);
        if ns < self.min_ns { self.min_ns = ns; }
        if ns > self.max_ns { self.max_ns = ns; }
    }

    fn record_failure(&mut self, label: &str, detail: &str) {
        if self.failures.len() < 10 {
            self.failures.push(format!("{}: {}", label, detail));
        }
    }

    fn record_error(&mut self, label: &str, detail: &str) {
        if self.errors.len() < 10 {
            self.errors.push(format!("{}: {}", label, detail));
        }
    }

    fn avg_ns(&self) -> u64 {
        if self.total == 0 { return 0; }
        self.sum_ns / self.total as u64
    }

    fn throughput_str(&self) -> String {
        let avg = self.avg_ns();
        if avg == 0 { return "N/A".to_string(); }
        let ops_sec = 1_000_000_000u64 / avg;
        if ops_sec >= 1_000_000 {
            format!("{:.1}M ops/s", ops_sec as f64 / 1_000_000.0)
        } else if ops_sec >= 1_000 {
            format!("{:.1}K ops/s", ops_sec as f64 / 1_000.0)
        } else {
            format!("{} ops/s", ops_sec)
        }
    }
}

// ════════════════════════════════════════════════════════════════════
// Include reference data (all via include! — zero new Python generation)
// ════════════════════════════════════════════════════════════════════

// Rational domain references (profile-independent)
include!("data/topology_sweep_decimal.rs");
include!("data/topology_sweep_symbolic.rs");
include!("data/topology_sweep_cross.rs");

// ════════════════════════════════════════════════════════════════════
// Rational validation runner (decimal, symbolic, cross-domain)
// ════════════════════════════════════════════════════════════════════

fn run_rational_validation(
    refs: &[(&'static str, &'static str, i128, i128, &'static str)],
    op: &str,
    name: &str,
) -> ValidationResult {
    let mut result = ValidationResult::new(name);
    for &(a_str, b_str, expected_num, expected_den, label) in refs {
        result.total += 1;
        let start = Instant::now();
        let eval_result = evaluate(&build_arith_expr(a_str, b_str, op));
        let elapsed_ns = start.elapsed().as_nanos() as u64;
        result.record_timing(elapsed_ns);
        match eval_result {
            Ok(val) => match val.to_rational() {
                Ok(rational) => {
                    if rational_exact_equal(&rational, expected_num, expected_den) {
                        result.passed += 1;
                    } else {
                        let actual_num = rational.numerator_i128();
                        let actual_den = rational.denominator_i128();
                        result.record_failure(label, &format!(
                            "{} {} {} = {:?}/{:?}, expected {}/{}",
                            a_str, op, b_str, actual_num, actual_den, expected_num, expected_den
                        ));
                    }
                }
                Err(e) => result.record_error(label, &format!("to_rational: {:?}", e)),
            },
            Err(e) => result.record_error(label, &format!("evaluate: {:?}", e)),
        }
    }
    result
}

// ════════════════════════════════════════════════════════════════════
// Ternary validation — reuses symbolic mpmath refs via ternary:ternary mode
// ════════════════════════════════════════════════════════════════════

/// Ternary sweep validation: routes symbolic reference data through ternary:ternary mode.
///
/// Reveals domain character:
/// - **Exact**: Integer results, power-of-3 denominators (1/3, 1/9, 1/27...)
/// - **Approximate**: Non-power-of-3 fractions (1/7, 1/13...) — within tolerance
/// - **Fail**: Outside tolerance — genuine domain limitation
fn run_ternary_sweep_validation(
    refs: &[(&'static str, &'static str, i128, i128, &'static str)],
    op: &str,
    name: &str,
) -> (ValidationResult, usize, usize) {
    let mut result = ValidationResult::new(name);
    let mut exact_count = 0usize;
    let mut approx_count = 0usize;
    // Ternary tolerance as rational: 1/10^N (generous, profile-aware)
    let tol_den: i128 = match FRAC_BITS {
        64 => 1_000_000_000_000,                // 10^12 — generous for 19-digit precision
        128 => 1_000_000_000_000_000_000_000,   // 10^21 — generous for 38-digit precision
        256 => 1_000_000_000_000_000_000_000,   // 10^21 — generous for 77-digit precision
        _ => 1_000_000,
    };

    set_gmath_mode("ternary:ternary").expect("set ternary mode");

    for &(a_str, b_str, expected_num, expected_den, label) in refs {
        result.total += 1;
        let start = Instant::now();
        let eval_result = evaluate(&build_arith_expr(a_str, b_str, op));
        let elapsed_ns = start.elapsed().as_nanos() as u64;
        result.record_timing(elapsed_ns);

        match eval_result {
            Ok(val) => match val.to_rational() {
                Ok(rational) => {
                    if rational_exact_equal(&rational, expected_num, expected_den) {
                        result.passed += 1;
                        exact_count += 1;
                    } else {
                        // Approximate check via I256 cross-multiplication
                        match rational_approx_equal(&rational, expected_num, expected_den, 1, tol_den) {
                            Some(true) => {
                                result.passed += 1;
                                approx_count += 1;
                            }
                            _ => {
                                result.record_failure(label, &format!(
                                    "{} {} {} outside tolerance (1/{})",
                                    a_str, op, b_str, tol_den
                                ));
                            }
                        }
                    }
                }
                Err(e) => result.record_error(label, &format!("to_rational: {:?}", e)),
            },
            Err(e) => result.record_error(label, &format!("evaluate: {:?}", e)),
        }
    }

    reset_gmath_mode();
    (result, exact_count, approx_count)
}

// ════════════════════════════════════════════════════════════════════
// Mode routing validation — 24 cases × 8 modes = 192 points
// ════════════════════════════════════════════════════════════════════

/// Mode routing test case: (a, b, op, expected_num, expected_den, label)
/// Covers integers, decimal-native, symbolic-native, and ternary-native inputs
/// to verify correctness across all compute:output mode combinations.
const MODE_ROUTING_CASES: &[(&str, &str, &str, i128, i128, &str)] = &[
    // Integer-result cases — MUST be exact in ALL modes
    ("3", "7", "add", 10, 1, "int_add_3_7"),
    ("5", "4", "mul", 20, 1, "int_mul_5_4"),
    ("100", "25", "sub", 75, 1, "int_sub_100_25"),
    ("20", "4", "div", 5, 1, "int_div_20_4"),
    ("1", "1", "add", 2, 1, "int_add_1_1"),
    ("50", "2", "div", 25, 1, "int_div_50_2"),

    // Decimal-native cases — exact in decimal/symbolic compute, approximate in binary/ternary
    ("0.1", "0.2", "add", 3, 10, "dec_add_01_02"),
    ("1.5", "0.5", "add", 2, 1, "dec_add_15_05"),         // integer result
    ("0.25", "4", "mul", 1, 1, "dec_mul_025_4"),           // integer result
    ("0.5", "0.25", "div", 2, 1, "dec_div_05_025"),        // integer result
    ("0.75", "0.25", "sub", 1, 2, "dec_sub_075_025"),
    ("2.5", "0.4", "mul", 1, 1, "dec_mul_25_04"),          // integer result

    // Symbolic-native cases — exact in symbolic compute, approximate in binary/decimal/ternary
    ("1/3", "2/3", "add", 1, 1, "sym_add_1_3_2_3"),       // integer result
    ("1/7", "1/13", "add", 20, 91, "sym_add_1_7_1_13"),
    ("22/7", "7/22", "mul", 1, 1, "sym_mul_reciprocal"),   // integer result
    ("5/6", "1/3", "sub", 1, 2, "sym_sub_5_6_1_3"),
    ("3/4", "2/3", "mul", 1, 2, "sym_mul_3_4_2_3"),
    ("1/3", "1/3", "div", 1, 1, "sym_div_1_3_1_3"),       // integer result

    // Ternary-native cases (power-of-3 denominators) — exact in ternary/symbolic compute
    ("1/3", "1/3", "add", 2, 3, "ter_add_1_3_1_3"),
    ("1/9", "4/9", "add", 5, 9, "ter_add_1_9_4_9"),
    ("1/3", "1/3", "mul", 1, 9, "ter_mul_1_3_1_3"),
    ("2/3", "1/3", "sub", 1, 3, "ter_sub_2_3_1_3"),
    ("1", "3", "div", 1, 3, "ter_div_1_3"),
    ("1/3", "3", "mul", 1, 1, "ter_mul_1_3_3"),           // integer result
];

struct ModeResult {
    mode: String,
    total: usize,
    exact: usize,
    approx: usize,
    lossy: usize,  // domain limitations, not bugs — proves auto:auto superiority
}

/// Mode routing tolerance as rational 1/tol_den (profile-aware)
fn mode_tolerance_den() -> i128 {
    match FRAC_BITS {
        64 => 1_000_000_000_000,                // 10^12
        128 => 1_000_000_000_000_000_000_000,   // 10^21
        256 => 1_000_000_000_000_000_000_000,   // 10^21
        _ => 1_000_000,
    }
}

fn run_mode_validation(mode_str: &str) -> ModeResult {
    if mode_str == "auto:auto" {
        reset_gmath_mode();
    } else {
        set_gmath_mode(mode_str).expect(&format!("set mode {}", mode_str));
    }

    let tol_den = mode_tolerance_den();

    let mut result = ModeResult {
        mode: mode_str.to_string(),
        total: 0,
        exact: 0,
        approx: 0,
        lossy: 0,
    };

    for &(a, b, op, expected_num, expected_den, label) in MODE_ROUTING_CASES {
        result.total += 1;
        let expr = build_arith_expr(a, b, op);
        match evaluate(&expr) {
            Ok(val) => match val.to_rational() {
                Ok(rational) => {
                    if rational_exact_equal(&rational, expected_num, expected_den) {
                        result.exact += 1;
                    } else {
                        // Approximate check via I256 cross-multiplication
                        match rational_approx_equal(&rational, expected_num, expected_den, 1, tol_den) {
                            Some(true) => {
                                result.approx += 1;
                            }
                            _ => {
                                result.lossy += 1;
                                let actual_num = rational.numerator_i128();
                                let actual_den = rational.denominator_i128();
                                eprintln!("║    LOSSY [{}] {}: {} {} {} = {:?}/{:?} expected {}/{}  ║",
                                    mode_str, label, a, op, b, actual_num, actual_den,
                                    expected_num, expected_den);
                            }
                        }
                    }
                }
                Err(e) => {
                    result.lossy += 1;
                    eprintln!("║    LOSSY [{}] {}: to_rational: {:?}  ║", mode_str, label, e);
                }
            },
            Err(e) => {
                result.lossy += 1;
                eprintln!("║    LOSSY [{}] {}: evaluate: {:?}  ║", mode_str, label, e);
            }
        }
    }

    reset_gmath_mode();
    result
}

// ════════════════════════════════════════════════════════════════════
// BENCHMARK: Throughput measurement helper
// ════════════════════════════════════════════════════════════════════

const BENCH_ITERATIONS: usize = 10_000;
const BENCH_WARMUP: usize = 100;

struct BenchResult {
    name: String,
    min_ns: u64,
    avg_ns: u64,
    p99_ns: u64,
    max_ns: u64,
    throughput_str: String,
}

fn bench_op<F: Fn() -> bool>(name: &str, f: F) -> BenchResult {
    let mut times = Vec::with_capacity(BENCH_ITERATIONS);

    // Warmup
    for _ in 0..BENCH_WARMUP {
        let _ = f();
    }

    // Timed iterations
    for _ in 0..BENCH_ITERATIONS {
        let start = Instant::now();
        let _ = f();
        let ns = start.elapsed().as_nanos() as u64;
        times.push(ns);
    }

    times.sort_unstable();
    let min = times[0];
    let max = times[times.len() - 1];
    let sum: u64 = times.iter().sum();
    let avg = sum / times.len() as u64;
    let p99_idx = ((times.len() as f64) * 0.99).ceil() as usize;
    let p99 = times[p99_idx.min(times.len() - 1)];

    let ops_sec = if avg > 0 { 1_000_000_000u64 / avg } else { 0 };
    let throughput = if ops_sec >= 1_000_000 {
        format!("{:.1}M ops/s", ops_sec as f64 / 1_000_000.0)
    } else if ops_sec >= 1_000 {
        format!("{:.1}K ops/s", ops_sec as f64 / 1_000.0)
    } else {
        format!("{} ops/s", ops_sec)
    };

    BenchResult {
        name: name.to_string(),
        min_ns: min,
        avg_ns: avg,
        p99_ns: p99,
        max_ns: max,
        throughput_str: throughput,
    }
}

// ════════════════════════════════════════════════════════════════════
// TEST 1: Arithmetic validation — all domains
// ════════════════════════════════════════════════════════════════════

#[test]
fn validate_arithmetic_all_domains() {
    eprintln!();
    eprintln!("╔═══════════════════════════════════════════════════════════════════════════╗");
    eprintln!("║  ARITHMETIC VALIDATION — {}                               ║", ACTIVE_PROFILE);
    eprintln!("╠═══════════════════════════════════════════════════════════════════════════╣");
    eprintln!("║  {:14} {:5} {:>7} {:>7} {:>10} {:>10} {:>14}  ║", "Domain", "Op", "Points", "Pass", "MaxULP", "ns/op", "Throughput");
    eprintln!("╠═══════════════════════════════════════════════════════════════════════════╣");

    let mut grand_total = 0usize;
    let mut grand_passed = 0usize;
    let mut all_pass = true;

    // ── Decimal domain ──
    let decimal_ops: [(&str, &[(&str, &str, i128, i128, &str)]); 4] = [
        ("add", TOPO_DECIMAL_ADD_REFS),
        ("sub", TOPO_DECIMAL_SUB_REFS),
        ("mul", TOPO_DECIMAL_MUL_REFS),
        ("div", TOPO_DECIMAL_DIV_REFS),
    ];
    for (op, refs) in &decimal_ops {
        let r = run_rational_validation(refs, op, &format!("decimal_{}", op));
        let min = if r.min_ns == u64::MAX { 0 } else { r.min_ns };
        eprintln!("║  {:14} {:5} {:>7} {:>7} {:>10} {:>8}ns {:>14}  ║",
            "decimal", op, r.total, r.passed, "exact", min, r.throughput_str());
        grand_total += r.total;
        grand_passed += r.passed;
        if r.passed != r.total { all_pass = false; }
        for f in &r.failures { eprintln!("║    FAIL: {}  ║", f); }
    }

    // ── Symbolic domain ──
    let symbolic_ops: [(&str, &[(&str, &str, i128, i128, &str)]); 4] = [
        ("add", TOPO_SYMBOLIC_ADD_REFS),
        ("sub", TOPO_SYMBOLIC_SUB_REFS),
        ("mul", TOPO_SYMBOLIC_MUL_REFS),
        ("div", TOPO_SYMBOLIC_DIV_REFS),
    ];
    for (op, refs) in &symbolic_ops {
        let r = run_rational_validation(refs, op, &format!("symbolic_{}", op));
        let min = if r.min_ns == u64::MAX { 0 } else { r.min_ns };
        eprintln!("║  {:14} {:5} {:>7} {:>7} {:>10} {:>8}ns {:>14}  ║",
            "symbolic", op, r.total, r.passed, "exact", min, r.throughput_str());
        grand_total += r.total;
        grand_passed += r.passed;
        if r.passed != r.total { all_pass = false; }
        for f in &r.failures { eprintln!("║    FAIL: {}  ║", f); }
    }

    // ── Cross-domain ──
    let cross_ops: [(&str, &[(&str, &str, i128, i128, &str)]); 4] = [
        ("add", TOPO_CROSS_ADD_REFS),
        ("sub", TOPO_CROSS_SUB_REFS),
        ("mul", TOPO_CROSS_MUL_REFS),
        ("div", TOPO_CROSS_DIV_REFS),
    ];
    for (op, refs) in &cross_ops {
        let r = run_rational_validation(refs, op, &format!("cross_{}", op));
        let min = if r.min_ns == u64::MAX { 0 } else { r.min_ns };
        eprintln!("║  {:14} {:5} {:>7} {:>7} {:>10} {:>8}ns {:>14}  ║",
            "cross-domain", op, r.total, r.passed, "exact", min, r.throughput_str());
        grand_total += r.total;
        grand_passed += r.passed;
        if r.passed != r.total { all_pass = false; }
        for f in &r.failures { eprintln!("║    FAIL: {}  ║", f); }
    }

    // ── Binary domain (profile-specific) ──
    let binary_results = run_binary_arithmetic_validation();
    for br in &binary_results {
        let min = if br.min_ns == u64::MAX { 0 } else { br.min_ns };
        let ulp_str = format!("{}", br.max_ulp);
        eprintln!("║  {:14} {:5} {:>7} {:>7} {:>10} {:>8}ns {:>14}  ║",
            "binary", &br.name[7..], br.total, br.passed, ulp_str, min, br.throughput_str());
        grand_total += br.total;
        grand_passed += br.passed;
        if br.passed != br.total { all_pass = false; }
        for f in &br.failures { eprintln!("║    FAIL: {}  ║", f); }
    }

    // ── Ternary domain (reuses symbolic refs via ternary:ternary mode) ──
    let ternary_sweep_ops: [(&str, &[(&str, &str, i128, i128, &str)]); 4] = [
        ("add", TOPO_SYMBOLIC_ADD_REFS),
        ("sub", TOPO_SYMBOLIC_SUB_REFS),
        ("mul", TOPO_SYMBOLIC_MUL_REFS),
        ("div", TOPO_SYMBOLIC_DIV_REFS),
    ];
    let mut tern_exact_total = 0usize;
    let mut tern_approx_total = 0usize;
    let mut tern_outside_total = 0usize;
    for (op, refs) in &ternary_sweep_ops {
        let (r, exact, approx) = run_ternary_sweep_validation(refs, op, &format!("ternary_{}", op));
        let ea_str = format!("{}/{}", exact, approx);
        let outside = r.total - r.passed;
        let min = if r.min_ns == u64::MAX { 0 } else { r.min_ns };
        eprintln!("║  {:14} {:5} {:>7} {:>7} {:>10} {:>8}ns {:>14}  ║",
            "ternary", op, r.total, r.passed, ea_str, min, r.throughput_str());
        // Ternary outside-tolerance = domain limitation, not bug.
        // Count toward grand total but don't fail the test.
        grand_total += r.total;
        grand_passed += r.passed;
        tern_exact_total += exact;
        tern_approx_total += approx;
        tern_outside_total += outside;
        // Do NOT set all_pass = false for ternary — these are characterization results
        if !r.failures.is_empty() {
            eprintln!("║  {:14} {:>3} outside tolerance (domain limitation, not bug)     ║",
                "", outside);
        }
    }

    eprintln!("╠═══════════════════════════════════════════════════════════════════════════╣");
    let tern_total = tern_exact_total + tern_approx_total + tern_outside_total;
    eprintln!("║  Ternary: {} exact, {} approx, {} outside tolerance (of {})  ║",
        tern_exact_total, tern_approx_total, tern_outside_total, tern_total);
    eprintln!("║  ARITHMETIC TOTAL: {}/{} passed                                        ║",
        grand_passed, grand_total);
    eprintln!("╚═══════════════════════════════════════════════════════════════════════════╝");

    assert!(all_pass,
        "Arithmetic validation: {}/{} passed ({} failures)",
        grand_passed, grand_total, grand_total - grand_passed);
}

/// Maximum allowed ULP for binary fixed-point arithmetic.
#[allow(dead_code)]
const MAX_ARITHMETIC_ULP: u128 = 1;

// Profile-specific binary arithmetic validation
#[cfg(table_format = "q64_64")]
mod binary_q64_64 {
    use super::*;
    include!("data/topology_sweep_q64_64.rs");

    fn ulp_i128(actual: i128, expected: i128) -> u128 {
        (actual - expected).unsigned_abs()
    }

    pub fn run() -> Vec<ValidationResult> {
        let ops: [(&str, &[(&str, &str, i128, &str)]); 4] = [
            ("add", TOPO_BINARY_ADD_REFS),
            ("sub", TOPO_BINARY_SUB_REFS),
            ("mul", TOPO_BINARY_MUL_REFS),
            ("div", TOPO_BINARY_DIV_REFS),
        ];

        let mut results = Vec::new();
        for (op, refs) in &ops {
            let mut result = ValidationResult::new(&format!("binary_{}", op));
            for &(a_str, b_str, expected, label) in *refs {
                result.total += 1;
                let start = Instant::now();
                let eval_result = evaluate(&build_arith_expr(a_str, b_str, op));
                let elapsed_ns = start.elapsed().as_nanos() as u64;
                result.record_timing(elapsed_ns);
                match eval_result {
                    Ok(val) => match val.as_binary_storage() {
                        Some(actual) => {
                            let ulp = ulp_i128(actual, expected);
                            if ulp > result.max_ulp { result.max_ulp = ulp; }
                            if ulp <= MAX_ARITHMETIC_ULP {
                                result.passed += 1;
                            } else {
                                result.record_failure(label, &format!(
                                    "{} {} {} = ULP {}", a_str, op, b_str, ulp
                                ));
                            }
                        }
                        None => result.record_error(label, "as_binary_storage() returned None"),
                    },
                    Err(e) => result.record_error(label, &format!("evaluate: {:?}", e)),
                }
            }
            results.push(result);
        }
        results
    }
}

#[cfg(table_format = "q64_64")]
fn run_binary_arithmetic_validation() -> Vec<ValidationResult> {
    binary_q64_64::run()
}

#[cfg(table_format = "q128_128")]
mod binary_q128_128 {
    use super::*;
    use g_math::fixed_point::domains::binary_fixed::i256::I256;
    include!("data/topology_sweep_q128_128.rs");

    fn ulp_i256(actual: I256, expected_words: [u64; 4]) -> u128 {
        let expected = I256 { words: expected_words };
        let diff = actual - expected;
        let is_neg = diff.words[3] >> 63 == 1;
        let abs_diff = if is_neg { -diff } else { diff };
        if abs_diff.words[2] != 0 || abs_diff.words[3] != 0 {
            return u128::MAX;
        }
        abs_diff.words[0] as u128 | ((abs_diff.words[1] as u128) << 64)
    }

    pub fn run() -> Vec<ValidationResult> {
        let ops: [(&str, &[(&str, &str, [u64; 4], &str)]); 4] = [
            ("add", TOPO_BINARY_ADD_REFS),
            ("sub", TOPO_BINARY_SUB_REFS),
            ("mul", TOPO_BINARY_MUL_REFS),
            ("div", TOPO_BINARY_DIV_REFS),
        ];

        let mut results = Vec::new();
        for (op, refs) in &ops {
            let mut result = ValidationResult::new(&format!("binary_{}", op));
            for &(a_str, b_str, expected_words, label) in *refs {
                result.total += 1;
                let start = Instant::now();
                let eval_result = evaluate(&build_arith_expr(a_str, b_str, op));
                let elapsed_ns = start.elapsed().as_nanos() as u64;
                result.record_timing(elapsed_ns);
                match eval_result {
                    Ok(val) => match val.as_binary_storage() {
                        Some(actual) => {
                            let ulp = ulp_i256(actual, expected_words);
                            if ulp > result.max_ulp { result.max_ulp = ulp; }
                            if ulp <= MAX_ARITHMETIC_ULP {
                                result.passed += 1;
                            } else {
                                result.record_failure(label, &format!(
                                    "{} {} {} = ULP {}", a_str, op, b_str, ulp
                                ));
                            }
                        }
                        None => result.record_error(label, "as_binary_storage() returned None"),
                    },
                    Err(e) => result.record_error(label, &format!("evaluate: {:?}", e)),
                }
            }
            results.push(result);
        }
        results
    }
}

#[cfg(table_format = "q128_128")]
fn run_binary_arithmetic_validation() -> Vec<ValidationResult> {
    binary_q128_128::run()
}

#[cfg(table_format = "q256_256")]
mod binary_q256_256 {
    use super::*;
    use g_math::fixed_point::domains::binary_fixed::i512::I512;
    include!("data/topology_sweep_q256_256.rs");

    fn ulp_i512(actual: I512, expected_words: [u64; 8]) -> u128 {
        let expected = I512 { words: expected_words };
        let diff = actual - expected;
        let is_neg = diff.words[7] >> 63 == 1;
        let abs_diff = if is_neg { -diff } else { diff };
        for i in 2..8 {
            if abs_diff.words[i] != 0 { return u128::MAX; }
        }
        abs_diff.words[0] as u128 | ((abs_diff.words[1] as u128) << 64)
    }

    pub fn run() -> Vec<ValidationResult> {
        let ops: [(&str, &[(&str, &str, [u64; 8], &str)]); 4] = [
            ("add", TOPO_BINARY_ADD_REFS),
            ("sub", TOPO_BINARY_SUB_REFS),
            ("mul", TOPO_BINARY_MUL_REFS),
            ("div", TOPO_BINARY_DIV_REFS),
        ];

        let mut results = Vec::new();
        for (op, refs) in &ops {
            let mut result = ValidationResult::new(&format!("binary_{}", op));
            for &(a_str, b_str, expected_words, label) in *refs {
                result.total += 1;
                let start = Instant::now();
                let eval_result = evaluate(&build_arith_expr(a_str, b_str, op));
                let elapsed_ns = start.elapsed().as_nanos() as u64;
                result.record_timing(elapsed_ns);
                match eval_result {
                    Ok(val) => match val.as_binary_storage() {
                        Some(actual) => {
                            let ulp = ulp_i512(actual, expected_words);
                            if ulp > result.max_ulp { result.max_ulp = ulp; }
                            if ulp <= MAX_ARITHMETIC_ULP {
                                result.passed += 1;
                            } else {
                                result.record_failure(label, &format!(
                                    "{} {} {} = ULP {}", a_str, op, b_str, ulp
                                ));
                            }
                        }
                        None => result.record_error(label, "as_binary_storage() returned None"),
                    },
                    Err(e) => result.record_error(label, &format!("evaluate: {:?}", e)),
                }
            }
            results.push(result);
        }
        results
    }
}

#[cfg(table_format = "q256_256")]
fn run_binary_arithmetic_validation() -> Vec<ValidationResult> {
    binary_q256_256::run()
}

// ════════════════════════════════════════════════════════════════════
// TEST 2: Transcendental validation — 18 functions
// ════════════════════════════════════════════════════════════════════

#[test]
fn validate_transcendentals_all() {
    eprintln!();
    eprintln!("╔═══════════════════════════════════════════════════════════════════════════╗");
    eprintln!("║  TRANSCENDENTAL VALIDATION — {}                           ║", ACTIVE_PROFILE);
    eprintln!("╠═══════════════════════════════════════════════════════════════════════════╣");
    eprintln!("║  {:10} {:>7} {:>7} {:>5} {:>6} {:>10} {:>14}       ║", "Function", "Points", "MaxULP", "p99", "Dec", "ns/op", "Throughput");
    eprintln!("╠═══════════════════════════════════════════════════════════════════════════╣");

    let results = run_transcendental_validation();

    let mut grand_points = 0usize;
    let mut grand_errors = 0usize;
    let mut overall_max_ulp: u128 = 0;
    let mut all_pass = true;

    for r in &results {
        let per_call_ns = if r.stats.count > 0 {
            r.elapsed_ns as f64 / r.stats.count as f64
        } else { 0.0 };
        let ops_per_sec = if r.elapsed_ns > 0 {
            r.stats.count as f64 / (r.elapsed_ns as f64 / 1e9)
        } else { 0.0 };
        let throughput = if ops_per_sec >= 1_000_000.0 {
            format!("{:.1}M ops/s", ops_per_sec / 1_000_000.0)
        } else if ops_per_sec >= 1_000.0 {
            format!("{:.1}K ops/s", ops_per_sec / 1_000.0)
        } else {
            format!("{:.0} ops/s", ops_per_sec)
        };
        let gd = r.stats.guaranteed_decimals();
        eprintln!("║  {:10} {:>7} {:>7} {:>5} {:>5}d {:>8.0}ns {:>14}       ║",
            r.stats.name, r.stats.count, r.stats.max_ulp, r.stats.p99(),
            gd, per_call_ns, throughput);

        grand_points += r.stats.count;
        grand_errors += r.stats.errors;
        if r.stats.max_ulp > overall_max_ulp { overall_max_ulp = r.stats.max_ulp; }
        if r.stats.max_ulp > r.max_allowed_ulp { all_pass = false; }
        if r.stats.errors > r.max_allowed_errors { all_pass = false; }
    }

    eprintln!("╠═══════════════════════════════════════════════════════════════════════════╣");
    eprintln!("║  TRANSCENDENTAL TOTAL: {} points | max_ulp={} | {} errors               ║",
        grand_points, overall_max_ulp, grand_errors);
    eprintln!("╚═══════════════════════════════════════════════════════════════════════════╝");

    for r in &results {
        assert!(
            r.stats.max_ulp <= r.max_allowed_ulp,
            "{} max ULP {} exceeds threshold {}",
            r.stats.name, r.stats.max_ulp, r.max_allowed_ulp
        );
        assert!(
            r.stats.errors <= r.max_allowed_errors,
            "{} had {} errors (max allowed: {})",
            r.stats.name, r.stats.errors, r.max_allowed_errors
        );
    }
    assert!(all_pass);
}

struct TranscendentalResult {
    pub stats: UlpStats,
    pub elapsed_ns: u64,
    pub max_allowed_ulp: u128,
    pub max_allowed_errors: usize,
}

// Profile-specific transcendental validation
#[cfg(table_format = "q64_64")]
mod trans_q64_64 {
    use super::*;
    include!("data/fasc_ulp_refs_q64_64.rs");

    fn eval_unary(input: &'static str, expected: i128, f: fn(LazyExpr) -> LazyExpr) -> Option<u128> {
        let expr = f(gmath_safe(input));
        let actual = evaluate(&expr).ok()?.as_binary_storage()?;
        Some((actual as i128 - expected).unsigned_abs())
    }

    fn eval_binary(a: &'static str, b: &'static str, expected: i128, f: fn(LazyExpr, LazyExpr) -> LazyExpr) -> Option<u128> {
        let expr = f(gmath_safe(a), gmath_safe(b));
        let actual = evaluate(&expr).ok()?.as_binary_storage()?;
        Some((actual as i128 - expected).unsigned_abs())
    }

    fn run_unary(name: &str, refs: &[(&'static str, i128, &'static str)], max_ulp: u128, max_err: usize, method: fn(LazyExpr) -> LazyExpr) -> TranscendentalResult {
        let mut stats = UlpStats::new(name);
        if let Some(&(input, _, _)) = refs.first() { let _ = evaluate(&method(gmath_safe(input))); }
        let start = Instant::now();
        for &(input, expected, label) in refs {
            match eval_unary(input, expected, method) {
                Some(ulp) => stats.record(ulp, label),
                None => stats.record_error(label),
            }
        }
        let elapsed = start.elapsed().as_nanos() as u64;
        TranscendentalResult { stats, elapsed_ns: elapsed, max_allowed_ulp: max_ulp, max_allowed_errors: max_err }
    }

    fn run_binary_fn(name: &str, refs: &[(&'static str, &'static str, i128, &'static str)], max_ulp: u128, f: fn(LazyExpr, LazyExpr) -> LazyExpr) -> TranscendentalResult {
        let mut stats = UlpStats::new(name);
        if let Some(&(a, b, _, _)) = refs.first() { let _ = evaluate(&f(gmath_safe(a), gmath_safe(b))); }
        let start = Instant::now();
        for &(a, b, expected, label) in refs {
            match eval_binary(a, b, expected, f) {
                Some(ulp) => stats.record(ulp, label),
                None => stats.record_error(label),
            }
        }
        let elapsed = start.elapsed().as_nanos() as u64;
        TranscendentalResult { stats, elapsed_ns: elapsed, max_allowed_ulp: max_ulp, max_allowed_errors: 0 }
    }

    pub fn run() -> Vec<TranscendentalResult> {
        vec![
            run_unary("exp",   FASC_EXP_REFS,   5, 0, LazyExpr::exp),
            run_unary("ln",    FASC_LN_REFS,    5, 0, LazyExpr::ln),
            run_unary("sqrt",  FASC_SQRT_REFS,  5, 0, LazyExpr::sqrt),
            run_unary("sin",   FASC_SIN_REFS,   5, 0, LazyExpr::sin),
            run_unary("cos",   FASC_COS_REFS,   5, 0, LazyExpr::cos),
            run_unary("tan",   FASC_TAN_REFS,   5, 0, LazyExpr::tan),
            run_unary("atan",  FASC_ATAN_REFS,  5, 0, LazyExpr::atan),
            run_unary("asin",  FASC_ASIN_REFS,  5, 0, LazyExpr::asin),
            run_unary("acos",  FASC_ACOS_REFS,  5, 0, LazyExpr::acos),
            run_unary("sinh",  FASC_SINH_REFS,  5, 0, LazyExpr::sinh),
            run_unary("cosh",  FASC_COSH_REFS,  5, 0, LazyExpr::cosh),
            run_unary("tanh",  FASC_TANH_REFS,  5, 0, LazyExpr::tanh),
            run_unary("asinh", FASC_ASINH_REFS, 5, 0, LazyExpr::asinh),
            run_unary("acosh", FASC_ACOSH_REFS, 5, 0, LazyExpr::acosh),
            run_unary("atanh", FASC_ATANH_REFS, 5, 0, LazyExpr::atanh),
            run_binary_fn("atan2",    FASC_ATAN2_REFS,          5, |a, b| a.atan2(b)),
            run_binary_fn("pow_int",  FASC_POW_INTEGER_REFS,    5, |a, b| a.pow(b)),
            run_binary_fn("pow_frac", FASC_POW_FRACTIONAL_REFS, 5, |a, b| a.pow(b)),
        ]
    }
}

#[cfg(table_format = "q64_64")]
fn run_transcendental_validation() -> Vec<TranscendentalResult> {
    trans_q64_64::run()
}

#[cfg(table_format = "q128_128")]
mod trans_q128_128 {
    use super::*;
    use g_math::fixed_point::domains::binary_fixed::i256::I256;
    include!("data/fasc_ulp_refs_q128_128.rs");

    fn ulp_i256(actual: I256, expected: I256) -> u128 {
        let diff = actual - expected;
        let is_neg = diff.words[3] >> 63 == 1;
        let abs_diff = if is_neg { -diff } else { diff };
        if abs_diff.words[2] != 0 || abs_diff.words[3] != 0 { return u128::MAX; }
        abs_diff.words[0] as u128 | ((abs_diff.words[1] as u128) << 64)
    }

    fn eval_unary(input: &'static str, expected_words: [u64; 4], f: fn(LazyExpr) -> LazyExpr) -> Option<u128> {
        let expr = f(gmath_safe(input));
        let actual = evaluate(&expr).ok()?.as_binary_storage()?;
        Some(ulp_i256(actual, I256 { words: expected_words }))
    }

    fn eval_binary(a: &'static str, b: &'static str, expected_words: [u64; 4], f: fn(LazyExpr, LazyExpr) -> LazyExpr) -> Option<u128> {
        let expr = f(gmath_safe(a), gmath_safe(b));
        let actual = evaluate(&expr).ok()?.as_binary_storage()?;
        Some(ulp_i256(actual, I256 { words: expected_words }))
    }

    fn run_unary(name: &str, refs: &[(&'static str, [u64; 4], &'static str)], max_ulp: u128, max_err: usize, method: fn(LazyExpr) -> LazyExpr) -> TranscendentalResult {
        let mut stats = UlpStats::new(name);
        if let Some(&(input, _, _)) = refs.first() { let _ = evaluate(&method(gmath_safe(input))); }
        let start = Instant::now();
        for &(input, expected, label) in refs {
            match eval_unary(input, expected, method) {
                Some(ulp) => stats.record(ulp, label),
                None => stats.record_error(label),
            }
        }
        let elapsed = start.elapsed().as_nanos() as u64;
        TranscendentalResult { stats, elapsed_ns: elapsed, max_allowed_ulp: max_ulp, max_allowed_errors: max_err }
    }

    fn run_binary_fn(name: &str, refs: &[(&'static str, &'static str, [u64; 4], &'static str)], max_ulp: u128, f: fn(LazyExpr, LazyExpr) -> LazyExpr) -> TranscendentalResult {
        let mut stats = UlpStats::new(name);
        if let Some(&(a, b, _, _)) = refs.first() { let _ = evaluate(&f(gmath_safe(a), gmath_safe(b))); }
        let start = Instant::now();
        for &(a, b, expected, label) in refs {
            match eval_binary(a, b, expected, f) {
                Some(ulp) => stats.record(ulp, label),
                None => stats.record_error(label),
            }
        }
        let elapsed = start.elapsed().as_nanos() as u64;
        TranscendentalResult { stats, elapsed_ns: elapsed, max_allowed_ulp: max_ulp, max_allowed_errors: 0 }
    }

    pub fn run() -> Vec<TranscendentalResult> {
        vec![
            run_unary("exp",   FASC_EXP_REFS,   5, 0, LazyExpr::exp),
            run_unary("ln",    FASC_LN_REFS,    5, 0, LazyExpr::ln),
            run_unary("sqrt",  FASC_SQRT_REFS,  5, 0, LazyExpr::sqrt),
            run_unary("sin",   FASC_SIN_REFS,   5, 0, LazyExpr::sin),
            run_unary("cos",   FASC_COS_REFS,   5, 0, LazyExpr::cos),
            run_unary("tan",   FASC_TAN_REFS,   5, 0, LazyExpr::tan),
            run_unary("atan",  FASC_ATAN_REFS,  5, 0, LazyExpr::atan),
            run_unary("asin",  FASC_ASIN_REFS,  5, 0, LazyExpr::asin),
            run_unary("acos",  FASC_ACOS_REFS,  5, 0, LazyExpr::acos),
            run_unary("sinh",  FASC_SINH_REFS,  5, 0, LazyExpr::sinh),
            run_unary("cosh",  FASC_COSH_REFS,  5, 0, LazyExpr::cosh),
            run_unary("tanh",  FASC_TANH_REFS,  5, 0, LazyExpr::tanh),
            run_unary("asinh", FASC_ASINH_REFS, 5, 0, LazyExpr::asinh),
            run_unary("acosh", FASC_ACOSH_REFS, 5, 0, LazyExpr::acosh),
            run_unary("atanh", FASC_ATANH_REFS, 5, 530, LazyExpr::atanh),
            run_binary_fn("atan2",    FASC_ATAN2_REFS,          5, |a, b| a.atan2(b)),
            run_binary_fn("pow_int",  FASC_POW_INTEGER_REFS,    5, |a, b| a.pow(b)),
            run_binary_fn("pow_frac", FASC_POW_FRACTIONAL_REFS, 5, |a, b| a.pow(b)),
        ]
    }
}

#[cfg(table_format = "q128_128")]
fn run_transcendental_validation() -> Vec<TranscendentalResult> {
    trans_q128_128::run()
}

#[cfg(table_format = "q256_256")]
mod trans_q256_256 {
    use super::*;
    use g_math::fixed_point::domains::binary_fixed::i512::I512;
    include!("data/fasc_ulp_refs_q256_256.rs");

    fn ulp_i512(actual: I512, expected: I512) -> u128 {
        let diff = actual - expected;
        let is_neg = diff.words[7] >> 63 == 1;
        let abs_diff = if is_neg { -diff } else { diff };
        for i in 2..8 { if abs_diff.words[i] != 0 { return u128::MAX; } }
        abs_diff.words[0] as u128 | ((abs_diff.words[1] as u128) << 64)
    }

    fn eval_unary(input: &'static str, expected_words: [u64; 8], f: fn(LazyExpr) -> LazyExpr) -> Option<u128> {
        let expr = f(gmath_safe(input));
        let actual = evaluate(&expr).ok()?.as_binary_storage()?;
        Some(ulp_i512(actual, I512 { words: expected_words }))
    }

    fn eval_binary(a: &'static str, b: &'static str, expected_words: [u64; 8], f: fn(LazyExpr, LazyExpr) -> LazyExpr) -> Option<u128> {
        let expr = f(gmath_safe(a), gmath_safe(b));
        let actual = evaluate(&expr).ok()?.as_binary_storage()?;
        Some(ulp_i512(actual, I512 { words: expected_words }))
    }

    fn run_unary(name: &str, refs: &[(&'static str, [u64; 8], &'static str)], max_ulp: u128, max_err: usize, method: fn(LazyExpr) -> LazyExpr) -> TranscendentalResult {
        let mut stats = UlpStats::new(name);
        if let Some(&(input, _, _)) = refs.first() { let _ = evaluate(&method(gmath_safe(input))); }
        let start = Instant::now();
        for &(input, expected, label) in refs {
            match eval_unary(input, expected, method) {
                Some(ulp) => stats.record(ulp, label),
                None => stats.record_error(label),
            }
        }
        let elapsed = start.elapsed().as_nanos() as u64;
        TranscendentalResult { stats, elapsed_ns: elapsed, max_allowed_ulp: max_ulp, max_allowed_errors: max_err }
    }

    fn run_binary_fn(name: &str, refs: &[(&'static str, &'static str, [u64; 8], &'static str)], max_ulp: u128, f: fn(LazyExpr, LazyExpr) -> LazyExpr) -> TranscendentalResult {
        let mut stats = UlpStats::new(name);
        if let Some(&(a, b, _, _)) = refs.first() { let _ = evaluate(&f(gmath_safe(a), gmath_safe(b))); }
        let start = Instant::now();
        for &(a, b, expected, label) in refs {
            match eval_binary(a, b, expected, f) {
                Some(ulp) => stats.record(ulp, label),
                None => stats.record_error(label),
            }
        }
        let elapsed = start.elapsed().as_nanos() as u64;
        TranscendentalResult { stats, elapsed_ns: elapsed, max_allowed_ulp: max_ulp, max_allowed_errors: 0 }
    }

    pub fn run() -> Vec<TranscendentalResult> {
        vec![
            run_unary("exp",   FASC_EXP_REFS,   5, 0, LazyExpr::exp),
            run_unary("ln",    FASC_LN_REFS,    5, 0, LazyExpr::ln),
            run_unary("sqrt",  FASC_SQRT_REFS,  5, 0, LazyExpr::sqrt),
            run_unary("sin",   FASC_SIN_REFS,   5, 0, LazyExpr::sin),
            run_unary("cos",   FASC_COS_REFS,   5, 0, LazyExpr::cos),
            run_unary("tan",   FASC_TAN_REFS,  15, 0, LazyExpr::tan),
            run_unary("atan",  FASC_ATAN_REFS,  5, 0, LazyExpr::atan),
            run_unary("asin",  FASC_ASIN_REFS,  5, 0, LazyExpr::asin),
            run_unary("acos",  FASC_ACOS_REFS,  5, 0, LazyExpr::acos),
            run_unary("sinh",  FASC_SINH_REFS,  5, 0, LazyExpr::sinh),
            run_unary("cosh",  FASC_COSH_REFS,  5, 0, LazyExpr::cosh),
            run_unary("tanh",  FASC_TANH_REFS,  5, 0, LazyExpr::tanh),
            run_unary("asinh", FASC_ASINH_REFS, 5, 0, LazyExpr::asinh),
            run_unary("acosh", FASC_ACOSH_REFS, 5, 0, LazyExpr::acosh),
            run_unary("atanh", FASC_ATANH_REFS, 5, 530, LazyExpr::atanh),
            run_binary_fn("atan2",    FASC_ATAN2_REFS,          5, |a, b| a.atan2(b)),
            run_binary_fn("pow_int",  FASC_POW_INTEGER_REFS,    5, |a, b| a.pow(b)),
            run_binary_fn("pow_frac", FASC_POW_FRACTIONAL_REFS, 5, |a, b| a.pow(b)),
        ]
    }
}

#[cfg(table_format = "q256_256")]
fn run_transcendental_validation() -> Vec<TranscendentalResult> {
    trans_q256_256::run()
}

// ── Q32.32 compact profile transcendentals ──
#[cfg(table_format = "q32_32")]
mod trans_q32_32 {
    use super::*;
    include!("data/fasc_ulp_refs_q32.rs");
    fn run_func(name: &str, method: fn(LazyExpr) -> LazyExpr) -> TranscendentalResult {
        let refs: Vec<_> = Q32_REFS.iter().filter(|r| r.3 == name).collect();
        let mut stats = UlpStats::new(name);
        if let Some(r) = refs.first() { let _ = evaluate(&method(gmath_safe(r.0))); }
        let start = Instant::now();
        for r in &refs {
            match evaluate(&method(gmath_safe(r.0))).ok().and_then(|v| v.as_binary_storage()) {
                Some(actual) => { stats.record(((actual as i64) - r.2).unsigned_abs() as u128, r.0); }
                None => stats.record_error(r.0),
            }
        }
        TranscendentalResult { stats, elapsed_ns: start.elapsed().as_nanos() as u64, max_allowed_ulp: 5, max_allowed_errors: 0 }
    }
    pub fn run() -> Vec<TranscendentalResult> {
        vec![
            run_func("exp", LazyExpr::exp), run_func("ln", LazyExpr::ln), run_func("sqrt", LazyExpr::sqrt),
            run_func("sin", LazyExpr::sin), run_func("cos", LazyExpr::cos), run_func("tan", LazyExpr::tan),
            run_func("atan", LazyExpr::atan), run_func("asin", LazyExpr::asin), run_func("acos", LazyExpr::acos),
            run_func("sinh", LazyExpr::sinh), run_func("cosh", LazyExpr::cosh), run_func("tanh", LazyExpr::tanh),
            run_func("asinh", LazyExpr::asinh), run_func("acosh", LazyExpr::acosh), run_func("atanh", LazyExpr::atanh),
        ]
    }
}
#[cfg(table_format = "q32_32")]
fn run_transcendental_validation() -> Vec<TranscendentalResult> { trans_q32_32::run() }

// ── Q16.16 realtime profile transcendentals ──
#[cfg(table_format = "q16_16")]
mod trans_q16_16 {
    use super::*;
    include!("data/fasc_ulp_refs_q16.rs");
    fn run_func(name: &str, method: fn(LazyExpr) -> LazyExpr) -> TranscendentalResult {
        let refs: Vec<_> = Q16_REFS.iter().filter(|r| r.3 == name).collect();
        let mut stats = UlpStats::new(name);
        if let Some(r) = refs.first() { let _ = evaluate(&method(gmath_safe(r.0))); }
        let start = Instant::now();
        for r in &refs {
            match evaluate(&method(gmath_safe(r.0))).ok().and_then(|v| v.as_binary_storage()) {
                Some(actual) => { stats.record(((actual as i32) - r.2).unsigned_abs() as u128, r.0); }
                None => stats.record_error(r.0),
            }
        }
        TranscendentalResult { stats, elapsed_ns: start.elapsed().as_nanos() as u64, max_allowed_ulp: 5, max_allowed_errors: 0 }
    }
    pub fn run() -> Vec<TranscendentalResult> {
        vec![
            run_func("exp", LazyExpr::exp), run_func("ln", LazyExpr::ln), run_func("sqrt", LazyExpr::sqrt),
            run_func("sin", LazyExpr::sin), run_func("cos", LazyExpr::cos), run_func("tan", LazyExpr::tan),
            run_func("atan", LazyExpr::atan), run_func("asin", LazyExpr::asin), run_func("acos", LazyExpr::acos),
            run_func("sinh", LazyExpr::sinh), run_func("cosh", LazyExpr::cosh), run_func("tanh", LazyExpr::tanh),
            run_func("asinh", LazyExpr::asinh), run_func("acosh", LazyExpr::acosh), run_func("atanh", LazyExpr::atanh),
        ]
    }
}
#[cfg(table_format = "q16_16")]
fn run_transcendental_validation() -> Vec<TranscendentalResult> { trans_q16_16::run() }

// ── Binary arithmetic stubs for compact/realtime (no topology sweep ref data) ──
#[cfg(table_format = "q32_32")]
fn run_binary_arithmetic_validation() -> Vec<ValidationResult> { vec![] }
#[cfg(table_format = "q16_16")]
fn run_binary_arithmetic_validation() -> Vec<ValidationResult> { vec![] }

// ════════════════════════════════════════════════════════════════════
// TEST 3: Throughput benchmark — iteration-based
// ════════════════════════════════════════════════════════════════════

#[test]
fn benchmark_throughput() {
    eprintln!();
    eprintln!("╔═══════════════════════════════════════════════════════════════════════════╗");
    eprintln!("║  THROUGHPUT BENCHMARK — {} ({} iterations)            ║", ACTIVE_PROFILE, BENCH_ITERATIONS);
    eprintln!("╠═══════════════════════════════════════════════════════════════════════════╣");
    eprintln!("║  {:20} {:>8} {:>8} {:>8} {:>8} {:>14}       ║", "Operation", "min", "avg", "p99", "max", "throughput");
    eprintln!("╠═══════════════════════════════════════════════════════════════════════════╣");

    let mut results: Vec<BenchResult> = Vec::new();

    // Pre-build all expressions — benchmark measures evaluate() only, not parsing
    let expr_bin_add = gmath("1.5") + gmath("2.5");
    let expr_bin_mul = gmath("3.14159") * gmath("2.71828");
    let expr_dec_add = gmath("0.1") + gmath("0.2");
    let expr_dec_mul = gmath("0.25") * gmath("0.75");
    let expr_sym_add = gmath("1/3") + gmath("1/7");
    let expr_sym_mul = gmath("22/7") * gmath("355/113");
    let expr_exp = gmath("1.5").exp();
    let expr_ln = gmath("2.5").ln();
    let expr_sqrt = gmath("2.0").sqrt();
    let expr_sin = gmath("1.0").sin();
    let expr_cos = gmath("1.0").cos();
    let expr_atan = gmath("0.5").atan();

    // Arithmetic benchmarks — 4 domains × 2 representative ops (add, mul)
    results.push(bench_op("binary_add", || {
        evaluate(&expr_bin_add).is_ok()
    }));
    results.push(bench_op("binary_mul", || {
        evaluate(&expr_bin_mul).is_ok()
    }));
    results.push(bench_op("decimal_add", || {
        evaluate(&expr_dec_add).is_ok()
    }));
    results.push(bench_op("decimal_mul", || {
        evaluate(&expr_dec_mul).is_ok()
    }));
    results.push(bench_op("symbolic_add", || {
        evaluate(&expr_sym_add).is_ok()
    }));
    results.push(bench_op("symbolic_mul", || {
        evaluate(&expr_sym_mul).is_ok()
    }));

    // Ternary benchmarks — must set mode before building AND evaluating
    set_gmath_mode("ternary:ternary").ok();
    let expr_tern_add = gmath("3") + gmath("7");
    let expr_tern_mul = gmath("5") * gmath("4");
    results.push(bench_op("ternary_add", || {
        evaluate(&expr_tern_add).is_ok()
    }));
    results.push(bench_op("ternary_mul", || {
        evaluate(&expr_tern_mul).is_ok()
    }));
    reset_gmath_mode();

    // Transcendental benchmarks — 6 representative functions
    results.push(bench_op("exp", || {
        evaluate(&expr_exp).is_ok()
    }));
    results.push(bench_op("ln", || {
        evaluate(&expr_ln).is_ok()
    }));
    results.push(bench_op("sqrt", || {
        evaluate(&expr_sqrt).is_ok()
    }));
    results.push(bench_op("sin", || {
        evaluate(&expr_sin).is_ok()
    }));
    results.push(bench_op("cos", || {
        evaluate(&expr_cos).is_ok()
    }));
    results.push(bench_op("atan", || {
        evaluate(&expr_atan).is_ok()
    }));

    // Mode routing throughput benchmarks — 11 modes (set once, benchmark evaluate)
    let mode_expr = gmath("1.5") + gmath("2.5");
    let mode_benchmarks: [&str; 11] = [
        "auto:binary", "auto:decimal", "auto:symbolic", "auto:ternary",
        "binary:binary", "binary:decimal", "decimal:decimal",
        "symbolic:decimal", "symbolic:symbolic",
        "ternary:decimal", "ternary:ternary",
    ];
    for mode in &mode_benchmarks {
        set_gmath_mode(mode).ok();
        let name = format!("mode:{}", mode);
        results.push(bench_op(&name, || {
            evaluate(&mode_expr).is_ok()
        }));
        reset_gmath_mode();
    }

    // Print results
    for r in &results {
        eprintln!("║  {:20} {:>6}ns {:>6}ns {:>6}ns {:>6}ns {:>14}       ║",
            r.name, r.min_ns, r.avg_ns, r.p99_ns, r.max_ns, r.throughput_str);
    }

    eprintln!("╚═══════════════════════════════════════════════════════════════════════════╝");
}

// ════════════════════════════════════════════════════════════════════
// TEST 4: Mode routing validation — 8 modes × 24 cases
// ════════════════════════════════════════════════════════════════════

#[test]
fn validate_mode_routing() {
    eprintln!();
    eprintln!("╔═══════════════════════════════════════════════════════════════════════════╗");
    eprintln!("║  MODE ROUTING VALIDATION — {}                               ║", ACTIVE_PROFILE);
    eprintln!("╠═══════════════════════════════════════════════════════════════════════════╣");
    eprintln!("║  {:22} {:>5} {:>5} {:>6} {:>5}                        ║",
        "Mode", "Total", "Exact", "Approx", "Lossy");
    eprintln!("╠═══════════════════════════════════════════════════════════════════════════╣");

    let modes = [
        // Auto compute — optimal precision, output conversion only
        "auto:auto",
        "auto:binary",
        "auto:decimal",
        "auto:symbolic",
        "auto:ternary",
        // Forced compute + matching/cross output
        "binary:binary",
        "binary:decimal",
        "decimal:decimal",
        "symbolic:decimal",
        "symbolic:symbolic",
        "ternary:decimal",
        "ternary:ternary",
    ];

    let mut grand_total = 0usize;
    let mut grand_passed = 0usize;
    let mut grand_lossy = 0usize;
    let mut auto_auto_exact = 0usize;
    let mut auto_auto_total = 0usize;

    for mode in &modes {
        let r = run_mode_validation(mode);
        let passed = r.exact + r.approx;
        // Highlight auto:auto with marker
        let marker = if *mode == "auto:auto" { " <--" } else { "" };
        eprintln!("║  {:22} {:>5} {:>5} {:>6} {:>5}{:4}                        ║",
            r.mode, r.total, r.exact, r.approx, r.lossy, marker);
        grand_total += r.total;
        grand_passed += passed;
        grand_lossy += r.lossy;
        if *mode == "auto:auto" {
            auto_auto_exact = r.exact + r.approx;
            auto_auto_total = r.total;
        }
    }

    eprintln!("╠═══════════════════════════════════════════════════════════════════════════╣");
    eprintln!("║  auto:auto: {}/{} exact — domain-optimal routing, 0 lossy               ║",
        auto_auto_exact, auto_auto_total);
    eprintln!("║  Other modes: {} lossy results = domain limitations (not bugs)           ║",
        grand_lossy);
    eprintln!("║  MODE ROUTING TOTAL: {}/{} validated ({} modes x {} cases)               ║",
        grand_passed + grand_lossy, grand_total, modes.len(), MODE_ROUTING_CASES.len());
    eprintln!("╚═══════════════════════════════════════════════════════════════════════════╝");

    // auto:auto MUST be perfect — this is the library's core value proposition
    assert_eq!(auto_auto_exact, auto_auto_total,
        "auto:auto must achieve 100% exact/approx results, got {}/{}",
        auto_auto_exact, auto_auto_total);
    // All modes must at least not crash — every evaluation must produce a result
    assert_eq!(grand_passed + grand_lossy, grand_total,
        "All evaluations must produce results (no crashes)");
}

// ════════════════════════════════════════════════════════════════════
// TEST 5: Validation summary
// ════════════════════════════════════════════════════════════════════

#[test]
fn validation_summary() {
    let decimal_total =
        TOPO_DECIMAL_ADD_REFS.len() + TOPO_DECIMAL_SUB_REFS.len() +
        TOPO_DECIMAL_MUL_REFS.len() + TOPO_DECIMAL_DIV_REFS.len();
    let symbolic_total =
        TOPO_SYMBOLIC_ADD_REFS.len() + TOPO_SYMBOLIC_SUB_REFS.len() +
        TOPO_SYMBOLIC_MUL_REFS.len() + TOPO_SYMBOLIC_DIV_REFS.len();
    let cross_total =
        TOPO_CROSS_ADD_REFS.len() + TOPO_CROSS_SUB_REFS.len() +
        TOPO_CROSS_MUL_REFS.len() + TOPO_CROSS_DIV_REFS.len();
    let ternary_total = symbolic_total; // Ternary reuses symbolic refs via ternary:ternary mode
    let rational_total = decimal_total + symbolic_total + cross_total;

    eprintln!();
    eprintln!("╔═══════════════════════════════════════════════════════════════════════════╗");
    eprintln!("║  VALIDATION SUMMARY — {}                                  ║", ACTIVE_PROFILE);
    eprintln!("╠═══════════════════════════════════════════════════════════════════════════╣");
    eprintln!("║  Decimal:      {:>6} mpmath-verified reference points                   ║", decimal_total);
    eprintln!("║  Symbolic:     {:>6} mpmath-verified reference points                   ║", symbolic_total);
    eprintln!("║  Cross-Domain: {:>6} mpmath-verified reference points                   ║", cross_total);
    eprintln!("║  Binary:       (profile-specific — see validate_arithmetic_all_domains)  ║");
    eprintln!("║  Ternary:      {:>6} mpmath-verified reference points (exact/approx)     ║", ternary_total);
    eprintln!("║  Rational:     {:>6} total arithmetic reference points                  ║", rational_total);
    eprintln!("║  Transcend.:   18 functions x 1000+ refs (profile-specific binary ULP)   ║");
    eprintln!("║  Mode routing: {:>6} points ({} modes x {} cases)                       ║",
        12 * MODE_ROUTING_CASES.len(), 12, MODE_ROUTING_CASES.len());
    eprintln!("╠═══════════════════════════════════════════════════════════════════════════╣");
    eprintln!("║  Domains: 4 (binary, decimal, symbolic, ternary)                         ║");
    eprintln!("║  Modes:   12 (5 auto:* + 7 explicit compute:output combos)               ║");
    eprintln!("║  Ops:     4 arithmetic + 18 transcendental = 22 total                    ║");
    eprintln!("║  Profile: {} | {} frac bits                              ║", ACTIVE_PROFILE, FRAC_BITS);
    eprintln!("╚═══════════════════════════════════════════════════════════════════════════╝");
}