g_math 0.4.2

//! Domain Arithmetic UGOD-FASC Validation Suite
//!
//! Validates basic arithmetic operations (add, sub, mul, div) across ALL domains
//! and ALL UGOD tiers through the full FASC pipeline, with performance profiling.
//!
//! Domains tested:
//! - Binary:       Integer inputs → as_binary_storage() → ULP compare
//! - Decimal:      Exact decimal inputs → to_rational() → cross-multiply compare
//! - Symbolic:     Fraction inputs → to_rational() → cross-multiply compare
//! - Cross-domain: Mixed inputs → to_rational() → cross-multiply compare
//!
//! Reference data generated by: scripts/generate_domain_arithmetic_refs.py
//!
//! Run per profile:
//!   GMATH_PROFILE=embedded   cargo test --test domain_arithmetic_validation -- --nocapture
//!   GMATH_PROFILE=balanced   cargo test --test domain_arithmetic_validation -- --nocapture
//!   GMATH_PROFILE=scientific cargo test --test domain_arithmetic_validation -- --nocapture

use g_math::fixed_point::canonical::{gmath, evaluate, LazyExpr};
use std::time::Instant;

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

/// Build a LazyExpr from an input string, correctly handling negative values.
/// Mirrors gmath_safe from fasc_ulp_validation.rs.
#[allow(dead_code)]
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 the expression for a binary arithmetic operation.
#[allow(dead_code)]
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),
    }
}

/// Compare two rationals via cross-multiplication: a_num/a_den == b_num/b_den
/// iff a_num * b_den == a_den * b_num.
///
/// Handles unreduced fractions correctly (e.g. 40/10 == 4/1).
#[allow(dead_code)]
fn rationals_equal(
    actual_num: Option<i128>, actual_den: Option<i128>,
    expected_num: i128, expected_den: i128,
) -> bool {
    let (a_num, a_den) = match (actual_num, actual_den) {
        (Some(n), Some(d)) => (n, d),
        _ => return false,
    };
    // Both zero
    if a_num == 0 && expected_num == 0 { return true; }
    // One zero
    if a_num == 0 || expected_num == 0 { return false; }
    // Cross-multiply with overflow check
    match (a_num.checked_mul(expected_den), a_den.checked_mul(expected_num)) {
        (Some(lhs), Some(rhs)) => lhs == rhs,
        _ => false, // overflow → report as mismatch (conservative)
    }
}

/// ULP statistics tracker (reused from fasc_ulp_validation pattern).
#[allow(dead_code)]
struct UlpStats {
    name: String,
    max_ulp: u128,
    sum_ulp: u128,
    count: usize,
    worst_label: String,
    errors: usize,
    ulps: Vec<u128>,
}

#[allow(dead_code)]
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;
        if self.errors <= 5 {
            eprintln!("  eval error: {}", label);
        }
    }

    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 mean_f64(&self) -> f64 {
        if self.count == 0 { return 0.0; }
        self.sum_ulp as f64 / self.count as f64
    }

    fn report(&self, profile: &str) {
        eprintln!(
            "  ULP  {:25} {}: max={:<8} mean={:<8.1} p99={:<8} pts={:<4} {}{}",
            self.name, profile, self.max_ulp, self.mean_f64(), self.p99(),
            self.count,
            if self.worst_label.is_empty() { String::new() } else { format!("worst={}", self.worst_label) },
            if self.errors > 0 { format!(" [{} errors]", self.errors) } else { String::new() }
        );
    }
}

/// Performance timing statistics.
#[allow(dead_code)]
struct TimingStats {
    name: String,
    min_ns: u64,
    max_ns: u64,
    sum_ns: u64,
    count: usize,
    times_ns: Vec<u64>,
}

#[allow(dead_code)]
impl TimingStats {
    fn new(name: &str) -> Self {
        Self {
            name: name.to_string(),
            min_ns: u64::MAX,
            max_ns: 0,
            sum_ns: 0,
            count: 0,
            times_ns: Vec::new(),
        }
    }

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

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

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

