mathr 0.1.9

//! Integration tests that exercise multiple modules together end-to-end.

use mathr::complex::Complex;
use mathr::eval::{eval, Context};
use mathr::expr::Expr;
use mathr::fft;
use mathr::interpolate;
use mathr::laurent;
use mathr::matrix::Matrix;
use mathr::notebook::Notebook;
use mathr::numtheory;
use mathr::parser::Parser;
use mathr::ode;
use mathr::rational::{parse_rational, Rational};
use mathr::simplify;
use mathr::special;
use mathr::solver;
use mathr::symbolic;
use mathr::taylor;

fn close(a: f64, b: f64, eps: f64) -> bool {
    (a - b).abs() < eps
}

// =========================================================================
// Expression pipeline: parse → differentiate → simplify → evaluate
// =========================================================================

#[test]
fn parse_diff_simplify_eval_pipeline() {
    // f(x) = x^3, f'(x) = 3x^2, f'(2) = 12
    let expr = Parser::parse("x^3").unwrap();
    let deriv = symbolic::differentiate(&expr, "x").unwrap();
    let simplified = simplify::simplify(&deriv);
    let mut ctx = Context::standard();
    ctx.set("x", 2.0);
    let val = eval(&simplified, &ctx).unwrap();
    assert!(close(val, 12.0, 1e-10));
}

fn square(args: &[f64]) -> mathr::error::Result<f64> {
    Ok(args[0] * args[0])
}

#[test]
fn parse_eval_with_user_function() {
    let mut ctx = Context::standard();
    ctx.insert_builtin("square", square);
    let expr = Parser::parse("square(3) + square(4)").unwrap();
    let val = eval(&expr, &ctx).unwrap();
    assert!(close(val, 25.0, 1e-10));
}

#[test]
fn taylor_series_evaluates_close_to_original() {
    // Taylor series of exp(x) around 0 with 10 terms should be very close at x=1
    let series = taylor::taylor_series_str("exp(x)", "x", 0.0, 10).unwrap();
    let mut ctx = Context::standard();
    ctx.set("x", 1.0);
    let val = eval(&series, &ctx).unwrap();
    assert!(close(val, std::f64::consts::E, 1e-6));
}

// =========================================================================
// FFT: convolution matches direct computation
// =========================================================================

#[test]
fn fft_convolution_matches_direct() {
    let a = vec![1.0, 2.0, 3.0, 4.0];
    let b = vec![5.0, 6.0, 7.0, 8.0];
    let fft_result = fft::convolve(&a, &b).unwrap();

    // Direct computation
    let n = a.len() + b.len() - 1;
    let mut direct = vec![0.0; n];
    for i in 0..a.len() {
        for j in 0..b.len() {
            direct[i + j] += a[i] * b[j];
        }
    }

    for i in 0..n {
        assert!(close(fft_result[i], direct[i], 1e-8));
    }
}

#[test]
fn fft_ifft_roundtrip_preserves_signal() {
    let input: Vec<Complex<f64>> = (0..16)
        .map(|i| Complex::new((i as f64 * 0.5).sin(), (i as f64 * 0.3).cos()))
        .collect();
    let freq = fft::fft(&input).unwrap();
    let recovered = fft::ifft(&freq).unwrap();
    for i in 0..16 {
        assert!(close(input[i].re, recovered[i].re, 1e-10));
        assert!(close(input[i].im, recovered[i].im, 1e-10));
    }
}

#[test]
fn window_function_normalization() {
    let n = 64;
    let signal = vec![1.0; n];
    let hann = fft::apply_window(&signal, fft::Window::Hann);
    let hamming = fft::apply_window(&signal, fft::Window::Hamming);
    let blackman = fft::apply_window(&signal, fft::Window::Blackman);

    // Coherent gain (sum / N) should be close to theoretical values
    // For symmetric windows with (N-1) denominator, exact only in limit
    let hann_gain = hann.iter().sum::<f64>() / n as f64;
    let hamming_gain = hamming.iter().sum::<f64>() / n as f64;
    let blackman_gain = blackman.iter().sum::<f64>() / n as f64;

    assert!(close(hann_gain, 0.5, 1e-2));
    assert!(close(hamming_gain, 0.54, 1e-2));
    assert!(close(blackman_gain, 0.42, 1e-2));
}

// =========================================================================
// Solver: find roots of expressions
// =========================================================================

#[test]
fn solve_polynomial_from_expression() {
    let expr = Parser::parse("x^2 - 4").unwrap();
    let ctx = Context::standard();
    let expr_clone = expr.clone();
    let ctx_clone = ctx.clone();
    let f = move |x: f64| {
        let mut c = ctx_clone.clone();
        c.set("x", x);
        eval(&expr_clone, &c).unwrap_or(f64::NAN)
    };
    let (root, residual) = solver::newton_central(f, 1.5, solver::SolveOptions::default()).unwrap();
    assert!(close(root, 2.0, 1e-8));
    assert!(residual.abs() < 1e-8);
}

#[test]
fn solve_bisection_from_expression() {
    let expr = Parser::parse("x^3 - x - 2").unwrap();
    let ctx = Context::standard();
    let expr_clone = expr.clone();
    let ctx_clone = ctx.clone();
    let f = move |x: f64| {
        let mut c = ctx_clone.clone();
        c.set("x", x);
        eval(&expr_clone, &c).unwrap_or(f64::NAN)
    };
    let (root, residual) = solver::bisect(f, 1.0, 2.0, solver::SolveOptions::default()).unwrap();
    assert!(close(root, 1.5213797068, 1e-6));
    assert!(residual.abs() < 1e-6);
}

// =========================================================================
// Interpolation: recovers polynomial exactly
// =========================================================================

#[test]
fn interpolation_recovers_polynomial() {
    // f(x) = 2x^2 - 3x + 1
    let f = |x: f64| 2.0 * x * x - 3.0 * x + 1.0;
    let points: Vec<(f64, f64)> = vec![(-1.0, f(-1.0)), (0.0, f(0.0)), (1.0, f(1.0)), (2.0, f(2.0))];

    let newton = interpolate::NewtonInterpolator::new(&points).unwrap();
    for x in [-0.5, 0.5, 1.5, 0.3, -0.7] {
        let interp = newton.eval(x);
        let exact = f(x);
        assert!(close(interp, exact, 1e-10));
    }

    // Lagrange should also match
    for x in [-0.5, 0.5, 1.5, 0.3, -0.7] {
        let interp = interpolate::lagrange_interp(&points, x).unwrap();
        let exact = f(x);
        assert!(close(interp, exact, 1e-10));
    }
}

// =========================================================================
// ODE: RK4 accuracy vs analytical solution
// =========================================================================

#[test]
fn rk4_exponential_accuracy() {
    let f = |_t: f64, y: f64| y;
    let result = ode::rk4(f, 0.0, 1.0, 1.0, 1000).unwrap();
    assert!(close(result, std::f64::consts::E, 1e-10));
}

#[test]
fn rk4_system_harmonic_oscillator() {
    let f = |_t: f64, y: &[f64]| vec![y[1], -y[0]];
    let result = ode::rk4_system(f, 0.0, std::f64::consts::FRAC_PI_2, &[1.0, 0.0], 1000).unwrap();
    assert!(close(result[0], 0.0, 1e-8));
    assert!(close(result[1], -1.0, 1e-8));
}

#[test]
fn rkf45_adaptive_accuracy() {
    let f = |_t: f64, y: f64| -3.0 * y;
    let result = ode::rkf45(f, 0.0, 2.0, 1.0, 1e-10).unwrap();
    let exact = (-6.0_f64).exp();
    assert!(close(result, exact, 1e-6));
}

// =========================================================================
// Matrix: solve and inverse consistency
// =========================================================================

#[test]
fn matrix_solve_and_inverse_consistency() {
    let a = Matrix::from_rows(&[
        vec![4.0, 3.0, 2.0],
        vec![1.0, 5.0, 3.0],
        vec![2.0, 1.0, 6.0],
    ]).unwrap();
    let b = vec![20.0, 14.0, 15.0];

    // Solve Ax = b
    let x = a.solve(&b).unwrap();

    // Verify Ax = b
    let ax = a.mul_vec(&x).unwrap();
    for i in 0..3 {
        assert!(close(ax[i], b[i], 1e-8));
    }

    // Inverse * b should give same x
    let inv = a.inverse().unwrap();
    let x2 = inv.mul_vec(&b).unwrap();
    for i in 0..3 {
        assert!(close(x[i], x2[i], 1e-8));
    }
}

#[test]
fn matrix_determinant_and_inverse_properties() {
    let a = Matrix::from_rows(&[
        vec![2.0, 1.0, 0.0],
        vec![1.0, 3.0, 1.0],
        vec![0.0, 1.0, 2.0],
    ]).unwrap();
    let det = a.determinant().unwrap();
    let inv = a.inverse().unwrap();
    let product = (&a * &inv).unwrap();

    // A * A⁻¹ should be identity
    for i in 0..3 {
        for j in 0..3 {
            let expected = if i == j { 1.0 } else { 0.0 };
            assert!(close(product.get(i, j), expected, 1e-8));
        }
    }
    // det(A⁻¹) = 1/det(A)
    let det_inv = inv.determinant().unwrap();
    assert!(close(det_inv, 1.0 / det, 1e-8));
}

// =========================================================================
// Special functions: known values and identities
// =========================================================================

#[test]
fn gamma_known_values() {
    assert!(close(special::gamma(1.0), 1.0, 1e-10));
    assert!(close(special::gamma(2.0), 1.0, 1e-10));
    assert!(close(special::gamma(3.0), 2.0, 1e-10));
    assert!(close(special::gamma(4.0), 6.0, 1e-10));
    assert!(close(special::gamma(5.0), 24.0, 1e-10));
    assert!(close(special::gamma(0.5), std::f64::consts::PI.sqrt(), 1e-10));
}

#[test]
fn gamma_reflection_identity() {
    // Γ(z)Γ(1-z) = π / sin(πz)
    for z in [0.1, 0.25, 0.5, 0.75, 0.9] {
        let lhs = special::gamma(z) * special::gamma(1.0 - z);
        let rhs = std::f64::consts::PI / (std::f64::consts::PI * z).sin();
        assert!(close(lhs, rhs, 1e-8));
    }
}

#[test]
fn beta_function_known_values() {
    assert!(close(special::beta(1.0, 1.0), 1.0, 1e-10));
    assert!(close(special::beta(0.5, 0.5), std::f64::consts::PI, 1e-9));
    // B(2, 3) = Γ(2)Γ(3)/Γ(5) = 1*2/24 = 1/12
    assert!(close(special::beta(2.0, 3.0), 1.0 / 12.0, 1e-10));
}

#[test]
fn erf_known_values() {
    assert!(close(special::erf(0.0), 0.0, 1e-15));
    assert!(close(special::erf(1.0), 0.8427007929, 1e-8));
    assert!(close(special::erf(-1.0), -0.8427007929, 1e-8));
    assert!(close(special::erf(0.5), 0.5204998778, 1e-8));
    // erfc(x) + erf(x) = 1
    for x in [0.0, 0.5, 1.0, 2.0, 3.0] {
        assert!(close(special::erf(x) + special::erfc(x), 1.0, 1e-10));
    }
}

#[test]
fn incomplete_gamma_chi_squared_cdf() {
    // P(1, x) = 1 - e^{-x} is the CDF of χ² with 2 dof
    for x in [0.5, 1.0, 2.0, 5.0, 10.0] {
        let p = special::incomplete_gamma_p(1.0, x);
        let expected = 1.0 - (-x).exp();
        assert!(close(p, expected, 1e-8));
    }
}

// =========================================================================
// Number theory: cross-validation
// =========================================================================

#[test]
fn miller_rabin_matches_trial_division() {
    for n in 2u64..1000 {
        let mr = numtheory::is_prime_miller_rabin(n, 20);
        let trial = numtheory::is_prime(n);
        assert_eq!(mr, trial, "Mismatch at n={}", n);
    }
}

#[test]
fn crt_and_mod_pow_consistency() {
    // Verify CRT solution satisfies all congruences
    let remainders = [3u64, 4, 5];
    let moduli = [5u64, 7, 11];
    let x = numtheory::chinese_remainder(&remainders, &moduli).unwrap();
    for i in 0..3 {
        assert_eq!(x % moduli[i], remainders[i]);
    }
}

#[test]
fn prime_factorization_product_check() {
    for n in [12u64, 60, 360, 1024, 999983, 1234567890] {
        let factors = numtheory::prime_factors(n);
        let product: u64 = factors.iter().product();
        assert_eq!(product, n);
        // All factors should be prime
        for &f in &factors {
            assert!(numtheory::is_prime(f), "{} is not prime", f);
        }
    }
}

#[test]
fn sieve_contains_all_primes() {
    let primes = numtheory::sieve_primes(10000);
    for &p in &primes {
        assert!(numtheory::is_prime_miller_rabin(p, 20), "{} in sieve but not prime", p);
    }
    // Count should match known value: 1229 primes below 10000
    assert_eq!(primes.len(), 1229);
}

// =========================================================================
// Expression equality: canonical form
// =========================================================================

#[test]
fn expr_equality_comprehensive() {
    // a + b == b + a
    let a = Expr::add(Expr::var("a"), Expr::var("b"));
    let b = Expr::add(Expr::var("b"), Expr::var("a"));
    assert!(a.equals(&b));

    // a * b == b * a
    let a = Expr::mul(Expr::var("a"), Expr::var("b"));
    let b = Expr::mul(Expr::var("b"), Expr::var("a"));
    assert!(a.equals(&b));

    // a - b == a + (-b)
    let a = Expr::sub(Expr::var("x"), Expr::var("y"));
    let b = Expr::add(Expr::var("x"), Expr::neg(Expr::var("y")));
    assert!(a.equals(&b));

    // 2 + 3 + x == 5 + x
    let a = Expr::add(Expr::add(Expr::num(2.0), Expr::num(3.0)), Expr::var("x"));
    let b = Expr::add(Expr::num(5.0), Expr::var("x"));
    assert!(a.equals(&b));

    // 2 * 3 * x == 6 * x
    let a = Expr::mul(Expr::mul(Expr::num(2.0), Expr::num(3.0)), Expr::var("x"));
    let b = Expr::mul(Expr::num(6.0), Expr::var("x"));
    assert!(a.equals(&b));

    // x + 1 != x + 2
    let a = Expr::add(Expr::var("x"), Expr::num(1.0));
    let b = Expr::add(Expr::var("x"), Expr::num(2.0));
    assert!(!a.equals(&b));

    // (a + b) + c == c + (a + b)  (flattening)
    let a = Expr::add(Expr::add(Expr::var("a"), Expr::var("b")), Expr::var("c"));
    let b = Expr::add(Expr::var("c"), Expr::add(Expr::var("a"), Expr::var("b")));
    assert!(a.equals(&b));
}

// =========================================================================
// Performance benchmarks (correctness + timing sanity)
// =========================================================================

#[test]
fn fft_large_signal_performance() {
    // 4096-point FFT should complete in well under 1 second
    let n = 4096;
    let samples: Vec<f64> = (0..n)
        .map(|i| (2.0 * std::f64::consts::PI * 64.0 * i as f64 / n as f64).sin())
        .collect();
    let start = std::time::Instant::now();
    let mags = fft::magnitude_spectrum(&samples).unwrap();
    let elapsed = start.elapsed();
    assert_eq!(mags.len(), n / 2 + 1);
    // Peak should be at bin 64
    let peak = mags[64];
    assert!(peak > 1000.0, "expected strong peak at bin 64, got {}", peak);
    assert!(elapsed.as_millis() < 500, "FFT took too long: {:?}", elapsed);
}

#[test]
fn matrix_inverse_large_performance() {
    // 50×50 matrix inverse should complete quickly
    let n = 50;
    let mut data = Vec::with_capacity(n * n);
    for i in 0..n {
        for j in 0..n {
            // Diagonally dominant matrix (guaranteed invertible)
            if i == j {
                data.push(10.0);
            } else {
                data.push(((i as f64 + 1.0) * (j as f64 + 1.0)).sin() * 0.1);
            }
        }
    }
    let a = Matrix::from_row_major(n, n, data).unwrap();
    let start = std::time::Instant::now();
    let inv = a.inverse().unwrap();
    let elapsed = start.elapsed();

    // Verify A * A⁻¹ ≈ I
    let product = (&a * &inv).unwrap();
    for i in 0..n {
        for j in 0..n {
            let expected = if i == j { 1.0 } else { 0.0 };
            assert!((product.get(i, j) - expected).abs() < 1e-6,
                "A*A⁻¹[{}][{}] = {}, expected {}", i, j, product.get(i, j), expected);
        }
    }
    assert!(elapsed.as_millis() < 500, "Matrix inverse took too long: {:?}", elapsed);
}

#[test]
fn matrix_symmetric_eig_decomposition() {
    let a = Matrix::from_rows(&[
        vec![4.0, 1.0, 2.0],
        vec![1.0, 3.0, 0.0],
        vec![2.0, 0.0, 5.0],
    ]).unwrap();
    let (vals, vecs) = a.symmetric_eig().unwrap();

    // Eigenvalues should be sorted ascending.
    assert!(vals[0] <= vals[1] && vals[1] <= vals[2]);

    // Verify A v = λ v for each eigenvector.
    for j in 0..3 {
        let v: Vec<f64> = (0..3).map(|i| vecs[(i, j)]).collect();
        let av = a.mul_vec(&v).unwrap();
        for i in 0..3 {
            assert!((av[i] - vals[j] * v[i]).abs() < 1e-8,
                "A v_{} != λ_{} v_{} at i={}", j, j, j, i);
        }
    }

    // Verify eigenvectors are orthonormal.
    for i in 0..3 {
        for j in 0..3 {
            let dot: f64 = (0..3).map(|k| vecs[(k, i)] * vecs[(k, j)]).sum();
            let expected = if i == j { 1.0 } else { 0.0 };
            assert!((dot - expected).abs() < 1e-8,
                "v_{}·v_{} = {}, expected {}", i, j, dot, expected);
        }
    }
}

#[test]
fn matrix_symmetric_eig_cli() {
    // Test the REPL dispatch for symlig
    let ctx = mathr::eval::Context::standard();
    let result = mathr::repl::dispatch_str("symlig 2 1 | 1 2", ctx).unwrap();
    assert!(result.is_some());
    let output = result.unwrap();
    assert!(output.contains("eigenvalues"));
    // eigenvalues of [[2,1],[1,2]] are 1 and 3
    assert!(output.contains(" 1") || output.contains("[1,") || output.contains("[1 "),
            "output should contain eigenvalue 1: {}", output);
    assert!(output.contains(" 3") || output.contains("[3,") || output.contains(", 3"),
            "output should contain eigenvalue 3: {}", output);
}

#[test]
fn matrix_hessenberg_cli() {
    let ctx = mathr::eval::Context::standard();
    let result = mathr::repl::dispatch_str("hessenberg 1 2 3 | 4 5 6 | 7 8 10", ctx).unwrap();
    assert!(result.is_some());
    let output = result.unwrap();
    assert!(output.contains("Hessenberg"));
    assert!(output.contains("orthogonal"));
}

#[test]
fn matrix_schur_cli() {
    let ctx = mathr::eval::Context::standard();
    let result = mathr::repl::dispatch_str("schur 4 1 2 | 1 3 0 | 2 0 5", ctx).unwrap();
    assert!(result.is_some());
    let output = result.unwrap();
    assert!(output.contains("triangular"));
    assert!(output.contains("orthogonal"));
}

// ===== QR decomposition =====

#[test]
fn qr_repl_square() {
    let ctx = mathr::eval::Context::standard();
    let result = mathr::repl::dispatch_str("qr 12 -51 4 | 6 167 -68 | -4 24 -41", ctx).unwrap().unwrap();
    assert!(result.contains("QR ok"));
    assert!(result.contains("Q ="));
    assert!(result.contains("R ="));
    // Reconstruction error should be tiny
    assert!(result.contains("error = 0.00e+00") || result.contains("error = "));
}

#[test]
fn qr_repl_rectangular_tall() {
    let ctx = mathr::eval::Context::standard();
    let result = mathr::repl::dispatch_str("qr 1 1 | 1 2 | 1 3 | 1 4", ctx).unwrap().unwrap();
    assert!(result.contains("QR ok"));
    assert!(result.contains("Q ="));
    assert!(result.contains("R ="));
}

#[test]
fn qr_repl_wide_matrix() {
    let ctx = mathr::eval::Context::standard();
    let result = mathr::repl::dispatch_str("qr 1 2 3 4 | 5 6 7 8", ctx).unwrap().unwrap();
    assert!(result.contains("QR ok"));
}

#[test]
fn qr_repl_identity() {
    let ctx = mathr::eval::Context::standard();
    let result = mathr::repl::dispatch_str("qr 1 0 0 | 0 1 0 | 0 0 1", ctx).unwrap().unwrap();
    assert!(result.contains("QR ok"));
}

#[test]
fn qr_api_reconstruct() {
    use mathr::matrix::Matrix;
    let a = Matrix::from_rows(&[
        vec![12.0, -51.0, 4.0],
        vec![6.0, 167.0, -68.0],
        vec![-4.0, 24.0, -41.0],
    ]).unwrap();
    let qr = a.qr().unwrap();
    let recon = qr.reconstruct();
    for i in 0..3 {
        for j in 0..3 {
            assert!((a[(i, j)] - recon[(i, j)]).abs() < 1e-10);
        }
    }
}

#[test]
fn qr_api_q_orthogonal() {
    use mathr::matrix::Matrix;
    let a = Matrix::from_rows(&[
        vec![1.0, 2.0],
        vec![3.0, 4.0],
        vec![5.0, 6.0],
    ]).unwrap();
    let qr = a.qr().unwrap();
    let q = qr.q();
    let qt = q.transpose();
    let qtq = (&qt * &q).unwrap();
    for i in 0..3 {
        for j in 0..3 {
            let expected = if i == j { 1.0 } else { 0.0 };
            assert!((qtq[(i, j)] - expected).abs() < 1e-10);
        }
    }
}

#[test]
fn qr_api_least_squares() {
    use mathr::matrix::Matrix;
    // Fit y = 1 + 2x: data (1,3), (2,5), (3,7)
    let a = Matrix::from_rows(&[
        vec![1.0, 1.0],
        vec![1.0, 2.0],
        vec![1.0, 3.0],
    ]).unwrap();
    let b = vec![3.0, 5.0, 7.0];
    let x = a.qr().unwrap().solve(&b).unwrap();
    assert!((x[0] - 1.0).abs() < 1e-10);
    assert!((x[1] - 2.0).abs() < 1e-10);
}

#[test]
fn qr_api_r_upper_triangular() {
    use mathr::matrix::Matrix;
    let a = Matrix::from_rows(&[
        vec![4.0, 1.0, 2.0],
        vec![3.0, 5.0, 1.0],
        vec![1.0, 2.0, 6.0],
    ]).unwrap();
    let qr = a.qr().unwrap();
    let r = qr.r();
    for i in 0..3 {
        for j in 0..i {
            assert!(r[(i, j)].abs() < 1e-10, "R[{},{}] = {} not zero", i, j, r[(i,j)]);
        }
    }
}

#[test]
fn qr_api_rank_deficient_errors() {
    use mathr::matrix::Matrix;
    let a = Matrix::from_rows(&[
        vec![1.0, 2.0],
        vec![1.0, 2.0],
        vec![1.0, 2.0],
    ]).unwrap();
    let b = vec![1.0, 2.0, 3.0];
    assert!(a.qr().unwrap().solve(&b).is_err());
}

#[test]
fn isolate_roots_cli() {
    let ctx = mathr::eval::Context::standard();
    // (x-1)(x-2)(x-3) = x^3 - 6x^2 + 11x - 6
    let result = mathr::repl::dispatch_str("isolate-roots 1 -6 11 -6", ctx).unwrap();
    assert!(result.is_some());
    let output = result.unwrap();
    assert!(output.contains("x"), "output should contain roots: {}", output);
    // Should find 3 roots
    assert!(output.lines().count() >= 3, "expected 3 roots, got: {}", output);
}

#[test]
fn isolate_roots_api() {
    // (x+2)(x-1)(x-3) = x^3 - 2x^2 - 5x + 6
    let intervals = mathr::solver::isolate_real_roots(&[1, -2, -5, 6]).unwrap();
    assert_eq!(intervals.len(), 3);
    let roots = [-2.0, 1.0, 3.0];
    for (i, r) in roots.iter().enumerate() {
        let (lo, hi) = intervals[i];
        assert!(lo <= *r && *r <= hi, "root {} not in ({}, {})", r, lo, hi);
    }
}

#[test]
fn gradient_repl() {
    let ctx = mathr::eval::Context::standard();
    let result = mathr::repl::dispatch_str("gradient x^2 + x*y + y^2", ctx).unwrap();
    assert!(result.is_some());
    let output = result.unwrap();
    assert!(output.contains("d/dx"), "output should contain d/dx: {}", output);
    assert!(output.contains("d/dy"), "output should contain d/dy: {}", output);
}

#[test]
fn pdiff_repl() {
    let ctx = mathr::eval::Context::standard();
    let result = mathr::repl::dispatch_str("pdiff x^2 * y + y^3 y", ctx).unwrap();
    assert!(result.is_some());
    let output = result.unwrap();
    // ∂/∂y (x^2*y + y^3) = x^2 + 3*y^2
    let ctx2 = mathr::eval::Context::standard();
    let mut ctx2 = ctx2;
    ctx2.set("x", 2.0);
    ctx2.set("y", 3.0);
    let e = mathr::parser::Parser::parse(&output).unwrap();
    let val = mathr::eval::eval(&e, &ctx2).unwrap();
    // x^2 + 3*y^2 = 4 + 27 = 31
    assert!((val - 31.0).abs() < 1e-9, "pdiff result evaluated to {} expected 31", val);
}

#[test]
fn fourier_repl() {
    let ctx = mathr::eval::Context::standard();
    // Fourier series of cos(x) on [-pi, pi] with 5 terms, eval at 0
    let result = mathr::repl::dispatch_str("fourier cos(x) 3.14159265358979 5 0", ctx).unwrap();
    assert!(result.is_some());
    let output = result.unwrap();
    assert!(output.contains("a0"), "output should contain a0: {}", output);
    assert!(output.contains("a1"), "output should contain a1: {}", output);
    // cos(x) is a pure Fourier mode, so a1 ≈ 1 and eval at 0 ≈ 1
    assert!(output.contains("f("), "output should contain evaluation: {}", output);
}

#[test]
fn mc_repl() {
    let ctx = mathr::eval::Context::standard();
    // Monte Carlo integral of x over [0, 1] = 0.5
    let result = mathr::repl::dispatch_str("mc x 0 1 100000 42", ctx).unwrap();
    assert!(result.is_some());
    let output = result.unwrap();
    assert!(output.contains("estimate"), "output should contain estimate: {}", output);
    assert!(output.contains("std_error"), "output should contain std_error: {}", output);
}

#[test]
fn sample_repl() {
    let result = mathr::repl::dispatch_str("sample normal 0 1 10000 42", mathr::eval::Context::standard()).unwrap();
    assert!(result.is_some());
    let output = result.unwrap();
    assert!(output.contains("n=10000"), "output should contain n=10000: {}", output);
    assert!(output.contains("mean="), "output should contain mean: {}", output);
}

#[test]
fn dist_repl() {
    let result = mathr::repl::dispatch_str("dist normal 0 0 1", mathr::eval::Context::standard()).unwrap();
    assert!(result.is_some());
    let output = result.unwrap();
    assert!(output.contains("pdf"), "output should contain pdf: {}", output);
    assert!(output.contains("cdf"), "output should contain cdf: {}", output);
    // Normal(0,0,1): pdf = 1/sqrt(2pi), cdf = 0.5
    assert!(output.contains("0.3989"), "pdf should be ~0.3989: {}", output);
    assert!(output.contains("0.5"), "cdf should be 0.5: {}", output);
}

#[test]
fn tikhonov_repl() {
    // Solve a simple 2x2 system with Tikhonov, lambda=0 → same as plain solve
    // A = [[2, 1], [1, 3]], b = [3, 4], x = [1, 1]
    let result = mathr::repl::dispatch_str(
        "tikhonov 2 1 | 1 3 | 3 4 0",
        mathr::eval::Context::standard(),
    )
    .unwrap();
    assert!(result.is_some());
    let output = result.unwrap();
    assert!(output.contains("x = ["), "output should contain x = [...]: {}", output);
    assert!(output.contains("1"), "output should contain solution: {}", output);
}

#[test]
fn spcg_repl() {
    // Solve a simple SPD system: A = [[2, 1], [1, 3]], b = [3, 4], x = [1, 1]
    let result = mathr::repl::dispatch_str(
        "spcg 2 1 | 1 3 | 3 4",
        mathr::eval::Context::standard(),
    )
    .unwrap();
    assert!(result.is_some());
    let output = result.unwrap();
    assert!(output.contains("x = ["), "output should contain x = [...]: {}", output);
    assert!(output.contains("iterations"), "output should report iterations: {}", output);
    assert!(output.contains("residual"), "output should report residual: {}", output);
}

#[test]
fn spcg_repl_nonsymmetric_error() {
    // A = [[2, 1], [0, 3]] is not symmetric — CG must reject it.
    let result = mathr::repl::dispatch_str(
        "spcg 2 1 | 0 3 | 3 4",
        mathr::eval::Context::standard(),
    );
    assert!(result.is_err(), "nonsymmetric input should error");
}

#[test]
fn spcg_jacobi_repl() {
    // Jacobi preconditioning solves the same SPD system.
    let result = mathr::repl::dispatch_str(
        "spcg jacobi 2 1 | 1 3 | 3 4",
        mathr::eval::Context::standard(),
    )
    .unwrap();
    assert!(result.is_some());
    let output = result.unwrap();
    assert!(output.contains("x = ["), "output should contain x = [...]: {}", output);
    assert!(output.contains("iterations"), "output should report iterations: {}", output);
}

#[test]
fn spbicg_repl() {
    // Nonsymmetric system [[2,1],[0,3]]·x = [3,4] has x = [5/6, 4/3].
    let result = mathr::repl::dispatch_str(
        "spbicg 2 1 | 0 3 | 3 4",
        mathr::eval::Context::standard(),
    )
    .unwrap();
    assert!(result.is_some());
    let output = result.unwrap();
    assert!(output.contains("x = ["), "output should contain x = [...]: {}", output);
    assert!(
        output.contains("0.83333"),
        "x[0] should be ~5/6: {}",
        output
    );
}