    fn report(&self) {
        if self.count == 0 { return; }
        eprintln!(
            "  TIME {:25} | ops={:>4} | min={:>8}ns | avg={:>8}ns | p99={:>8}ns | max={:>8}ns",
            self.name, self.count,
            if self.min_ns == u64::MAX { 0 } else { self.min_ns },
            self.avg(), self.p99(), self.max_ns,
        );
    }
}

/// Exact-match statistics for rational comparison.
#[allow(dead_code)]
struct MatchStats {
    name: String,
    matches: usize,
    mismatches: usize,
    errors: usize,
    first_mismatch: Option<String>,
}

#[allow(dead_code)]
impl MatchStats {
    fn new(name: &str) -> Self {
        Self {
            name: name.to_string(),
            matches: 0,
            mismatches: 0,
            errors: 0,
            first_mismatch: None,
        }
    }

    fn record_match(&mut self) {
        self.matches += 1;
    }

    fn record_mismatch(&mut self, label: &str, detail: &str) {
        self.mismatches += 1;
        if self.first_mismatch.is_none() {
            self.first_mismatch = Some(format!("{}: {}", label, detail));
        }
    }

    fn record_error(&mut self, label: &str) {
        self.errors += 1;
        if self.errors <= 5 {
            eprintln!("  eval error: {}", label);
        }
    }

    fn report(&self, profile: &str) {
        let total = self.matches + self.mismatches;
        eprintln!(
            "  MATCH {:25} {}: {}/{} exact {}{}",
            self.name, profile, self.matches, total,
            if let Some(ref m) = self.first_mismatch { format!("(first mismatch: {})", m) } else { String::new() },
            if self.errors > 0 { format!(" [{} errors]", self.errors) } else { String::new() }
        );
    }
}

// ════════════════════════════════════════════════════════════════════
// Common reference data (decimal, symbolic, cross-domain)
// ════════════════════════════════════════════════════════════════════

include!("data/domain_arith_refs_common.rs");

// ════════════════════════════════════════════════════════════════════
// Rational comparison test runner
// ════════════════════════════════════════════════════════════════════