#[test]
fn spbicg_repl_singular_error() {
    // Singular inconsistent system must break down, not return garbage.
    let result = mathr::repl::dispatch_str(
        "spbicg 1 1 | 1 1 | 1 2",
        mathr::eval::Context::standard(),
    );
    assert!(result.is_err(), "singular inconsistent input should error");
}

#[test]
fn bspline_repl() {
    // Cubic through (0,0) (1,1) (2,4) (3,9): y = x² is reproduced exactly.
    let result = mathr::repl::dispatch_str(
        "bspline 0 0 1 1 2 4 3 9 1.5",
        mathr::eval::Context::standard(),
    )
    .unwrap();
    assert!(result.is_some());
    let output = result.unwrap();
    assert!(
        output.contains("bspline(1.5) = 2.25"),
        "cubic should reproduce x² exactly: {}",
        output
    );
    assert!(output.contains("degree 3"), "should report degree: {}", output);
}

#[test]
fn hermite_repl() {
    // Slopes matching y = x² on [0, 1] give exact reproduction.
    let result = mathr::repl::dispatch_str(
        "hermite 0 0 0 1 1 2 0.5",
        mathr::eval::Context::standard(),
    )
    .unwrap();
    assert!(result.is_some());
    let output = result.unwrap();
    assert!(
        output.contains("hermite(0.5) = 0.25"),
        "hermite should reproduce x² exactly: {}",
        output
    );
}

#[test]
fn hermite_repl_bad_input() {
    // Triples + x_at: wrong token count must error.
    let result = mathr::repl::dispatch_str(
        "hermite 0 0 0 1 1 0.5",
        mathr::eval::Context::standard(),
    );
    assert!(result.is_err(), "wrong token count should error");
}

#[test]
fn spbicg_ilu_repl() {
    // ILU(0) preconditioning solves the same nonsymmetric system.
    let result = mathr::repl::dispatch_str(
        "spbicg ilu 2 1 | 0 3 | 3 4",
        mathr::eval::Context::standard(),
    )
    .unwrap();
    assert!(result.is_some());
    let output = result.unwrap();
    assert!(output.contains("x = ["), "output should contain x = [...]: {}", output);
    assert!(output.contains("0.83333"), "x[0] should be ~5/6: {}", output);
}

#[test]
fn spbicg_ilu_repl_zero_pivot_error() {
    // [[1,1],[1,1]] elimination produces a zero pivot — ILU(0) must error.
    let result = mathr::repl::dispatch_str(
        "spbicg ilu 1 1 | 1 1 | 1 2",
        mathr::eval::Context::standard(),
    );
    assert!(result.is_err(), "zero-pivot input should error");
}

#[test]
fn laurent_repl() {
    // Laurent series of 1/x around x=0, pole order 1, 3 positive terms
    let result = mathr::repl::dispatch_str(
        "laurent 1/x 0 1 3",
        mathr::eval::Context::standard(),
    )
    .unwrap();
    assert!(result.is_some());
    let output = result.unwrap();
    assert!(output.contains("1/x"), "output should contain 1/x: {}", output);
}

#[test]
fn rat_repl() {
    // 1/2 + 1/3 = 5/6
    let result = mathr::repl::dispatch_str(
        "rat 1/2 + 1/3",
        mathr::eval::Context::standard(),
    )
    .unwrap();
    assert!(result.is_some());
    let output = result.unwrap();
    assert!(output.contains("5/6"), "output should contain 5/6: {}", output);
}

#[test]
fn next_pow2_correctness() {
    assert_eq!(fft::next_pow2(0), 1);
    assert_eq!(fft::next_pow2(1), 1);
    assert_eq!(fft::next_pow2(2), 2);
    assert_eq!(fft::next_pow2(3), 4);
    assert_eq!(fft::next_pow2(4), 4);
    assert_eq!(fft::next_pow2(5), 8);
    assert_eq!(fft::next_pow2(7), 8);
    assert_eq!(fft::next_pow2(8), 8);
    assert_eq!(fft::next_pow2(9), 16);
    assert_eq!(fft::next_pow2(1023), 1024);
    assert_eq!(fft::next_pow2(1024), 1024);
    assert_eq!(fft::next_pow2(1025), 2048);
}

#[test]
fn convolution_large_signal_performance() {
    // Convolve two 1024-length signals via FFT
    let a: Vec<f64> = (0..1024).map(|i| (i as f64 * 0.01).sin()).collect();
    let b: Vec<f64> = (0..1024).map(|i| (i as f64 * 0.02).cos()).collect();
    let start = std::time::Instant::now();
    let result = fft::convolve(&a, &b).unwrap();
    let elapsed = start.elapsed();
    assert_eq!(result.len(), 2047);
    assert!(elapsed.as_millis() < 200, "Convolution took too long: {:?}", elapsed);
}

// =========================================================================
// TeX / Markdown input: parse → eval, parse → diff, parse → solve
// =========================================================================

#[test]
fn tex_eval_fraction() {
    let e = Parser::parse(r"\frac{1}{2} + \frac{3}{4}").unwrap();
    let ctx = Context::standard();
    let v = eval(&e, &ctx).unwrap();
    assert!(close(v, 1.25, 1e-10));
}

#[test]
fn tex_eval_sqrt() {
    let e = Parser::parse(r"\sqrt{16}").unwrap();
    let ctx = Context::standard();
    let v = eval(&e, &ctx).unwrap();
    assert!(close(v, 4.0, 1e-10));
}

#[test]
fn tex_eval_sin_pi() {
    let e = Parser::parse(r"\sin(\pi / 4)").unwrap();
    let ctx = Context::standard();
    let v = eval(&e, &ctx).unwrap();
    assert!(close(v, std::f64::consts::FRAC_PI_4.sin(), 1e-10));
}

#[test]
fn tex_eval_cdot() {
    let e = Parser::parse(r"2 \cdot 3 + 4").unwrap();
    let ctx = Context::standard();
    let v = eval(&e, &ctx).unwrap();
    assert!(close(v, 10.0, 1e-10));
}

#[test]
fn tex_eval_left_right() {
    let e = Parser::parse(r"\left( 1 + 2 \right) \cdot 3").unwrap();
    let ctx = Context::standard();
    let v = eval(&e, &ctx).unwrap();
    assert!(close(v, 9.0, 1e-10));
}

#[test]
fn tex_eval_gamma() {
    let e = Parser::parse(r"\Gamma{0.5}").unwrap();
    let ctx = Context::standard();
    let v = eval(&e, &ctx).unwrap();
    assert!(close(v, std::f64::consts::PI.sqrt(), 1e-8));
}

#[test]
fn tex_eval_log_subscript() {
    let e = Parser::parse(r"\log_2{8}").unwrap();
    let ctx = Context::standard();
    let v = eval(&e, &ctx).unwrap();
    assert!(close(v, 3.0, 1e-10));
}

#[test]
fn tex_eval_operatorname() {
    let e = Parser::parse(r"\operatorname{erf}(1.0)").unwrap();
    let ctx = Context::standard();
    let v = eval(&e, &ctx).unwrap();
    assert!(close(v, 0.8427007929, 1e-8));
}

#[test]
fn tex_eval_nested_frac() {
    let e = Parser::parse(r"\frac{\frac{1}{2}}{\frac{3}{4}}").unwrap();
    let ctx = Context::standard();
    let v = eval(&e, &ctx).unwrap();
    assert!(close(v, 2.0 / 3.0, 1e-10));
}

#[test]
fn tex_eval_implicit_mult() {
    let e = Parser::parse(r"2\pi").unwrap();
    let ctx = Context::standard();
    let v = eval(&e, &ctx).unwrap();
    assert!(close(v, 2.0 * std::f64::consts::PI, 1e-10));
}

#[test]
fn markdown_inline_eval() {
    let e = Parser::parse(r"$\sin(\pi / 4)$").unwrap();
    let ctx = Context::standard();
    let v = eval(&e, &ctx).unwrap();
    assert!(close(v, std::f64::consts::FRAC_PI_4.sin(), 1e-10));
}

#[test]
fn markdown_display_eval() {
    let e = Parser::parse(r"$$\frac{1}{2} + \frac{3}{4}$$").unwrap();
    let ctx = Context::standard();
    let v = eval(&e, &ctx).unwrap();
    assert!(close(v, 1.25, 1e-10));
}

#[test]
fn latex_bracket_eval() {
    let e = Parser::parse(r"\[\sqrt{16} + \cos(0)\]").unwrap();
    let ctx = Context::standard();
    let v = eval(&e, &ctx).unwrap();
    assert!(close(v, 5.0, 1e-10));
}

#[test]
fn latex_paren_eval() {
    let e = Parser::parse(r"\(\log_2{8} + 1\)").unwrap();
    let ctx = Context::standard();
    let v = eval(&e, &ctx).unwrap();
    assert!(close(v, 4.0, 1e-10));
}

#[test]
fn tex_diff_matches_plain() {
    let tex_expr = Parser::parse(r"\sin(x^2)").unwrap();
    let plain_expr = Parser::parse("sin(x^2)").unwrap();
    let d_tex = symbolic::differentiate(&tex_expr, "x").unwrap();
    let d_plain = symbolic::differentiate(&plain_expr, "x").unwrap();
    assert_eq!(d_tex.canonicalize(), d_plain.canonicalize());
}

#[test]
fn tex_diff_frac() {
    // d/dx (x^2 + 1)/(x - 1) — parse as TeX, differentiate, evaluate
    let e = Parser::parse(r"\frac{x^2 + 1}{x - 1}").unwrap();
    let d = symbolic::differentiate(&e, "x").unwrap();
    let s = simplify::simplify(&d);
    // At x=2: f(x) = (4+1)/(2-1) = 5, f'(x) = (2x(x-1) - (x^2+1)) / (x-1)^2
    // f'(2) = (4*1 - 5) / 1 = -1
    let mut ctx = Context::standard();
    ctx.set("x", 2.0);
    let v = eval(&s, &ctx).unwrap();
    assert!(close(v, -1.0, 1e-8));
}

#[test]
fn tex_solve_from_frac() {
    let e = Parser::parse(r"\frac{x^2 - 4}{1}").unwrap();
    let ctx = Context::standard();
    let expr_clone = e.clone();
    let ctx_clone = ctx.clone();
    let f = move |x: f64| {
        let mut c = ctx_clone.clone();
        c.set("x", x);
        eval(&expr_clone, &c).unwrap_or(f64::NAN)
    };
    let (root, residual) = solver::newton_central(f, 1.0, solver::SolveOptions::default()).unwrap();
    assert!(close(root, 2.0, 1e-8));
    assert!(residual.abs() < 1e-8);
}

#[test]
fn tex_simplify_matches_plain() {
    let tex_e = Parser::parse(r"\frac{x^2 - 1}{x - 1}").unwrap();
    let plain_e = Parser::parse("(x^2 - 1)/(x - 1)").unwrap();
    let s_tex = simplify::simplify(&tex_e);
    let s_plain = simplify::simplify(&plain_e);
    assert_eq!(s_tex.canonicalize(), s_plain.canonicalize());
}

#[test]
fn tex_taylor_matches_plain() {
    let tex_series = taylor::taylor_series_str(r"\exp(x)", "x", 0.0, 5).unwrap();
    let plain_series = taylor::taylor_series_str("exp(x)", "x", 0.0, 5).unwrap();
    assert_eq!(tex_series.canonicalize(), plain_series.canonicalize());
}

#[test]
fn tex_taylor_evaluates_correctly() {
    let series = taylor::taylor_series_str(r"\sin(x)", "x", 0.0, 7).unwrap();
    let mut ctx = Context::standard();
    ctx.set("x", 0.5);
    let v = eval(&series, &ctx).unwrap();
    assert!(close(v, 0.5_f64.sin(), 1e-4));
}

#[test]
fn tex_pow_brace_eval() {
    let e = Parser::parse(r"2^{10}").unwrap();
    let ctx = Context::standard();
    let v = eval(&e, &ctx).unwrap();
    assert!(close(v, 1024.0, 1e-10));
}

#[test]
fn tex_mixed_expression_eval() {
    // 3\sqrt{16} + \frac{1}{2} = 12.5
    let e = Parser::parse(r"3\sqrt{16} + \frac{1}{2}").unwrap();
    let ctx = Context::standard();
    let v = eval(&e, &ctx).unwrap();
    assert!(close(v, 12.5, 1e-10));
}

#[test]
fn tex_left_right_with_frac() {
    let e = Parser::parse(r"\left( \frac{1}{2} + \frac{1}{2} \right) \cdot 10").unwrap();
    let ctx = Context::standard();
    let v = eval(&e, &ctx).unwrap();
    assert!(close(v, 10.0, 1e-10));
}

#[test]
fn tex_all_delimiters_produce_same_ast() {
    let exprs = [
        Parser::parse(r"\frac{1}{2} + \frac{3}{4}").unwrap(),
        Parser::parse(r"$\frac{1}{2} + \frac{3}{4}$").unwrap(),
        Parser::parse(r"$$\frac{1}{2} + \frac{3}{4}$$").unwrap(),
        Parser::parse(r"\(\frac{1}{2} + \frac{3}{4}\)").unwrap(),
        Parser::parse(r"\[\frac{1}{2} + \frac{3}{4}\]").unwrap(),
    ];
    for i in 1..exprs.len() {
        assert_eq!(exprs[0].canonicalize(), exprs[i].canonicalize());
    }
}

#[test]
fn tex_text_variable_eval() {
    let e = Parser::parse(r"\text{alpha} + 1").unwrap();
    let mut ctx = Context::standard();
    ctx.set("alpha", 2.0);
    let v = eval(&e, &ctx).unwrap();
    assert!(close(v, 3.0, 1e-10));
}

// =========================================================================
// Laurent series: parse → expand → evaluate
// =========================================================================

#[test]
fn laurent_simple_pole_eval() {
    // f(x) = 1/x, Laurent series should give c_-1 = 1, rest = 0
    let ls = laurent::laurent_series_str("1/x", "x", 0.0, 1, 3).unwrap();
    assert!(close(ls.coeff(-1), 1.0, 1e-6));
    assert!(close(ls.coeff(0), 0.0, 1e-6));
    // eval at x=5: 1/5 = 0.2
    assert!(close(ls.eval(5.0), 0.2, 1e-6));
}

#[test]
fn laurent_double_pole_eval() {
    // f(x) = 1/x^2, c_-2 = 1, rest = 0
    let ls = laurent::laurent_series_str("1/x^2", "x", 0.0, 2, 3).unwrap();
    assert!(close(ls.coeff(-2), 1.0, 1e-4));
    // eval at x=4: 1/16 = 0.0625
    assert!(close(ls.eval(4.0), 0.0625, 1e-4));
}

#[test]
fn laurent_exp_over_x_pipeline() {
    // f(x) = exp(x)/x = 1/x + 1 + x/2 + x^2/6 + ...
    let ls = laurent::laurent_series_str("exp(x)/x", "x", 0.0, 1, 6).unwrap();
    assert!(close(ls.coeff(-1), 1.0, 1e-6));
    assert!(close(ls.coeff(0), 1.0, 1e-6));
    assert!(close(ls.coeff(1), 0.5, 1e-6));
    assert!(close(ls.coeff(2), 1.0 / 6.0, 1e-6));
    // eval at x=1: e ≈ 2.71828
    assert!(close(ls.eval(1.0), std::f64::consts::E, 1e-2));
}

#[test]
fn laurent_around_nonzero_eval() {
    // f(x) = 1/(x-2), pole at x=2
    let ls = laurent::laurent_series_str("1/(x-2)", "x", 2.0, 1, 3).unwrap();
    assert!(close(ls.coeff(-1), 1.0, 1e-6));
    // eval at x=5: 1/3
    assert!(close(ls.eval(5.0), 1.0 / 3.0, 1e-4));
}

#[test]
fn laurent_no_pole_matches_taylor() {
    // f(x) = exp(x), no pole → Laurent with pole_order=0 should match Taylor
    let ls = laurent::laurent_series_str("exp(x)", "x", 0.0, 0, 5).unwrap();
    assert_eq!(ls.pole_order, 0);
    assert!(close(ls.coeff(0), 1.0, 1e-6));
    assert!(close(ls.coeff(1), 1.0, 1e-6));
    assert!(close(ls.coeff(2), 0.5, 1e-6));
}

#[test]
fn laurent_rational_function() {
    // f(x) = 1/(x(1-x)) = 1/x + 1 + x + x^2 + ...
    let ls = laurent::laurent_series_str("1/(x*(1-x))", "x", 0.0, 1, 5).unwrap();
    assert!(close(ls.coeff(-1), 1.0, 1e-6));
    assert!(close(ls.coeff(0), 1.0, 1e-6));
    assert!(close(ls.coeff(1), 1.0, 1e-6));
    assert!(close(ls.coeff(2), 1.0, 1e-6));
}

#[test]
fn laurent_to_string_contains_principal() {
    let ls = laurent::laurent_series_str("1/x + 2 + x", "x", 0.0, 1, 2).unwrap();
    let s = ls.to_string();
    assert!(s.contains("1/x"), "string should contain 1/x: {}", s);
}

// =========================================================================
// Rational arithmetic: parse → arithmetic → verify exactness
// =========================================================================

#[test]
fn rational_add_exact() {
    let a = parse_rational("1/2").unwrap();
    let b = parse_rational("1/3").unwrap();
    let c = a + b;
    assert_eq!(c.num(), 5);
    assert_eq!(c.den(), 6);
}

#[test]
fn rational_sub_exact() {
    let a = parse_rational("1/2").unwrap();
    let b = parse_rational("1/3").unwrap();
    let c = a - b;
    assert_eq!(c.num(), 1);
    assert_eq!(c.den(), 6);
}

#[test]
fn rational_mul_exact() {
    let a = Rational::new(2, 3).unwrap();
    let b = Rational::new(3, 4).unwrap();
    let c = a * b;
    assert_eq!(c.num(), 1);
    assert_eq!(c.den(), 2);
}

#[test]
fn rational_div_exact() {
    let a = Rational::new(2, 3).unwrap();
    let b = Rational::new(4, 5).unwrap();
    let c = a / b;
    assert_eq!(c.num(), 5);
    assert_eq!(c.den(), 6);
}

#[test]
fn rational_powi_exact() {
    let a = Rational::new(2, 3).unwrap();
    assert_eq!(a.powi(3), Rational::new(8, 27).unwrap());
    assert_eq!(a.powi(-2), Rational::new(9, 4).unwrap());
    assert_eq!(a.powi(0), Rational::from_int(1));
}

#[test]
fn rational_parse_decimal() {
    let r = parse_rational("0.5").unwrap();
    assert_eq!(r.num(), 1);
    assert_eq!(r.den(), 2);

    let r = parse_rational("-1.25").unwrap();
    assert_eq!(r.num(), -5);
    assert_eq!(r.den(), 4);
}

#[test]
fn rational_parse_fraction() {
    let r = parse_rational("3/4").unwrap();
    assert_eq!(r.num(), 3);
    assert_eq!(r.den(), 4);

    let r = parse_rational("-3/4").unwrap();
    assert_eq!(r.num(), -3);
    assert_eq!(r.den(), 4);
}

#[test]
fn rational_reduction() {
    let r = Rational::new(6, 8).unwrap();
    assert_eq!(r.num(), 3);
    assert_eq!(r.den(), 4);
}

#[test]
fn rational_equality_after_reduction() {
    let a = Rational::new(1, 2).unwrap();
    let b = Rational::new(2, 4).unwrap();
    assert_eq!(a, b);
}

#[test]
fn rational_ordering() {
    let a = Rational::new(1, 3).unwrap();
    let b = Rational::new(1, 2).unwrap();
    assert!(a < b);
    assert!(b > a);
}

#[test]
fn rational_large_arithmetic() {
    // i128 intermediates should handle large denominators
    let a = Rational::new(1, 1_000_000_000).unwrap();
    let b = Rational::new(1, 1_000_000_000).unwrap();
    let c = a + b;
    assert_eq!(c.num(), 1);
    assert_eq!(c.den(), 500_000_000);
}

#[test]
fn rational_chained_arithmetic() {
    // (1/2 + 1/3) * (1/4) = (5/6) * (1/4) = 5/24
    let result = (parse_rational("1/2").unwrap() + parse_rational("1/3").unwrap())
        * parse_rational("1/4").unwrap();
    assert_eq!(result.num(), 5);
    assert_eq!(result.den(), 24);
}

#[test]
fn rational_to_f64_accuracy() {
    let r = Rational::new(22, 7).unwrap();
    assert!((r.to_f64() - 22.0 / 7.0).abs() < 1e-15);
}

#[test]
fn rational_reciprocal() {
    let a = Rational::new(3, 4).unwrap();
    let r = a.recip().unwrap();
    assert_eq!(r.num(), 4);
    assert_eq!(r.den(), 3);
}

#[test]
fn rational_abs_and_neg() {
    let a = Rational::new(-3, 4).unwrap();
    assert_eq!(a.abs(), Rational::new(3, 4).unwrap());
    assert_eq!(-a, Rational::new(3, 4).unwrap());
}

#[test]
fn rational_display() {
    assert_eq!(Rational::from_int(5).to_string(), "5");
    assert_eq!(Rational::new(3, 4).unwrap().to_string(), "3/4");
    assert_eq!(Rational::new(-3, 4).unwrap().to_string(), "-3/4");
}

#[test]
fn rational_zero_den_errors() {
    assert!(Rational::new(1, 0).is_err());
    assert!(parse_rational("1/0").is_err());
}

// =========================================================================
// Cross-module pipelines
// =========================================================================

#[test]
fn rational_to_f64_then_eval() {
    // Convert rational to f64 and use in expression evaluation
    let r = Rational::new(1, 4).unwrap();
    let val = r.to_f64();
    let mut ctx = Context::standard();
    ctx.set("x", val);
    let expr = Parser::parse("x * 4 + 1").unwrap();
    let result = eval(&expr, &ctx).unwrap();
    assert!(close(result, 2.0, 1e-15));
}

#[test]
fn laurent_then_solve_near_pole() {
    // Compute Laurent series of 1/(x-1), then use solver to find root of f(x) - 0.5
    // 1/(x-1) = 0.5 → x = 3
    let ls = laurent::laurent_series_str("1/(x-1)", "x", 1.0, 1, 3).unwrap();
    // The series evaluated at x=3 should be close to 0.5
    let val = ls.eval(3.0);
    assert!(close(val, 0.5, 1e-4));
}

#[test]
fn taylor_then_laurent_consistency() {
    // For a function with no pole, Taylor and Laurent(pole_order=0) should agree
    let taylor_series = taylor::taylor_series_str("cos(x)", "x", 0.0, 5).unwrap();
    let laurent_series = laurent::laurent_series_str("cos(x)", "x", 0.0, 0, 5).unwrap();
    // Compare first few coefficients
    assert!(close(laurent_series.coeff(0), 1.0, 1e-6));
    // Taylor c_0 should also be 1
    let mut ctx = Context::standard();
    ctx.set("x", 0.0);
    let t0 = eval(&taylor_series, &ctx).unwrap();
    assert!(close(t0, 1.0, 1e-10));
    // Both should agree at x=0.5
    ctx.set("x", 0.5);
    let tv = eval(&taylor_series, &ctx).unwrap();
    let lv = laurent_series.eval(0.5);
    assert!(close(tv, lv, 1e-3));
}

#[test]
fn rational_arithmetic_repl_all_ops() {
    // Test all four operators via REPL
    let ctx = mathr::eval::Context::standard();

    let r1 = mathr::repl::dispatch_str("rat 1/2 + 1/4", ctx.clone()).unwrap().unwrap();
    assert!(r1.contains("3/4"), "1/2 + 1/4 should give 3/4: {}", r1);

    let r2 = mathr::repl::dispatch_str("rat 1/2 - 1/4", ctx.clone()).unwrap().unwrap();
    assert!(r2.contains("1/4"), "1/2 - 1/4 should give 1/4: {}", r2);

    let r3 = mathr::repl::dispatch_str("rat 2/3 * 3/4", ctx.clone()).unwrap().unwrap();
    assert!(r3.contains("1/2"), "2/3 * 3/4 should give 1/2: {}", r3);

    let r4 = mathr::repl::dispatch_str("rat 2/3 / 4/5", ctx.clone()).unwrap().unwrap();
    assert!(r4.contains("5/6"), "2/3 / 4/5 should give 5/6: {}", r4);
}

#[test]
fn rational_decimal_repl() {
    let ctx = mathr::eval::Context::standard();
    let r = mathr::repl::dispatch_str("rat 0.5 + 0.25", ctx).unwrap().unwrap();
    assert!(r.contains("3/4"), "0.5 + 0.25 should give 3/4: {}", r);
}

#[test]
fn laurent_repl_nonzero_center() {
    let ctx = mathr::eval::Context::standard();
    let r = mathr::repl::dispatch_str("laurent 1/(x-2) 2 1 3", ctx).unwrap().unwrap();
    assert!(r.contains("1/(x - 2)"), "output should contain 1/(x - 2): {}", r);
}

// =========================================================================
// Notebook: .mnb file format, cell evaluation, save/load
// =========================================================================

#[test]
fn notebook_create_eval_cell() {
    let mut nb = Notebook::new();
    let id = nb.add_cell("sin(pi/4)");
    nb.eval_cell(id, &mut Context::standard()).unwrap();
    assert!(nb.cells[id].output.contains("0.707"), "output: {}", nb.cells[id].output);
}

#[test]
fn notebook_eval_tex_cell() {
    let mut nb = Notebook::new();
    let id = nb.add_cell(r"\frac{1}{2} + \frac{3}{4}");
    nb.eval_cell(id, &mut Context::standard()).unwrap();
    assert!(nb.cells[id].output.contains("1.25"), "output: {}", nb.cells[id].output);
}

#[test]
fn notebook_eval_all_cells() {
    let mut nb = Notebook::new();
    nb.add_cell("1 + 2");
    nb.add_cell("3 * 4");
    nb.add_cell("sin(0)");
    let mut ctx = Context::standard();
    nb.eval_all(&mut ctx).unwrap();
    assert!(nb.cells[0].output.contains("3"));
    assert!(nb.cells[1].output.contains("12"));
    assert!(nb.cells[2].output.contains("0"));
}

#[test]
fn notebook_eval_diff_cell() {
    let mut nb = Notebook::new();
    let id = nb.add_cell("diff x^3");
    nb.eval_cell(id, &mut Context::standard()).unwrap();
    assert!(nb.cells[id].output.contains("3") && nb.cells[id].output.contains("x"),
        "output should contain derivative: {}", nb.cells[id].output);
}

#[test]
fn notebook_eval_solve_cell() {
    let mut nb = Notebook::new();
    let id = nb.add_cell("solve x^2 - 4");
    nb.eval_cell(id, &mut Context::standard()).unwrap();
    assert!(nb.cells[id].output.contains("2") || nb.cells[id].output.contains("root"),
        "output should contain root: {}", nb.cells[id].output);
}

#[test]
fn notebook_json_roundtrip() {
    let mut nb = Notebook::new();
    nb.add_cell("sin(pi/4)");
    nb.add_cell(r"\frac{1}{2}");
    nb.cells[0].output = "0.707...".to_string();
    nb.cells[1].output = "0.5".to_string();
    let json = nb.to_json();
    let nb2 = mathr::notebook::parse_notebook_json(&json).unwrap();
    assert_eq!(nb2.cells.len(), 2);
    assert_eq!(nb2.cells[0].input, "sin(pi/4)");
    assert_eq!(nb2.cells[0].output, "0.707...");
    assert_eq!(nb2.cells[1].input, r"\frac{1}{2}");
    assert_eq!(nb2.cells[1].output, "0.5");
}

#[test]
fn notebook_save_load_file() {
    let path = std::env::temp_dir().join("mathr_integration_test.mnb");
    let mut nb = Notebook::new();
    nb.add_cell("1 + 2");
    nb.add_cell("sin(pi/4)");
    nb.cells[0].output = "3".to_string();
    nb.save(&path).unwrap();
    let nb2 = Notebook::load(&path).unwrap();
    assert_eq!(nb2.cells.len(), 2);
    assert_eq!(nb2.cells[0].input, "1 + 2");
    assert_eq!(nb2.cells[0].output, "3");
    let _ = std::fs::remove_file(&path);
}

#[test]
fn notebook_set_input_clears_output() {
    let mut nb = Notebook::new();
    let id = nb.add_cell("1 + 1");
    nb.eval_cell(id, &mut Context::standard()).unwrap();
    assert!(!nb.cells[id].output.is_empty());
    nb.set_input(id, "2 + 2").unwrap();
    assert_eq!(nb.cells[id].input, "2 + 2");
    assert!(nb.cells[id].output.is_empty());
}

#[test]
fn notebook_remove_cell_reindexes() {
    let mut nb = Notebook::new();
    nb.add_cell("a");
    nb.add_cell("b");
    nb.add_cell("c");
    nb.remove_cell(1).unwrap();
    assert_eq!(nb.cells.len(), 2);
    assert_eq!(nb.cells[0].id, 0);
    assert_eq!(nb.cells[1].id, 1);
    assert_eq!(nb.cells[1].input, "c");
}

#[test]
fn notebook_parse_empty_cells() {
    let json = r#"{"cells": []}"#;
    let nb = mathr::notebook::parse_notebook_json(json).unwrap();
    assert_eq!(nb.cells.len(), 0);
}