/// Run a rational comparison test suite with timing.
/// Uses cross-multiplication to handle unreduced fractions.
fn run_rational_test(
    refs: &[(&'static str, &'static str, i128, i128, &'static str)],
    op: &str,
    test_name: &str,
    profile: &str,
) -> (MatchStats, TimingStats) {
    let mut match_stats = MatchStats::new(test_name);
    let mut timing = TimingStats::new(test_name);

    for &(a_str, b_str, expected_num, expected_den, label) in refs.iter() {
        let expr = build_arith_expr(a_str, b_str, op);
        let start = Instant::now();
        let result = evaluate(&expr);
        let elapsed = start.elapsed();
        timing.record(elapsed.as_nanos() as u64);

        match result {
            Ok(val) => {
                match val.to_rational() {
                    Ok(rational) => {
                        let actual_num = rational.numerator_i128();
                        let actual_den = rational.denominator_i128();
                        if rationals_equal(actual_num, actual_den, expected_num, expected_den) {
                            match_stats.record_match();
                        } else {
                            let detail = format!(
                                "expected {}/{} got {:?}/{:?} for {} {} {}",
                                expected_num, expected_den,
                                actual_num, actual_den,
                                a_str, op, b_str,
                            );
                            match_stats.record_mismatch(label, &detail);
                            if match_stats.mismatches <= 3 {
                                eprintln!("  MISMATCH {}: {}", label, detail);
                            }
                        }
                    }
                    Err(e) => {
                        match_stats.record_error(label);
                        if match_stats.errors <= 3 {
                            eprintln!("  to_rational error {}: {:?}", label, e);
                        }
                    }
                }
            }
            Err(e) => {
                match_stats.record_error(label);
                if match_stats.errors <= 3 {
                    eprintln!("  evaluate error {}: {:?}", label, e);
                }
            }
        }
    }

    match_stats.report(profile);
    timing.report();

    (match_stats, timing)
}

// ════════════════════════════════════════════════════════════════════
// Macros for generating test functions
// ════════════════════════════════════════════════════════════════════

/// Generate a rational comparison test (decimal, symbolic, cross-domain).
macro_rules! rational_test {
    ($test_name:ident, $display_name:expr, $refs:ident, $op:expr, $profile:expr,
     $max_mismatches:expr, $max_errors:expr) => {
        #[test]
        fn $test_name() {
            let (match_stats, _timing) = run_rational_test(&$refs, $op, $display_name, $profile);
            assert!(
                match_stats.mismatches <= $max_mismatches,
                "{} {}: {} mismatches (max allowed: {}) (first: {:?})",
                $display_name, $profile,
                match_stats.mismatches, $max_mismatches,
                match_stats.first_mismatch,
            );
            assert!(
                match_stats.errors <= $max_errors,
                "{} {}: {} errors (max allowed: {})",
                $display_name, $profile, match_stats.errors, $max_errors,
            );
        }
    };
}

/// Generate a binary ULP test (profile-specific).
#[allow(unused_macros)]
macro_rules! binary_ulp_test {
    ($test_name:ident, $display_name:expr, $refs:ident, $op:expr, $max_ulp:expr, $profile:expr, $eval_fn:ident) => {
        #[test]
        fn $test_name() {
            let mut ulp_stats = UlpStats::new($display_name);
            let mut timing = TimingStats::new($display_name);

            for &(a_str, b_str, expected, label) in $refs.iter() {
                let expr = build_arith_expr(a_str, b_str, $op);
                let start = Instant::now();
                let result = evaluate(&expr);
                let elapsed = start.elapsed();
                timing.record(elapsed.as_nanos() as u64);

                match result {
                    Ok(val) => {
                        match val.as_binary_storage() {
                            Some(actual) => {
                                let ulp = $eval_fn(actual, expected);
                                ulp_stats.record(ulp, label);
                            }
                            None => {
                                ulp_stats.record_error(label);
                            }
                        }
                    }
                    Err(_) => ulp_stats.record_error(label),
                }
            }

            ulp_stats.report($profile);
            timing.report();

            assert!(
                ulp_stats.max_ulp <= $max_ulp,
                "{} {} max ULP {} exceeds threshold {} (worst: {})",
                $display_name, $profile, ulp_stats.max_ulp, $max_ulp, ulp_stats.worst_label,
            );
        }
    };
}

// ════════════════════════════════════════════════════════════════════
// Detect active profile for reporting
// ════════════════════════════════════════════════════════════════════

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

// ════════════════════════════════════════════════════════════════════
// Decimal domain tests (all profiles — cross-multiply rational compare)
// ════════════════════════════════════════════════════════════════════

// Decimal arithmetic through FASC: result may be Decimal or Symbolic.
// Cross-multiplication handles unreduced fractions (40/10 == 4/1).
// Mixed decimal places may route through rational (→ PrecisionLimit).
// Known issues:
//   - Decimal mul: keeps decimals=1 instead of doubling (e.g. 1.5*2.5 → 3.8 not 3.75)
//   - Decimal div: mixed decimal places may compute wrong result
rational_test!(test_decimal_add, "decimal_add", DECIMAL_ADD_REFS, "add", ACTIVE_PROFILE, 0, 5);
rational_test!(test_decimal_sub, "decimal_sub", DECIMAL_SUB_REFS, "sub", ACTIVE_PROFILE, 0, 5);
rational_test!(test_decimal_mul, "decimal_mul", DECIMAL_MUL_REFS, "mul", ACTIVE_PROFILE, 5, 5);
rational_test!(test_decimal_div, "decimal_div", DECIMAL_DIV_REFS, "div", ACTIVE_PROFILE, 2, 5);

// ════════════════════════════════════════════════════════════════════
// Symbolic domain tests (all profiles — cross-multiply rational compare)
// ════════════════════════════════════════════════════════════════════

// Symbolic (fraction) arithmetic. Cross-multiplication handles unreduced results.
// Allow some mismatches: pre-existing tier arithmetic issues for certain denominators.
// Allow some errors: TierOverflow for large denominators (1/997 + 1/991).
rational_test!(test_symbolic_add, "symbolic_add", SYMBOLIC_ADD_REFS, "add", ACTIVE_PROFILE, 2, 10);
rational_test!(test_symbolic_sub, "symbolic_sub", SYMBOLIC_SUB_REFS, "sub", ACTIVE_PROFILE, 2, 10);
rational_test!(test_symbolic_mul, "symbolic_mul", SYMBOLIC_MUL_REFS, "mul", ACTIVE_PROFILE, 2, 5);
rational_test!(test_symbolic_div, "symbolic_div", SYMBOLIC_DIV_REFS, "div", ACTIVE_PROFILE, 2, 5);

// ════════════════════════════════════════════════════════════════════
// Cross-domain tests (all profiles — cross-multiply rational compare)
// ════════════════════════════════════════════════════════════════════

// Cross-domain arithmetic (Binary+Decimal, Binary+Symbolic, Decimal+Symbolic).
// May hit PrecisionLimit for many combinations — cross-domain ops convert to rational
// internally, and the tier arithmetic may overflow for mixed-domain operands.
rational_test!(test_cross_domain_add, "cross_add", CROSS_ADD_REFS, "add", ACTIVE_PROFILE, 2, 11);
rational_test!(test_cross_domain_sub, "cross_sub", CROSS_SUB_REFS, "sub", ACTIVE_PROFILE, 2, 10);
rational_test!(test_cross_domain_mul, "cross_mul", CROSS_MUL_REFS, "mul", ACTIVE_PROFILE, 2, 10);
rational_test!(test_cross_domain_div, "cross_div", CROSS_DIV_REFS, "div", ACTIVE_PROFILE, 2, 10);

// ════════════════════════════════════════════════════════════════════
// Q64.64 Binary tests (embedded profile)
// ════════════════════════════════════════════════════════════════════

#[cfg(table_format = "q64_64")]
mod q64_64_binary {
    use super::*;

    include!("data/domain_arith_refs_q64_64.rs");

    const PROFILE: &str = "Q64.64";

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

    // Integer add/sub/mul: exact (0 ULP expected)
    binary_ulp_test!(test_binary_add, "binary_add", BINARY_ADD_REFS, "add", 0, PROFILE, ulp_i128);
    binary_ulp_test!(test_binary_sub, "binary_sub", BINARY_SUB_REFS, "sub", 0, PROFILE, ulp_i128);
    binary_ulp_test!(test_binary_mul, "binary_mul", BINARY_MUL_REFS, "mul", 0, PROFILE, ulp_i128);
    // Division: exact for integer-divisible and power-of-2 fractions,
    // potentially large ULP for others due to fixed-point truncation.
    // Report stats but do not assert (binary division precision is a known limitation).
    #[test]
    fn test_binary_div() {
        let mut ulp_stats = UlpStats::new("binary_div");
        let mut timing = TimingStats::new("binary_div");

        for &(a_str, b_str, expected, label) in BINARY_DIV_REFS.iter() {
            let expr = build_arith_expr(a_str, b_str, "div");
            let start = Instant::now();
            let result = evaluate(&expr);
            let elapsed = start.elapsed();
            timing.record(elapsed.as_nanos() as u64);

            match result {
                Ok(val) => {
                    match val.as_binary_storage() {
                        Some(actual) => {
                            let ulp = ulp_i128(actual, expected);
                            ulp_stats.record(ulp, label);
                        }
                        None => ulp_stats.record_error(label),
                    }
                }
                Err(_) => ulp_stats.record_error(label),
            }
        }

        ulp_stats.report(PROFILE);
        timing.report();

        // Report but use generous threshold — binary division precision
        // depends on implementation details (intermediate bit width).
        // Exact divisions (10/2, 1/4, etc.) should be 0 ULP.
        // Inexact (7/3, 1/7) may have significant ULP.
        eprintln!("  NOTE: binary_div ULP reflects fixed-point truncation for inexact results");
    }
}

// ════════════════════════════════════════════════════════════════════
// Q128.128 Binary tests (balanced profile)
// ════════════════════════════════════════════════════════════════════

#[cfg(table_format = "q128_128")]
mod q128_128_binary {
    use super::*;
    use g_math::fixed_point::domains::binary_fixed::i256::I256;

    include!("data/domain_arith_refs_q128_128.rs");

    const PROFILE: &str = "Q128.128";

    fn ulp_i256(actual: I256, expected: [u64; 4]) -> u128 {
        let expected = I256 { words: expected };
        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)
    }

    binary_ulp_test!(test_binary_add, "binary_add", BINARY_ADD_REFS, "add", 0, PROFILE, ulp_i256);
    binary_ulp_test!(test_binary_sub, "binary_sub", BINARY_SUB_REFS, "sub", 0, PROFILE, ulp_i256);
    binary_ulp_test!(test_binary_mul, "binary_mul", BINARY_MUL_REFS, "mul", 0, PROFILE, ulp_i256);

    #[test]
    fn test_binary_div() {
        let mut ulp_stats = UlpStats::new("binary_div");
        let mut timing = TimingStats::new("binary_div");

        for &(a_str, b_str, expected, label) in BINARY_DIV_REFS.iter() {
            let expr = build_arith_expr(a_str, b_str, "div");
            let start = Instant::now();
            let result = evaluate(&expr);
            let elapsed = start.elapsed();
            timing.record(elapsed.as_nanos() as u64);

            match result {
                Ok(val) => {
                    match val.as_binary_storage() {
                        Some(actual) => {
                            let ulp = ulp_i256(actual, expected);
                            ulp_stats.record(ulp, label);
                        }
                        None => ulp_stats.record_error(label),
                    }
                }
                Err(_) => ulp_stats.record_error(label),
            }
        }

        ulp_stats.report(PROFILE);
        timing.report();
        eprintln!("  NOTE: binary_div ULP reflects fixed-point truncation for inexact results");
    }
}

// ════════════════════════════════════════════════════════════════════
// Q256.256 Binary tests (scientific profile)
// ════════════════════════════════════════════════════════════════════

#[cfg(table_format = "q256_256")]
mod q256_256_binary {
    use super::*;
    use g_math::fixed_point::domains::binary_fixed::i512::I512;

    include!("data/domain_arith_refs_q256_256.rs");

    const PROFILE: &str = "Q256.256";

    fn ulp_i512(actual: I512, expected: [u64; 8]) -> u128 {
        let expected = I512 { words: expected };
        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)
    }

    binary_ulp_test!(test_binary_add, "binary_add", BINARY_ADD_REFS, "add", 0, PROFILE, ulp_i512);
    binary_ulp_test!(test_binary_sub, "binary_sub", BINARY_SUB_REFS, "sub", 0, PROFILE, ulp_i512);
    binary_ulp_test!(test_binary_mul, "binary_mul", BINARY_MUL_REFS, "mul", 0, PROFILE, ulp_i512);

    #[test]
    fn test_binary_div() {
        let mut ulp_stats = UlpStats::new("binary_div");
        let mut timing = TimingStats::new("binary_div");

        for &(a_str, b_str, expected, label) in BINARY_DIV_REFS.iter() {
            let expr = build_arith_expr(a_str, b_str, "div");
            let start = Instant::now();
            let result = evaluate(&expr);
            let elapsed = start.elapsed();
            timing.record(elapsed.as_nanos() as u64);

            match result {
                Ok(val) => {
                    match val.as_binary_storage() {
                        Some(actual) => {
                            let ulp = ulp_i512(actual, expected);
                            ulp_stats.record(ulp, label);
                        }
                        None => ulp_stats.record_error(label),
                    }
                }
                Err(_) => ulp_stats.record_error(label),
            }
        }

        ulp_stats.report(PROFILE);
        timing.report();
        eprintln!("  NOTE: binary_div ULP reflects fixed-point truncation for inexact results");
    }
}