#[test]
fn notebook_parse_bad_json_errors() {
    assert!(mathr::notebook::parse_notebook_json(r#"{"foo": "bar"}"#).is_err());
}

#[test]
fn notebook_json_escape_special_chars() {
    let mut nb = Notebook::new();
    nb.add_cell("a\nb\tc");
    let json = nb.to_json();
    assert!(json.contains("\\n"));
    assert!(json.contains("\\t"));
    let nb2 = mathr::notebook::parse_notebook_json(&json).unwrap();
    assert_eq!(nb2.cells[0].input, "a\nb\tc");
}

#[test]
fn notebook_load_example_file() {
    let path = std::path::Path::new(env!("CARGO_MANIFEST_DIR")).join("examples/notebooks/demo.mnb");
    let nb = Notebook::load(&path).unwrap();
    assert_eq!(nb.cells.len(), 4);
    assert_eq!(nb.cells[0].input, "sin(pi/4)");
    assert_eq!(nb.cells[1].input, r"\frac{1}{2} + \frac{3}{4}");
    assert_eq!(nb.cells[2].input, "diff x^3");
}

#[test]
fn notebook_eval_loaded_example() {
    let path = std::path::Path::new(env!("CARGO_MANIFEST_DIR")).join("examples/notebooks/demo.mnb");
    let mut nb = Notebook::load(&path).unwrap();
    let mut ctx = Context::standard();
    nb.eval_all(&mut ctx).unwrap();
    assert!(nb.cells[0].output.contains("0.707"));
    assert!(nb.cells[1].output.contains("1.25"));
    assert!(nb.cells[2].output.contains("3"));
}

// =========================================================================
// Notebook: cell types, reordering, shared context
// =========================================================================

#[test]
fn notebook_cell_types_and_reordering() {
    use mathr::notebook::{CellType, Notebook};
    let mut nb = Notebook::new();
    nb.add_cell("1 + 2");
    nb.add_cell_with_type("# Notes", CellType::Text);
    nb.add_cell("sin(0)");

    // Move the text cell up
    nb.move_cell_up(1).unwrap();
    assert_eq!(nb.cells[0].cell_type, CellType::Text);
    assert_eq!(nb.cells[1].cell_type, CellType::Math);

    // Duplicate
    let new_id = nb.duplicate_cell(0).unwrap();
    assert_eq!(nb.cells.len(), 4);
    assert_eq!(nb.cells[new_id].cell_type, CellType::Text);
}

#[test]
fn notebook_shared_context_across_cells() {
    use mathr::notebook::Notebook;
    let mut nb = Notebook::new();
    nb.add_cell("let x = 7");
    nb.add_cell("x * 3");
    let mut ctx = Context::standard();
    nb.eval_all(&mut ctx).unwrap();
    // Second cell should see x=7 from the first cell
    assert!(nb.cells[1].output.contains("21"), "output: {}", nb.cells[1].output);
}

#[test]
fn notebook_text_cell_skipped_in_eval() {
    use mathr::notebook::{CellType, Notebook};
    let mut nb = Notebook::new();
    nb.add_cell_with_type("This is a note", CellType::Text);
    nb.add_cell("1 + 1");
    let mut ctx = Context::standard();
    nb.eval_all(&mut ctx).unwrap();
    // Text cell output is the input itself
    assert_eq!(nb.cells[0].output, "This is a note");
    // Math cell still evaluates
    assert!(nb.cells[1].output.contains("2"));
}

#[test]
fn notebook_mathml_cell_eval() {
    use mathr::notebook::{CellType, Notebook};
    let mut nb = Notebook::new();
    // MathML cell: <mfrac><mn>1</mn><mn>2</mn></mfrac> → 1/2 → 0.5
    nb.add_cell_with_type(
        "<mfrac><mn>1</mn><mn>2</mn></mfrac>",
        CellType::MathML,
    );
    let mut ctx = Context::standard();
    nb.eval_all(&mut ctx).unwrap();
    assert!(nb.cells[0].output.contains("0.5"), "output: {}", nb.cells[0].output);
}

#[test]
fn notebook_mathml_cell_pow() {
    use mathr::notebook::{CellType, Notebook};
    let mut nb = Notebook::new();
    // <msup><mi>x</mi><mn>2</mn></msup> with x=3 → 9
    nb.add_cell_with_type("let x = 3", CellType::Math);
    nb.add_cell_with_type("<msup><mi>x</mi><mn>2</mn></msup>", CellType::MathML);
    let mut ctx = Context::standard();
    nb.eval_all(&mut ctx).unwrap();
    assert!(nb.cells[1].output.contains("9"), "output: {}", nb.cells[1].output);
}

#[test]
fn notebook_mathml_cell_type_roundtrip() {
    use mathr::notebook::{CellType, Notebook, parse_notebook_json};
    let mut nb = Notebook::new();
    nb.add_cell_with_type("<mn>42</mn>", CellType::MathML);
    let json = nb.to_json();
    let parsed = parse_notebook_json(&json).unwrap();
    assert_eq!(parsed.cells[0].cell_type, CellType::MathML);
}

#[test]
fn notebook_json_roundtrip_with_cell_types() {
    use mathr::notebook::{CellType, Notebook, parse_notebook_json};
    let mut nb = Notebook::new();
    nb.add_cell("sin(pi/4)");
    nb.add_cell_with_type("# Markdown note", CellType::Text);
    let json = nb.to_json();
    let nb2 = parse_notebook_json(&json).unwrap();
    assert_eq!(nb2.cells[0].cell_type, CellType::Math);
    assert_eq!(nb2.cells[1].cell_type, CellType::Text);
}

// =========================================================================
// Fast math: Chebyshev-based approximations
// =========================================================================

#[test]
fn fastmath_chebyshev_approx_custom_function() {
    use mathr::fastmath::ChebyshevApprox;
    let approx = ChebyshevApprox::new(|x| x * x + 1.0, -5.0, 5.0, 8);
    for x in [-5.0, -2.0, 0.0, 1.5, 5.0] {
        assert!(close(approx.eval(x), x * x + 1.0, 1e-10), "at x={}", x);
    }
}

#[test]
fn fastmath_sin_cos_accuracy() {
    use mathr::fastmath::{fast_cos, fast_sin};
    for x in (-500..=500).step_by(7) {
        let x = x as f64 * 0.01;
        assert!((fast_sin(x) - x.sin()).abs() < 1e-13, "sin({})", x);
        assert!((fast_cos(x) - x.cos()).abs() < 1e-13, "cos({})", x);
    }
}

#[test]
fn fastmath_exp_log_accuracy() {
    use mathr::fastmath::{fast_exp, fast_log};
    for x in (-50..=50).step_by(3) {
        let x = x as f64 * 0.1;
        let rel = (fast_exp(x) - x.exp()).abs() / x.exp().abs();
        assert!(rel < 1e-13, "exp({})", x);
    }
    for i in 1..=100_000 {
        let x = i as f64 * 0.001;
        assert!((fast_log(x) - x.ln()).abs() < 1e-13, "log({})", x);
    }
}

#[test]
fn fastmath_repl_dispatch() {
    let ctx = mathr::eval::Context::standard();
    let result = mathr::repl::dispatch_str("fast sin 1.5", ctx).unwrap().unwrap();
    assert!(result.contains("fast sin(1.5)"));
    assert!(result.contains("err:"));

    let ctx = mathr::eval::Context::standard();
    let result = mathr::repl::dispatch_str("fast exp 2.0", ctx).unwrap().unwrap();
    assert!(result.contains("fast exp(2)"));
}

// =========================================================================
// Inline plot rendering (in-memory PNG bytes)
// =========================================================================

#[test]
fn plot_function_to_bytes_produces_valid_png() {
    use mathr::plot::plot_function_to_bytes;
    let expr = Parser::parse("sin(x)").unwrap();
    let bytes = plot_function_to_bytes(&expr, "x", 0.0, std::f64::consts::PI, 100, "y = sin(x)").unwrap();
    // PNG magic bytes
    assert_eq!(&bytes[..8], &[137, 80, 78, 71, 13, 10, 26, 10]);
    assert!(bytes.len() > 1000, "PNG should be substantial: {} bytes", bytes.len());
}

#[test]
fn plot_scatter_to_bytes_produces_valid_png() {
    use mathr::plot::plot_scatter_to_bytes;
    let points = vec![(0.0, 0.0), (1.0, 1.0), (2.0, 4.0), (3.0, 9.0)];
    let bytes = plot_scatter_to_bytes(&points, "quadratic", "x", "y").unwrap();
    assert_eq!(&bytes[..8], &[137, 80, 78, 71, 13, 10, 26, 10]);
}

#[test]
fn plot_function_to_bytes_bad_range_errors() {
    use mathr::plot::plot_function_to_bytes;
    let expr = Parser::parse("x").unwrap();
    let result = plot_function_to_bytes(&expr, "x", 5.0, 1.0, 100, "y = x");
    assert!(result.is_err());
}

#[test]
fn notebook_features_example_loads() {
    let path = std::path::Path::new(env!("CARGO_MANIFEST_DIR"))
        .join("examples/notebooks/notebook_features.mnb");
    let nb = Notebook::load(&path).unwrap();
    assert!(nb.cells.len() > 10, "should have many cells");
    // First cell should be a text cell
    assert_eq!(nb.cells[0].cell_type, mathr::notebook::CellType::Text);
    // Should have math cells too
    assert!(nb.cells.iter().any(|c| c.cell_type == mathr::notebook::CellType::Math));
}

// =========================================================================
// Big integer arithmetic (arbitrary precision)
// =========================================================================

#[test]
fn bigint_repl_prime() {
    let ctx = mathr::eval::Context::standard();
    let result = mathr::repl::dispatch_str("big prime 1000000007", ctx).unwrap().unwrap();
    assert!(result.contains("prime"), "result: {}", result);

    let ctx = mathr::eval::Context::standard();
    let result = mathr::repl::dispatch_str("big prime 1000000008", ctx).unwrap().unwrap();
    assert!(result.contains("composite"), "result: {}", result);
}

#[test]
fn bigint_repl_factor() {
    let ctx = mathr::eval::Context::standard();
    let result = mathr::repl::dispatch_str("big factor 360", ctx).unwrap().unwrap();
    assert!(result.contains("2^3"), "result: {}", result);
    assert!(result.contains("3^2"), "result: {}", result);
    assert!(result.contains("·"), "result: {}", result); // middle dot notation
}

#[test]
fn bigint_repl_factorial() {
    let ctx = mathr::eval::Context::standard();
    let result = mathr::repl::dispatch_str("fact 25", ctx).unwrap().unwrap();
    // 25! = 15511210043330985984000000 — auto-upgraded to BigInt
    assert!(result.contains("15511210043330985984000000"), "result: {}", result);
}

#[test]
fn bigint_repl_fibonacci() {
    let ctx = mathr::eval::Context::standard();
    let result = mathr::repl::dispatch_str("fib 100", ctx).unwrap().unwrap();
    // F_100 = 354224848179261915075 — auto-upgraded to BigInt
    assert!(result.contains("354224848179261915075"), "result: {}", result);
}

#[test]
fn bigint_repl_binomial() {
    let ctx = mathr::eval::Context::standard();
    let result = mathr::repl::dispatch_str("binom 100 50", ctx).unwrap().unwrap();
    // C(100,50) = 100891344545564193334812497256 — auto-upgraded to BigInt
    assert!(result.contains("100891344545564193334812497256"), "result: {}", result);
}

#[test]
fn bigint_repl_gcd() {
    let ctx = mathr::eval::Context::standard();
    let result = mathr::repl::dispatch_str("big gcd 1234567890123456 987654321098765", ctx)
        .unwrap().unwrap();
    assert!(result.contains("gcd("), "result: {}", result);
}

#[test]
fn bigint_repl_modpow() {
    let ctx = mathr::eval::Context::standard();
    let result = mathr::repl::dispatch_str("big modpow 2 100 1000000007", ctx)
        .unwrap().unwrap();
    assert!(result.contains("≡"), "result: {}", result);
    assert!(result.contains("(mod"), "result: {}", result);
}

#[test]
fn bigint_repl_totient() {
    let ctx = mathr::eval::Context::standard();
    let result = mathr::repl::dispatch_str("big totient 360", ctx).unwrap().unwrap();
    assert!(result.contains("φ("), "result: {}", result);
    assert!(result.contains("96"), "result: {}", result); // φ(360) = 96
}

#[test]
fn bigint_factorize_large_semiprime() {
    use mathr::bigint;
    let n = num_bigint::BigInt::from(1000000007u64) * num_bigint::BigInt::from(1000000009u64);
    let factors = bigint::factorize(&n);
    assert_eq!(factors.len(), 2);
    assert_eq!(factors[0].0, num_bigint::BigInt::from(1000000007u64));
    assert_eq!(factors[1].0, num_bigint::BigInt::from(1000000009u64));
}

#[test]
fn bigint_is_prime_mersenne() {
    use mathr::bigint;
    use num_traits::One;
    // 2^61 - 1 is a Mersenne prime
    let m = num_bigint::BigInt::from(1u64 << 61) - num_bigint::BigInt::one();
    assert!(bigint::is_prime(&m, 20));
}

// =========================================================================
// Automatic differentiation (dual numbers)
// =========================================================================

#[test]
fn autodiff_repl_derivative() {
    let ctx = mathr::eval::Context::standard();
    let result = mathr::repl::dispatch_str("ad x^3 + 2*x^2 - x + 5 at x=2", ctx)
        .unwrap().unwrap();
    // f(2) = 8 + 8 - 2 + 5 = 19, f'(2) = 12 + 8 - 1 = 19
    assert!(result.contains("19"), "result: {}", result);
    assert!(result.contains("f'"), "result: {}", result);
}

#[test]
fn autodiff_repl_trig() {
    let ctx = mathr::eval::Context::standard();
    let result = mathr::repl::dispatch_str("ad sin(x) at x=0", ctx)
        .unwrap().unwrap();
    // f(0) = 0, f'(0) = cos(0) = 1
    assert!(result.contains("f(x) = 0"), "result: {}", result);
    assert!(result.contains("f'(x) = 1"), "result: {}", result);
}

#[test]
fn autodiff_repl_gradient() {
    let ctx = mathr::eval::Context::standard();
    let result = mathr::repl::dispatch_str("ad grad x^2 + y^3 with x=2,y=3", ctx)
        .unwrap().unwrap();
    // df/dx = 2x = 4, df/dy = 3y^2 = 27
    assert!(result.contains("4"), "result: {}", result);
    assert!(result.contains("27"), "result: {}", result);
}

#[test]
fn autodiff_repl_jacobian() {
    let ctx = mathr::eval::Context::standard();
    let result = mathr::repl::dispatch_str(
        "ad jacobian x^2 + y, x * y^2 with x=2,y=3",
        ctx,
    ).unwrap().unwrap();
    // J = [[4, 1], [9, 12]]
    assert!(result.contains("4"), "result: {}", result);
    assert!(result.contains("9"), "result: {}", result);
    assert!(result.contains("12"), "result: {}", result);
}

#[test]
fn autodiff_derivative_polynomial() {
    use mathr::autodiff;
    use mathr::parser::Parser;
    use mathr::eval::Context;
    let expr = Parser::parse("x^3 + 2*x^2 - x + 5").unwrap();
    let ctx = Context::standard();
    let d = autodiff::derivative(&expr, "x", 2.0, &ctx).unwrap();
    assert!((d.val - 19.0).abs() < 1e-10);
    assert!((d.deriv - 19.0).abs() < 1e-10);
}

#[test]
fn autodiff_gradient_multivariate() {
    use mathr::autodiff;
    use mathr::parser::Parser;
    use mathr::eval::Context;
    let expr = Parser::parse("x^2 + y^3").unwrap();
    let mut ctx = Context::standard();
    ctx.set("x", 2.0);
    ctx.set("y", 3.0);
    let grad = autodiff::gradient(&expr, &ctx).unwrap();
    let x_grad = grad.iter().find(|(n, _)| n == "x").unwrap().1;
    let y_grad = grad.iter().find(|(n, _)| n == "y").unwrap().1;
    assert!((x_grad - 4.0).abs() < 1e-10);
    assert!((y_grad - 27.0).abs() < 1e-10);
}

#[test]
fn autodiff_jacobian_matrix() {
    use mathr::autodiff;
    use mathr::parser::Parser;
    use mathr::eval::Context;
    let f1 = Parser::parse("x^2 + y").unwrap();
    let f2 = Parser::parse("x * y^2").unwrap();
    let mut ctx = Context::standard();
    ctx.set("x", 2.0);
    ctx.set("y", 3.0);
    let jac = autodiff::jacobian(&[f1, f2], &ctx).unwrap();
    assert!((jac[0][0] - 4.0).abs() < 1e-10); // 2x = 4
    assert!((jac[0][1] - 1.0).abs() < 1e-10);
    assert!((jac[1][0] - 9.0).abs() < 1e-10); // y^2 = 9
    assert!((jac[1][1] - 12.0).abs() < 1e-10); // 2xy = 12
}

// =========================================================================
// Auto-upgrade: REPL commands that auto-use BigInt on overflow
// =========================================================================

#[test]
fn autoupgrade_factorial() {
    // 25! overflows u64 — should auto-upgrade to BigInt
    let ctx = mathr::eval::Context::standard();
    let result = mathr::repl::dispatch_str("fact 25", ctx).unwrap().unwrap();
    assert_eq!(result, "15511210043330985984000000");
}

#[test]
fn autoupgrade_fibonacci() {
    // F(100) overflows u64 — should auto-upgrade to BigInt
    let ctx = mathr::eval::Context::standard();
    let result = mathr::repl::dispatch_str("fib 100", ctx).unwrap().unwrap();
    assert_eq!(result, "354224848179261915075");
}

#[test]
fn autoupgrade_binomial() {
    // C(100, 50) overflows u64 — should auto-upgrade to BigInt
    let ctx = mathr::eval::Context::standard();
    let result = mathr::repl::dispatch_str("binom 100 50", ctx).unwrap().unwrap();
    assert_eq!(result, "100891344545564193334812497256");
}

#[test]
fn autoupgrade_factorial_small_stays_u64() {
    // 10! fits in u64 — should use the fast u64 path
    let ctx = mathr::eval::Context::standard();
    let result = mathr::repl::dispatch_str("fact 10", ctx).unwrap().unwrap();
    assert_eq!(result, "3628800");
}

// =========================================================================
// Postfix factorial notation (n!)
// =========================================================================

#[test]
fn postfix_factorial_eval() {
    let ctx = mathr::eval::Context::standard();
    let result = mathr::repl::dispatch_str("5!", ctx).unwrap().unwrap();
    assert_eq!(result, "120");
}

#[test]
fn postfix_factorial_large() {
    // 20! = 2432902008176640000 — largest factorial exact in f64
    let ctx = mathr::eval::Context::standard();
    let result = mathr::repl::dispatch_str("20!", ctx).unwrap().unwrap();
    assert_eq!(result, "2432902008176640000");
}

#[test]
fn postfix_factorial_in_expression() {
    let ctx = mathr::eval::Context::standard();
    let result = mathr::repl::dispatch_str("2 * 3!", ctx).unwrap().unwrap();
    assert_eq!(result, "12");
}

#[test]
fn postfix_factorial_parens() {
    let ctx = mathr::eval::Context::standard();
    let result = mathr::repl::dispatch_str("(2+3)!", ctx).unwrap().unwrap();
    assert_eq!(result, "120");
}

// =========================================================================
// Standard math notation: |x|, a mod b, gcd(a,b), lcm(a,b), C(n,k), \binom
// =========================================================================

#[test]
fn abs_bars_eval() {
    let ctx = mathr::eval::Context::standard();
    assert_eq!(mathr::repl::dispatch_str("|-5|", ctx.clone()).unwrap().unwrap(), "5");
    assert_eq!(mathr::repl::dispatch_str("|3|", ctx.clone()).unwrap().unwrap(), "3");
    assert_eq!(mathr::repl::dispatch_str("|-3|", ctx).unwrap().unwrap(), "3");
}

#[test]
fn abs_bars_expression() {
    let ctx = mathr::eval::Context::standard();
    let result = mathr::repl::dispatch_str("|sin(pi)|", ctx).unwrap().unwrap();
    let v: f64 = result.parse().unwrap();
    assert!(v.abs() < 1e-10);
}

#[test]
fn infix_mod_eval() {
    let ctx = mathr::eval::Context::standard();
    assert_eq!(mathr::repl::dispatch_str("7 mod 3", ctx.clone()).unwrap().unwrap(), "1");
    assert_eq!(mathr::repl::dispatch_str("10 mod 4", ctx).unwrap().unwrap(), "2");
}

#[test]
fn gcd_lcm_as_functions_eval() {
    let ctx = mathr::eval::Context::standard();
    assert_eq!(mathr::repl::dispatch_str("gcd(12, 8)", ctx.clone()).unwrap().unwrap(), "4");
    assert_eq!(mathr::repl::dispatch_str("lcm(4, 6)", ctx).unwrap().unwrap(), "12");
}

#[test]
fn binomial_c_function_eval() {
    let ctx = mathr::eval::Context::standard();
    assert_eq!(mathr::repl::dispatch_str("C(5, 2)", ctx.clone()).unwrap().unwrap(), "10");
    assert_eq!(mathr::repl::dispatch_str("C(10, 3)", ctx).unwrap().unwrap(), "120");
}

#[test]
fn tex_binom_eval() {
    let ctx = mathr::eval::Context::standard();
    let result = mathr::repl::dispatch_str("\\binom{5}{2}", ctx).unwrap().unwrap();
    assert_eq!(result, "10");
}

#[test]
fn tex_gcd_lcm_eval() {
    let ctx = mathr::eval::Context::standard();
    let result = mathr::repl::dispatch_str("\\gcd(12, 8)", ctx.clone()).unwrap().unwrap();
    assert_eq!(result, "4");
    let result = mathr::repl::dispatch_str("\\lcm(4, 6)", ctx).unwrap().unwrap();
    assert_eq!(result, "12");
}

// =========================================================================
// MathML export and import
// =========================================================================

#[test]
fn mathml_export_repl() {
    let ctx = mathr::eval::Context::standard();
    let result = mathr::repl::dispatch_str("mathml x^2 + 1", ctx).unwrap().unwrap();
    assert!(result.contains("<math"), "result: {}", result);
    assert!(result.contains("<msup>"), "result: {}", result);
    assert!(result.contains("<mi>x</mi>"), "result: {}", result);
}

#[test]
fn mathml_export_frac() {
    let ctx = mathr::eval::Context::standard();
    let result = mathr::repl::dispatch_str("mathml 1/2", ctx).unwrap().unwrap();
    assert!(result.contains("<mfrac>"), "result: {}", result);
}

#[test]
fn mathml_export_sqrt() {
    let ctx = mathr::eval::Context::standard();
    let result = mathr::repl::dispatch_str("mathml sqrt(x)", ctx).unwrap().unwrap();
    assert!(result.contains("<msqrt>"), "result: {}", result);
}

#[test]
fn mathml_import_repl() {
    let ctx = mathr::eval::Context::standard();
    let ml = "<math xmlns=\"http://www.w3.org/1998/Math/MathML\"><mfrac><mn>1</mn><mn>2</mn></mfrac></math>";
    let result = mathr::repl::dispatch_str(&format!("mathml import {}", ml), ctx)
        .unwrap().unwrap();
    assert_eq!(result, "1/2");
}

#[test]
fn mathml_import_pow() {
    let ctx = mathr::eval::Context::standard();
    let ml = "<msup><mi>x</mi><mn>2</mn></msup>";
    let result = mathr::repl::dispatch_str(&format!("mathml import {}", ml), ctx)
        .unwrap().unwrap();
    assert_eq!(result, "x^2");
}

#[test]
fn mathml_roundtrip_api() {
    use mathr::mathml;
    use mathr::parser::Parser;
    let original = Parser::parse("2*x + sin(x)").unwrap();
    let ml = mathml::to_mathml(&original);
    let parsed = mathml::from_mathml(&ml).unwrap();
    assert_eq!(parsed, original);
}

// =========================================================================
// Expression serialization (S-expr, JSON, RPN)
// =========================================================================

#[test]
fn serialize_sexpr_export_repl() {
    let ctx = mathr::eval::Context::standard();
    let result = mathr::repl::dispatch_str("serialize sexpr 2*x + 1", ctx)
        .unwrap().unwrap();
    assert_eq!(result, "(add (mul (num 2) (var x)) (num 1))");
}

#[test]
fn serialize_sexpr_import_repl() {
    let ctx = mathr::eval::Context::standard();
    let result = mathr::repl::dispatch_str(
        "serialize sexpr import (add (mul (num 2) (var x)) (num 1))",
        ctx,
    )
    .unwrap().unwrap();
    assert_eq!(result, "2*x + 1");
}

#[test]
fn serialize_json_export_repl() {
    let ctx = mathr::eval::Context::standard();
    let result = mathr::repl::dispatch_str("serialize json x^2", ctx)
        .unwrap().unwrap();
    assert!(result.contains("\"t\":\"pow\""), "result: {}", result);
    assert!(result.contains("\"v\":\"x\""), "result: {}", result);
}

#[test]
fn serialize_json_import_repl() {
    let ctx = mathr::eval::Context::standard();
    let j = r#"{"t":"pow","a":{"t":"var","v":"x"},"b":{"t":"num","v":2}}"#;
    let result = mathr::repl::dispatch_str(
        &format!("serialize json import {}", j),
        ctx,
    )
    .unwrap().unwrap();
    assert_eq!(result, "x^2");
}

#[test]
fn serialize_rpn_export_repl() {
    let ctx = mathr::eval::Context::standard();
    let result = mathr::repl::dispatch_str("serialize rpn 2*x + 1", ctx)
        .unwrap().unwrap();
    assert_eq!(result, "2 x * 1 +");
}

#[test]
fn serialize_rpn_import_repl() {
    let ctx = mathr::eval::Context::standard();
    let result = mathr::repl::dispatch_str(
        "serialize rpn import 2 x * 1 +",
        ctx,
    )
    .unwrap().unwrap();
    assert_eq!(result, "2*x + 1");
}

#[test]
fn serialize_rpn_func_export() {
    let ctx = mathr::eval::Context::standard();
    let result = mathr::repl::dispatch_str("serialize rpn sin(x)", ctx)
        .unwrap().unwrap();
    assert_eq!(result, "x sin:1");
}

#[test]
fn serialize_roundtrip_api_all_formats() {
    use mathr::parser::Parser;
    use mathr::serialize::{from_json, from_rpn, from_sexpr, to_json, to_rpn, to_sexpr};
    let original = Parser::parse("2*x + sin(x) - 1").unwrap();
    let cases: [(String, fn(&str) -> mathr::error::Result<mathr::expr::Expr>); 3] = [
        (to_sexpr(&original), from_sexpr),
        (to_json(&original), from_json),
        (to_rpn(&original), from_rpn),
    ];
    for (text, parse) in cases {
        let parsed = parse(&text).unwrap();
        assert!(original.equals(&parsed), "round-trip failed for: {}", text);
    }
}

#[test]
fn serialize_bad_format_repl() {
    let ctx = mathr::eval::Context::standard();
    let result = mathr::repl::dispatch_str("serialize xml x", ctx);
    assert!(result.is_err());
}

// =========================================================================
// Interval arithmetic
// =========================================================================

// Complex evaluation
// =========================================================================

#[test]
fn cval_repl_euler_identity() {
    let ctx = mathr::eval::Context::standard();
    let result = mathr::repl::dispatch_str("cval exp(i*pi)", ctx).unwrap().unwrap();
    assert_eq!(result, "-1");
}

#[test]
fn cval_repl_arithmetic_display() {
    let ctx = mathr::eval::Context::standard();
    let result = mathr::repl::dispatch_str("cval (1 + 2i) * (3 - i)", ctx).unwrap().unwrap();
    assert_eq!(result, "5 + 5i");
}

#[test]
fn cval_repl_sqrt_negative() {
    let ctx = mathr::eval::Context::standard();
    let result = mathr::repl::dispatch_str("cval sqrt(-1)", ctx).unwrap().unwrap();
    assert_eq!(result, "i");
}

#[test]
fn cval_repl_variables() {
    let ctx = mathr::eval::Context::standard();
    let result = mathr::repl::dispatch_str("cval x*i - y with x=2, y=-1", ctx).unwrap().unwrap();
    assert_eq!(result, "1 + 2i");
}

#[test]
fn cval_repl_error_on_unknown_var() {
    let ctx = mathr::eval::Context::standard();
    assert!(mathr::repl::dispatch_str("cval i + nosuchvar", ctx).is_err());
}

// Symbolic quadratic solve
// =========================================================================

#[test]
fn qsolve_repl_integer_roots() {
    let ctx = mathr::eval::Context::standard();
    let result = mathr::repl::dispatch_str("qsolve x^2 - 4", ctx).unwrap().unwrap();
    assert_eq!(result, "x = 2, x = -2");
}

#[test]
fn qsolve_repl_rhs_form() {
    let ctx = mathr::eval::Context::standard();
    let result = mathr::repl::dispatch_str("qsolve x^2 = 2x + 3", ctx).unwrap().unwrap();
    assert_eq!(result, "x = 3, x = -1");
}

#[test]
fn qsolve_repl_complex_roots() {
    let ctx = mathr::eval::Context::standard();
    let result = mathr::repl::dispatch_str("qsolve x^2 - 2x + 5", ctx).unwrap().unwrap();
    assert_eq!(result, "x = 1 + 2*i, x = 1 - 2*i");
}

#[test]
fn qsolve_repl_linear() {
    let ctx = mathr::eval::Context::standard();
    let result = mathr::repl::dispatch_str("qsolve 2x - 6", ctx).unwrap().unwrap();
    assert_eq!(result, "x = 3");
}

#[test]
fn qsolve_repl_steps_have_discriminant() {
    let ctx = mathr::eval::Context::standard();
    let steps = mathr::repl::dispatch_steps("qsolve x^2 - 2x + 5", ctx).unwrap();
    assert!(steps.len() >= 4);
    assert!(steps.iter().any(|s| s.contains("discriminant")));
    assert!(steps.last().unwrap().contains("1 + 2*i"));
}

#[test]
fn qsolve_repl_rejects_non_polynomial() {
    let ctx = mathr::eval::Context::standard();
    assert!(mathr::repl::dispatch_str("qsolve sin(x) = 0", ctx).is_err());
}

// Summation / product
// =========================================================================

#[test]
fn sum_repl_arithmetic_series() {
    let ctx = mathr::eval::Context::standard();
    let result = mathr::repl::dispatch_str("sum k 1 100", ctx).unwrap().unwrap();
    assert_eq!(result, "5050");
}

#[test]
fn sum_repl_huge_range_exact() {
    let ctx = mathr::eval::Context::standard();
    let result = mathr::repl::dispatch_str("sum x 1 1000000", ctx).unwrap().unwrap();
    assert_eq!(result, "500000500000");
}

#[test]
fn sum_repl_exact_rational() {
    let ctx = mathr::eval::Context::standard();
    let result = mathr::repl::dispatch_str("sum 1/k 1 10", ctx).unwrap().unwrap();
    assert_eq!(result, "7381/2520");
}

#[test]
fn prod_repl_basic() {
    let ctx = mathr::eval::Context::standard();
    let result = mathr::repl::dispatch_str("prod x 1 5", ctx).unwrap().unwrap();
    assert_eq!(result, "120");
}

#[test]
fn sum_repl_rejects_huge_range() {
    let ctx = mathr::eval::Context::standard();
    assert!(mathr::repl::dispatch_str("sum x 1 2000000", ctx).is_err());
}

// Interval arithmetic
// =========================================================================

#[test]
fn interval_repl_polynomial() {
    let ctx = mathr::eval::Context::standard();
    let result = mathr::repl::dispatch_str(
        "interval x^2 + 1 with x=[-2,3]",
        ctx,
    )
    .unwrap().unwrap();
    assert_eq!(result, "[1, 10]");
}

#[test]
fn interval_repl_sin_full_range() {
    let ctx = mathr::eval::Context::standard();
    let result = mathr::repl::dispatch_str(
        "interval sin(x) with x=[0,6.283185307179586]",
        ctx,
    )
    .unwrap().unwrap();
    assert_eq!(result, "[-1, 1]");
}

#[test]
fn interval_repl_multivar() {
    let ctx = mathr::eval::Context::standard();
    let result = mathr::repl::dispatch_str(
        "interval x*y with x=[1,2],y=[3,4]",
        ctx,
    )
    .unwrap().unwrap();
    assert_eq!(result, "[3, 8]");
}

#[test]
fn interval_repl_dependency_problem() {
    // x - x over [1,2] → [-1, 1] (dependency problem, not [0,0])
    let ctx = mathr::eval::Context::standard();
    let result = mathr::repl::dispatch_str(
        "interval x - x with x=[1,2]",
        ctx,
    )
    .unwrap().unwrap();
    assert_eq!(result, "[-1, 1]");
}

#[test]
fn interval_repl_div_by_zero() {
    let ctx = mathr::eval::Context::standard();
    let result = mathr::repl::dispatch_str(
        "interval 1/x with x=[-1,1]",
        ctx,
    )
    .unwrap().unwrap();
    assert_eq!(result, "[-∞, ∞]");
}

#[test]
fn interval_repl_exp() {
    let ctx = mathr::eval::Context::standard();
    let result = mathr::repl::dispatch_str(
        "interval exp(x) with x=[0,1]",
        ctx,
    )
    .unwrap().unwrap();
    // exp(0)=1, exp(1)=e≈2.71828...
    assert!(result.starts_with("[1, 2.718"));
}

#[test]
fn interval_repl_odd_power() {
    // x^3 over [-2,1] → [-8, 1] (odd, monotonic)
    let ctx = mathr::eval::Context::standard();
    let result = mathr::repl::dispatch_str(
        "interval x^3 with x=[-2,1]",
        ctx,
    )
    .unwrap().unwrap();
    assert_eq!(result, "[-8, 1]");
}

#[test]
fn interval_repl_missing_with_errors() {
    let ctx = mathr::eval::Context::standard();
    let result = mathr::repl::dispatch_str("interval x^2", ctx);
    assert!(result.is_err());
}

#[test]
fn interval_repl_unknown_var_errors() {
    let ctx = mathr::eval::Context::standard();
    let result = mathr::repl::dispatch_str("interval y + 1 with x=[1,2]", ctx);
    assert!(result.is_err());
}

#[test]
fn interval_api_eval() {
    use mathr::interval::{eval_interval, Interval};
    use mathr::parser::Parser;
    use std::collections::HashMap;
    let e = Parser::parse("x^2 + 2*x + 1").unwrap();
    let mut vars = HashMap::new();
    vars.insert("x".to_string(), Interval::new(-1.0, 1.0));
    let r = eval_interval(&e, &vars).unwrap();
    // (x+1)^2 over [-1, 1] → [0, 4]
    assert!(r.contains(0.0));
    assert!(r.contains(4.0));
}

#[test]
fn interval_api_empty_on_negative_sqrt() {
    use mathr::interval::{eval_interval, Interval};
    use mathr::parser::Parser;
    use std::collections::HashMap;
    let e = Parser::parse("sqrt(x)").unwrap();
    let mut vars = HashMap::new();
    vars.insert("x".to_string(), Interval::new(-4.0, -1.0));
    let r = eval_interval(&e, &vars).unwrap();
    assert!(r.is_empty);
}




// =========================================================================
// Arbitrary-precision decimal arithmetic (bigdec)
// =========================================================================

#[test]
fn bigdec_pi_and_e_reference_digits() {
    let pi = mathr::bigdec::pi(50).unwrap();
    assert_eq!(
        pi.to_string(),
        "3.1415926535897932384626433832795028841971693993751"
    );
    let e = mathr::bigdec::e(30).unwrap();
    assert!(e.to_string().starts_with("2.71828182845904523536028747135"));
}

#[test]
fn bigdec_sqrt_exp_ln_round_trip() {
    use mathr::bigdec::BigDecimal;
    let two = BigDecimal::from(2);
    let s = mathr::bigdec::sqrt(&two, 30).unwrap();
    assert!(s.to_string().starts_with("1.4142135623730950488016887242"));
    let five = BigDecimal::from(5);
    let l = mathr::bigdec::ln(&five, 30).unwrap();
    let back = mathr::bigdec::exp(&l, 30).unwrap();
    let diff = (mathr::bigdec::to_f64(&back) - 5.0).abs();
    assert!(diff < 1e-25);
}

#[test]
fn bigdec_trig_matches_f64() {
    use mathr::bigdec::BigDecimal;
    let one = BigDecimal::from(1);
    let s = mathr::bigdec::sin(&one, 30).unwrap();
    assert!((mathr::bigdec::to_f64(&s) - 1.0f64.sin()).abs() < 1e-15);
    let t = mathr::bigdec::tanh(&one, 30).unwrap();
    assert!((mathr::bigdec::to_f64(&t) - 1.0f64.tanh()).abs() < 1e-15);
}

#[test]
fn bigdec_eval_decimal_with_vars() {
    use mathr::bigdec::{eval_decimal_rounded, parse as dec_parse};
    use std::collections::HashMap;
    let e = mathr::parser::Parser::parse("x^2 + 1/x").unwrap();
    let mut vars = HashMap::new();
    vars.insert("x".to_string(), dec_parse("1.5").unwrap());
    let r = eval_decimal_rounded(&e, &vars, 20).unwrap();
    // 2.25 + 2/3 = 2.9166666666666666667 (20 sig digits)
    assert!(r.to_string().starts_with("2.9166666666666666667"));
}

#[test]
fn bigdec_repl_dec_command() {
    let ctx = mathr::eval::Context::standard();
    let result = mathr::repl::dispatch_str("dec sqrt(2) prec 30", ctx).unwrap().unwrap();
    assert!(result.starts_with("1.4142135623730950488016887242"));

    let ctx = mathr::eval::Context::standard();
    let result = mathr::repl::dispatch_str("dec 1/3 prec 10", ctx).unwrap().unwrap();
    assert_eq!(result, "0.3333333333");

    let ctx = mathr::eval::Context::standard();
    let result = mathr::repl::dispatch_str("dec pi prec 50", ctx).unwrap().unwrap();
    assert!(result.starts_with("3.1415926535897932384626433832795028841971693993751"));
}

#[test]
fn bigdec_repl_dec_with_assignments() {
    let ctx = mathr::eval::Context::standard();
    let result = mathr::repl::dispatch_str("dec x*2 + 1 with x=1.5, prec 10", ctx)
        .unwrap()
        .unwrap();
    assert_eq!(result, "4.000000000");
}

#[test]
fn bigdec_repl_dec_domain_errors() {
    let ctx = mathr::eval::Context::standard();
    assert!(mathr::repl::dispatch_str("dec ln(-1) prec 10", ctx).is_err());
    let ctx = mathr::eval::Context::standard();
    assert!(mathr::repl::dispatch_str("dec 1 prec 0", ctx).is_err());
}

// =========================================================================
// Symbolic limits (limit)
// =========================================================================

#[test]
fn limit_pipeline_lhopital_and_probe() {
    use mathr::limit::{limit, LimitValue};
    use mathr::parser::Parser;

    // Removable singularity resolved by L'Hôpital
    let e = Parser::parse("(x^2 - 1)/(x - 1)").unwrap();
    assert_eq!(limit(&e, "x", 1.0).unwrap(), LimitValue::Finite(2.0));

    // sin(x)/x at 0 via L'Hôpital
    let e = Parser::parse("sin(x)/x").unwrap();
    assert_eq!(limit(&e, "x", 0.0).unwrap(), LimitValue::Finite(1.0));

    // Pole at a finite point
    let e = Parser::parse("1/x^2").unwrap();
    assert_eq!(limit(&e, "x", 0.0).unwrap(), LimitValue::PosInfinity);

    // Opposite-side divergence
    let e = Parser::parse("1/x").unwrap();
    assert_eq!(limit(&e, "x", 0.0).unwrap(), LimitValue::DoesNotExist);

    // Limits at infinity
    let e = Parser::parse("(2*x + 1)/(x + 5)").unwrap();
    assert_eq!(limit(&e, "x", f64::INFINITY).unwrap(), LimitValue::Finite(2.0));
    let e = Parser::parse("x/exp(x)").unwrap();
    assert_eq!(limit(&e, "x", f64::INFINITY).unwrap(), LimitValue::Finite(0.0));
}

#[test]
fn limit_repl_dispatch() {
    let ctx = mathr::eval::Context::standard();
    let r = mathr::repl::dispatch_str("limit sin(x)/x 0", ctx).unwrap().unwrap();
    assert!(r.contains("= 1"), "got: {}", r);

    let ctx = mathr::eval::Context::standard();
    let r = mathr::repl::dispatch_str("limit (x^2-1)/(x-1) x 1", ctx).unwrap().unwrap();
    assert!(r.contains("= 2"), "got: {}", r);

    let ctx = mathr::eval::Context::standard();
    let r = mathr::repl::dispatch_str("limit exp(-x) x inf", ctx).unwrap().unwrap();
    assert!(r.contains("= 0"), "got: {}", r);
}

#[test]
fn limit_steps_via_dispatch() {
    let steps = mathr::repl::dispatch_steps("limit sin(x)/x x 0", mathr::eval::Context::standard()).unwrap();
    assert!(steps[0].contains("limit of"));
    assert!(steps.last().unwrap().contains("limit = 1"));
}

// =========================================================================
// Polynomial expansion (poly)
// =========================================================================

#[test]
fn poly_expand_pipeline() {
    use mathr::poly::expand;
    let e = mathr::parser::Parser::parse("(x+1)^3").unwrap();
    let expanded = expand(&e);
    assert!(expanded.to_string().contains("x^3"));
    assert!(expanded.to_string().contains("3*x^2"));

    // Expanded polynomial evaluates identically to the source
    let mut ctx = mathr::eval::Context::standard();
    ctx.set("x", 2.5);
    let a = mathr::eval::eval(&e, &ctx).unwrap();
    let b = mathr::eval::eval(&expanded, &ctx).unwrap();
    assert!((a - b).abs() < 1e-10);
}

#[test]
fn poly_expand_multivariate_and_functions() {
    use mathr::poly::expand;
    let e = mathr::parser::Parser::parse("(x+y)*(x-y)").unwrap();
    assert!(expand(&e).to_string().contains("x^2 - y^2"));

    // Functions stay opaque but polynomial sides distribute
    let e = mathr::parser::Parser::parse("sin(x)*(x + 1)").unwrap();
    let out = expand(&e).to_string();
    assert!(out.contains("sin(x)"), "got: {}", out);
}

#[test]
fn poly_expand_repl_dispatch() {
    let ctx = mathr::eval::Context::standard();
    let r = mathr::repl::dispatch_str("expand (x+2)*(x+3)", ctx).unwrap().unwrap();
    assert!(r.contains("x^2 + 5*x + 6"), "got: {}", r);
}

// =========================================================================
// Partial fraction decomposition (apart)
// =========================================================================

#[test]
fn apart_numeric_equivalence_pipeline() {
    use mathr::apart::apart;
    use mathr::eval::{eval, Context};
    use mathr::parser::Parser;

    for (src, points) in [
        ("1/(x^2 - 1)", vec![0.3, 0.9, 2.2, -3.1]),
        ("1/(x*(x+1)^2)", vec![0.5, 2.0, -3.0]),
        ("(x^2 + 1)/(x - 1)", vec![0.5, 2.0, -1.5]),
        ("(3*x + 5)/((x^2 + x + 1)*(x - 2))", vec![0.5, 1.0, -2.0, 3.0]),
    ] {
        let e = Parser::parse(src).unwrap();
        let out = apart(&e, "x").unwrap();
        for &x in &points {
            let mut ctx = Context::standard();
            ctx.set("x", x);
            let a = eval(&e, &ctx).unwrap();
            let b = eval(&out, &ctx).unwrap();
            assert!(
                (a - b).abs() < 1e-6 * a.abs().max(1.0),
                "{}: mismatch at x={}: {} vs {}",
                src, x, a, b
            );
        }
    }
}

#[test]
fn apart_repl_dispatch() {
    let ctx = mathr::eval::Context::standard();
    let r = mathr::repl::dispatch_str("apart 1/(x*(x+1))", ctx).unwrap().unwrap();
    assert!(r.contains("1/x"), "got: {}", r);
    assert!(r.contains("1/(x + 1)"), "got: {}", r);

    let ctx = mathr::eval::Context::standard();
    let r = mathr::repl::dispatch_str("apart (x^2+1)/(x-1) x", ctx).unwrap().unwrap();
    assert!(r.contains("x + 1"), "got: {}", r);
}

#[test]
fn apart_steps_via_dispatch() {
    let steps = mathr::repl::dispatch_steps(
        "apart 1/(x*(x+1))",
        mathr::eval::Context::standard(),
    )
    .unwrap();
    assert!(steps[0].contains("apart"), "got: {:?}", steps);
    assert!(steps.last().unwrap().contains("1/x"), "got: {:?}", steps);
}

// =========================================================================
// Condition number & nullspace (matrix)
// =========================================================================

#[test]
fn matrix_condition_number_pipeline() {
    use mathr::matrix::Matrix;
    // Well-conditioned
    let d = Matrix::from_rows(&[vec![3.0, 0.0], vec![0.0, 1.0]]).unwrap();
    assert!((d.condition_number().unwrap() - 3.0).abs() < 1e-8);
    // Singular
    let s = Matrix::from_rows(&[vec![1.0, 1.0], vec![1.0, 1.0]]).unwrap();
    assert_eq!(s.condition_number().unwrap(), f64::INFINITY);
    // Hilbert: ill-conditioned
    let h = Matrix::hilbert(6);
    let k = h.condition_number().unwrap();
    assert!(k > 1e6 && k < 1e9);
}

#[test]
fn matrix_nullspace_pipeline() {
    use mathr::matrix::Matrix;
    // Full rank → empty
    let i = Matrix::identity(2);
    assert!(i.nullspace(0.0).unwrap().is_empty());
    // Rank-deficient → basis vectors satisfy A·v = 0
    let a = Matrix::from_rows(&[
        vec![1.0, 2.0, 0.0, 3.0],
        vec![2.0, 4.0, 0.0, 6.0],
        vec![1.0, 1.0, 1.0, 0.0],
    ]).unwrap();
    let basis = a.nullspace(0.0).unwrap();
    assert!(!basis.is_empty());
    for v in &basis {
        let av = a.mul_vec(v).unwrap();
        assert!(av.iter().all(|&x| x.abs() < 1e-8));
    }
}

#[test]
fn matrix_cond_null_repl_dispatch() {
    let ctx = mathr::eval::Context::standard();
    let r = mathr::repl::dispatch_str("cond 1 2 | 3 4", ctx).unwrap().unwrap();
    assert!(r.contains("cond ="), "got: {}", r);

    let ctx = mathr::eval::Context::standard();
    let r = mathr::repl::dispatch_str("cond 1 1 | 1 1", ctx).unwrap().unwrap();
    assert!(r.contains("inf"), "got: {}", r);

    let ctx = mathr::eval::Context::standard();
    let r = mathr::repl::dispatch_str("null 1 1 | 1 1", ctx).unwrap().unwrap();
    assert!(r.contains("dim 1"), "got: {}", r);

    let ctx = mathr::eval::Context::standard();
    let r = mathr::repl::dispatch_str("null 1 0 | 0 1", ctx).unwrap().unwrap();
    assert!(r.contains("full column rank"), "got: {}", r);
}

#[test]
fn dists_repl_dispatch() {
    // t distribution: t_0.975(10) = 2.228
    let ctx = mathr::eval::Context::standard();
    let r = mathr::repl::dispatch_str("dist t 2.228 10", ctx).unwrap().unwrap();
    assert!(r.contains("cdf = 0.97"), "got: {}", r);

    // chi-squared: k=2 closed form cdf = 1 - e^{-x/2}
    let ctx = mathr::eval::Context::standard();
    let r = mathr::repl::dispatch_str("dist chi2 2 2", ctx).unwrap().unwrap();
    assert!(r.contains("0.63212"), "got: {}", r);

    // binomial pmf: B(10, 0.5) at k=5 = 252/1024
    let ctx = mathr::eval::Context::standard();
    let r = mathr::repl::dispatch_str("dist binom 5 10 0.5", ctx).unwrap().unwrap();
    assert!(r.contains("0.24609"), "got: {}", r);

    // poisson pmf: P(2 | lambda=2) = 2 e^{-2}
    let ctx = mathr::eval::Context::standard();
    let r = mathr::repl::dispatch_str("dist poisson 2 2", ctx).unwrap().unwrap();
    assert!(r.contains("0.27067"), "got: {}", r);
}

#[test]
fn ttest_repl_dispatch() {
    // one-sample: data 1..5 vs mu = 0
    let ctx = mathr::eval::Context::standard();
    let r = mathr::repl::dispatch_str("ttest 0 1 2 3 4 5", ctx).unwrap().unwrap();
    assert!(r.contains("one-sample t-test"), "got: {}", r);
    assert!(r.contains("df=4"), "got: {}", r);

    // two-sample Welch: [1,2,3] vs [4,5,6]
    let ctx = mathr::eval::Context::standard();
    let r = mathr::repl::dispatch_str("ttest 1 2 3 | 4 5 6", ctx).unwrap().unwrap();
    assert!(r.contains("welch t-test"), "got: {}", r);

    // paired
    let ctx = mathr::eval::Context::standard();
    let r = mathr::repl::dispatch_str("ttest paired 8 7 6 9 10 | 6 7 5 7 8", ctx)
        .unwrap()
        .unwrap();
    assert!(r.contains("paired t-test"), "got: {}", r);
}

#[test]
fn chitest_anova_repl_dispatch() {
    // uniform chi-squared: 6 cells, total 88 -> chi2 = 2.0
    let ctx = mathr::eval::Context::standard();
    let r = mathr::repl::dispatch_str("chitest 16 18 16 14 12 12", ctx).unwrap().unwrap();
    assert!(r.contains("stat=2"), "got: {}", r);
    assert!(r.contains("df=5"), "got: {}", r);

    // anova: [1,2,3] [4,5,6] [7,8,9] -> F = 27, df = (2, 6)
    let ctx = mathr::eval::Context::standard();
    let r = mathr::repl::dispatch_str("anova 1 2 3 | 4 5 6 | 7 8 9", ctx).unwrap().unwrap();
    assert!(r.contains("one-way ANOVA"), "got: {}", r);
    assert!(r.contains("df=(2, 6)"), "got: {}", r);
}

#[test]
fn dists_library_pipeline() {
    // beta_inc CDF identity: binomial tail via incomplete beta
    let p = mathr::dists::binomial_cdf(5, 10, 0.3);
    let sum: f64 = (0..=5).map(|k| mathr::dists::binomial_pmf(k, 10, 0.3)).sum();
    assert!(close(p, sum, 1e-12));

    // F(1, nu) cdf at t^2 equals two-sided central t probability
    let p_f = mathr::dists::f_cdf(6.25, 1.0, 12.0);
    let p_t = 2.0 * mathr::dists::student_t_cdf(2.5, 12.0) - 1.0;
    assert!(close(p_f, p_t, 1e-12));

    // normal quantile round-trip
    let x = mathr::dists::normal_ppf(0.975);
    assert!(close(mathr::stats::normal_cdf(x, 0.0, 1.0), 0.975, 1e-10));
}

#[test]
fn qtile_repl_dispatch() {
    // standard normal critical value
    let ctx = mathr::eval::Context::standard();
    let r = mathr::repl::dispatch_str("qtile normal 0.975", ctx).unwrap().unwrap();
    assert!(r.contains("1.95996"), "got: {}", r);

    // t critical value: t_{0.975}(10) = 2.2281388...
    let ctx = mathr::eval::Context::standard();
    let r = mathr::repl::dispatch_str("qtile t 0.975 10", ctx).unwrap().unwrap();
    assert!(r.contains("2.22813"), "got: {}", r);

    // chi2 critical value: chi2_{0.95}(2) = -2 ln(0.05) = 5.99146...
    let ctx = mathr::eval::Context::standard();
    let r = mathr::repl::dispatch_str("qtile chi2 0.95 2", ctx).unwrap().unwrap();
    assert!(r.contains("5.99146"), "got: {}", r);

    // normal with mu/sigma
    let ctx = mathr::eval::Context::standard();
    let r = mathr::repl::dispatch_str("qtile normal 0.975 10 2", ctx).unwrap().unwrap();
    assert!(r.contains("13.9199"), "got: {}", r);

    // discrete dists rejected
    let ctx = mathr::eval::Context::standard();
    let r = mathr::repl::dispatch_str("qtile poisson 0.95 2", ctx);
    assert!(r.is_err());
}

#[test]
fn qtile_cdf_round_trip_pipeline() {
    // p-values <-> critical values agree across the dists module
    let mut ctx = mathr::eval::Context::standard();
    let stat = 2.5_f64;
    let r = mathr::repl::dispatch_str("dist t 2.5 12", ctx.clone()).unwrap().unwrap();
    let cdf: f64 = r
        .lines()
        .find(|l| l.starts_with("cdf"))
        .and_then(|l| l.split('=').nth(1))
        .and_then(|v| v.trim().parse().ok())
        .unwrap();
    let q = mathr::dists::student_t_ppf(cdf, 12.0);
    // cdf passes through the REPL's 12-significant-digit display snapping,
    // so the round trip can only recover ~1e-9, not full f64 precision
    assert!(close(q, stat, 1e-8), "q={q}");
    let _ = &mut ctx;
}

#[test]
fn fit_repl_dispatch() {
    // exact linear fit with explicit guesses
    let ctx = mathr::eval::Context::standard();
    let r = mathr::repl::dispatch_str("fit a*x + b with a=1, b=1 0 1 1 3 2 5", ctx)
        .unwrap()
        .unwrap();
    assert!(r.contains("a = 2"), "got: {}", r);
    assert!(r.contains("b = 1"), "got: {}", r);
    assert!(r.contains("converged"), "got: {}", r);

    // inferred parameters (all variables except the independent one), `in t`
    let ctx = mathr::eval::Context::standard();
    let r = mathr::repl::dispatch_str("fit a*t + b in t 0 1 1 3 2 5", ctx).unwrap().unwrap();
    assert!(r.contains("a = 2"), "got: {}", r);
    assert!(r.contains("b = 1"), "got: {}", r);

    // nonlinear model recovers rounded-off data parameters
    let ctx = mathr::eval::Context::standard();
    let r = mathr::repl::dispatch_str(
        "fit a*exp(-b*x) in x with a=3, b=0.5 0 3 1 1.4957 2 0.7458 3 0.3717 4 0.1852",
        ctx,
    )
    .unwrap()
    .unwrap();
    assert!(r.contains("a = 3.0000"), "got: {}", r);
    assert!(r.contains("b = 0.696"), "got: {}", r);
    assert!(r.contains("converged"), "got: {}", r);
}

#[test]
fn curve_fit_library_matches_regression() {
    // library-level: LM equals closed-form linear regression
    let data = [(0.0, 1.0), (1.0, 2.1), (2.0, 2.9), (3.0, 4.2)];
    let fit = mathr::curvefit::curve_fit(
        |x, p| Ok(p[0] * x + p[1]),
        &data,
        &[0.5, 0.5],
        &mathr::curvefit::LmOptions::default(),
    )
    .unwrap();
    let (slope, intercept) =
        mathr::stats::linear_regression(&[0.0, 1.0, 2.0, 3.0], &[1.0, 2.1, 2.9, 4.2]).unwrap();
    assert!(close(fit.params[0], slope, 1e-8));
    assert!(close(fit.params[1], intercept, 1e-8));
    assert!(fit.converged);
}

#[test]
fn nonparametric_repl_dispatch() {
    // Mann-Whitney: perfect separation of 8v8
    let ctx = mathr::eval::Context::standard();
    let r = mathr::repl::dispatch_str("mwu 1 2 3 4 5 6 7 8 | 9 10 11 12 13 14 15 16", ctx)
        .unwrap()
        .unwrap();
    assert!(r.contains("mann-whitney U"), "got: {}", r);
    assert!(r.contains("stat=0"), "got: {}", r);
    assert!(r.contains("p=0.00093"), "got: {}", r);

    // Wilcoxon paired: a = [8,7,6,9,10] vs b = [6,7,5,7,8] -> W+ = 10
    let ctx = mathr::eval::Context::standard();
    let r = mathr::repl::dispatch_str("wilcoxon 8 7 6 9 10 | 6 7 5 7 8", ctx).unwrap().unwrap();
    assert!(r.contains("wilcoxon signed-rank"), "got: {}", r);
    assert!(r.contains("stat=10"), "got: {}", r);

    // Kruskal-Wallis textbook triple -> H = 7.2
    let ctx = mathr::eval::Context::standard();
    let r = mathr::repl::dispatch_str("kw 2 4 3 | 5 6 7 | 8 10 9", ctx).unwrap().unwrap();
    assert!(r.contains("kruskal-wallis"), "got: {}", r);
    assert!(r.contains("stat=7.2"), "got: {}", r);
    assert!(r.contains("df=2"), "got: {}", r);
}

#[test]
fn spearman_bootstrap_dispatch() {
    // Spearman via library (monotone data -> 1)
    let rho = mathr::dists::spearman_corr(&[1.0, 2.0, 3.0, 4.0], &[2.0, 4.0, 6.0, 8.0]).unwrap();
    assert!(close(rho, 1.0, 1e-14));

    // bootstrap REPL: median CI of 1..9 is within [2, 8] and contains 5
    let ctx = mathr::eval::Context::standard();
    let r = mathr::repl::dispatch_str("boot median 1 2 3 4 5 6 7 8 9", ctx).unwrap().unwrap();
    assert!(r.contains("median = 5"), "got: {}", r);
    assert!(r.contains("95% ci = ["), "got: {}", r);
    assert!(r.contains("iters = 10000"), "got: {}", r);
}

#[test]
fn logit_repl_dispatch() {
    // balanced grouped design: closed form b0=-3ln2, b1=2ln2
    let ctx = mathr::eval::Context::standard();
    let r = mathr::repl::dispatch_str("logit 0 0 1 0 1 1 with 1 1 1 2 2 2", ctx)
        .unwrap()
        .unwrap();
    assert!(r.contains("intercept = -2.0794415"), "got: {}", r);
    assert!(r.contains("b1 = 1.38629435"), "got: {}", r);
    assert!(r.contains("converged"), "got: {}", r);

    // multivariate with '|' column separator: x1 has no effect
    let ctx = mathr::eval::Context::standard();
    let r = mathr::repl::dispatch_str(
        "logit 0 1 0 0 1 0 0 1 1 1 0 1 with 0 0 1 1 2 2 3 3 4 4 5 5 | 0 0 0 0 0 0 1 1 1 1 1 1",
        ctx,
    )
    .unwrap()
    .unwrap();
    assert!(r.contains("b1 = 0."), "got: {}", r);
    assert!(r.contains("b2 = 1.3862943603"), "got: {}", r);

    // library-level closed form + predict_proba round trip
    let x = [&[1.0, 1.0, 1.0, 2.0, 2.0, 2.0][..]];
    let y = [0.0, 0.0, 1.0, 0.0, 1.0, 1.0];
    let fit = mathr::logit::logistic_regression(&x, &y, &mathr::logit::LogitOptions::default())
        .unwrap();
    assert!(close(
        fit.coefficients[1],
        2.0 * std::f64::consts::LN_2,
        1e-8
    ));
    let p = mathr::logit::predict_proba(&fit.coefficients, &[2.0]).unwrap();
    assert!(close(p, 2.0 / 3.0, 1e-8));
}

#[test]
fn pchip_repl_dispatch() {
    // linear data reproduced exactly
    let ctx = mathr::eval::Context::standard();
    let r = mathr::repl::dispatch_str("pchip 0 1 1 3 2 5 0.5", ctx).unwrap().unwrap();
    assert!(r.contains("pchip(0.5) = 2"), "got: {}", r);

    // monotone data with a steep rise: value stays within the local knot range
    let ctx = mathr::eval::Context::standard();
    let r = mathr::repl::dispatch_str("pchip 0 0 1 0 2 0.3 3 4 4 4.2 2.5", ctx)
        .unwrap()
        .unwrap();
    assert!(r.contains("pchip(2.5) = 2.17"), "got: {}", r);

    // library-level: interpolates knots exactly, clamps outside
    let p = mathr::pchip::Pchip::new(&[0.0, 1.0, 2.0], &[1.0, 3.0, 5.0]).unwrap();
    assert!(close(p.eval(1.0), 3.0, 1e-14));
    assert_eq!(p.eval(9.0), 5.0);
}

#[test]
fn minimize_repl_dispatch() {
    // quadratic with closed-form vertex x = 1.5
    let ctx = mathr::eval::Context::standard();
    let r = mathr::repl::dispatch_str("minimize x^2 - 3*x + 2 x 0 5", ctx).unwrap().unwrap();
    assert!(r.contains("x* = 1.5"), "got: {}", r);
    assert!(r.contains("= -0.25"), "got: {}", r);

    // sin on [3, 5]: minimum at 3pi/2 ~ 4.712
    let ctx = mathr::eval::Context::standard();
    let r = mathr::repl::dispatch_str("minimize sin(x) x 3 5", ctx).unwrap().unwrap();
    assert!(r.contains("4.7123"), "got: {}", r);

    // library-level Nelder-Mead on Rosenbrock
    let rosen = |p: &[f64]| {
        let (a, b) = (p[0] - 1.0, p[0] * p[0] - p[1]);
        a * a + 100.0 * b * b
    };
    let res = mathr::optim::nelder_mead(rosen, &[-1.2, 1.0], &mathr::optim::OptOptions::default())
        .unwrap();
    assert!(res.converged);
    assert!(close(res.x[0], 1.0, 1e-4) && close(res.x[1], 1.0, 1e-4));
}