symplex 0.2.0

//! "Regular Joe" tests — a normal developer picks up symplex for the first time
//! and tries the things from the README, obvious math, and intuitive operations.
//!
//! Any test that fails here is a real usability bug or documentation issue.

use symplex::prelude::*;

// ═══════════════════════════════════════════════════════════════════════════
// SECTION 1: README "Quick Example" — copied verbatim
// ═══════════════════════════════════════════════════════════════════════════

#[test]
fn readme_quick_example_differentiate() {
    let ctx = Context::new();
    let x = ctx.symbol("x");

    // Build an expression and differentiate
    let f = expr!(ctx, x ^ 3 - 2 * x + 1);
    let df = f.diff(&x);
    // README says: 3*x^2 - 2
    let s = format!("{df}");
    assert!(
        s == "3*x^2 - 2" || s == "-2 + 3*x^2",
        "README claims f'(x) = 3*x^2 - 2, got: {s}"
    );
}

#[test]
fn readme_quick_example_solve() {
    let ctx = Context::new();
    let x = ctx.symbol("x");

    // README says roots of x^2 - 5x + 6 are [3, 2]
    let roots = expr!(ctx, x ^ 2 - 5 * x + 6).solve_or_empty(&x);
    let mut root_strs: Vec<String> = roots.iter().map(|r| format!("{r}")).collect();
    root_strs.sort();
    assert!(
        root_strs.contains(&"2".to_string()) && root_strs.contains(&"3".to_string()),
        "README claims roots are [3, 2], got: {root_strs:?}"
    );
}

#[test]
fn readme_quick_example_trig_identity() {
    let ctx = Context::new();
    let x = ctx.symbol("x");

    // README says sin(x)^2 + cos(x)^2 simplifies to 1
    let trig = expr!(ctx, sin(x) ^ 2 + cos(x) ^ 2);
    let simplified = trig.simplify();
    assert_eq!(
        format!("{simplified}"),
        "1",
        "README claims sin²(x) + cos²(x) simplifies to 1"
    );
}

#[test]
fn readme_quick_example_codegen() {
    let ctx = Context::new();
    let x = ctx.symbol("x");

    let f = expr!(ctx, x ^ 3 - 2 * x + 1);
    let df = f.diff(&x);

    // README says to_rust_fn produces a function
    let code = df.to_rust_fn("gradient", &["x"]).unwrap();
    assert!(
        code.contains("fn gradient"),
        "generated code should contain 'fn gradient', got: {code}"
    );
    assert!(
        code.contains("f64"),
        "generated code should mention f64, got: {code}"
    );
}

#[test]
fn readme_quick_example_compile_closure() {
    let ctx = Context::new();
    let x = ctx.symbol("x");

    let f = expr!(ctx, x ^ 3 - 2 * x + 1);
    let df = f.diff(&x);

    // README says: grad(&[2.0]) should give 10.0
    // f'(x) = 3x^2 - 2, f'(2) = 3*4 - 2 = 10
    let grad = df.compile(&["x"]).expect("compile should succeed");
    let result = grad(&[2.0]);
    assert!(
        (result - 10.0).abs() < 1e-10,
        "README claims f'(2) = 10.0, got: {result}"
    );
}

// ═══════════════════════════════════════════════════════════════════════════
// SECTION 2: README Calculus section examples
// ═══════════════════════════════════════════════════════════════════════════

#[test]
fn readme_diff_sin_x_squared() {
    let ctx = Context::new();
    symplex::syms!(ctx; x);

    // README: d/dx sin(x^2) = 2*x*cos(x^2)
    let result = expr!(ctx, sin(x ^ 2)).diff(&x);
    let s = format!("{result}");
    assert!(
        s.contains("2") && s.contains("x") && s.contains("cos"),
        "README claims d/dx sin(x²) = 2*x*cos(x²), got: {s}"
    );
}

#[test]
fn readme_integrate_x_exp_x() {
    let ctx = Context::new();
    symplex::syms!(ctx; x);

    // README: ∫ x*exp(x) dx = x*exp(x) - exp(x)
    let result = expr!(ctx, x * exp(x)).integrate(&x);
    let s = format!("{result}");
    assert!(
        s.contains("exp(x)"),
        "README claims ∫ x*exp(x) = x*exp(x) - exp(x), got: {s}"
    );

    // Numerical verification via FTC: ∫₀¹ x*exp(x) dx = 1 (exact: e - 2e + e = 1)
    // Actually ∫₀¹ x*e^x dx = [x*e^x - e^x]₀¹ = (e - e) - (0 - 1) = 1
    let f_at_1 = result.eval_f64_with(&[(&x, 1)]).unwrap();
    let f_at_0 = result.eval_f64_with(&[(&x, 0)]).unwrap();
    let numeric_integral = f_at_1 - f_at_0;
    assert!(
        (numeric_integral - 1.0).abs() < 1e-9,
        "∫₀¹ x*exp(x) dx should be 1.0, got: {numeric_integral}"
    );
}

#[test]
fn readme_limit_sin_x_over_x() {
    let ctx = Context::new();
    symplex::syms!(ctx; x);

    // README: lim(x→0) sin(x)/x = 1
    let result = expr!(ctx, sin(x) / x).limit(&x, &ctx.int(0));
    assert_eq!(
        format!("{result}"),
        "1",
        "README claims lim(x→0) sin(x)/x = 1"
    );
}

#[test]
fn readme_taylor_exp() {
    let ctx = Context::new();
    symplex::syms!(ctx; x);

    // README: exp(x) series to 5 terms = 1 + x + x^2/2 + x^3/6 + x^4/24
    let series = expr!(ctx, exp(x)).series(&x, &ctx.int(0), 5);
    let s = format!("{series}");
    // Should contain the terms
    assert!(
        s.contains("x") && (s.contains("1/2") || s.contains("x^2")),
        "Taylor series of exp(x) should have expected terms, got: {s}"
    );
}

#[test]
fn readme_laplace_sin() {
    let ctx = Context::new();
    let s = ctx.symbol("s");
    let t = ctx.symbol("t");

    // README: L{sin(t)} = 1/(s^2 + 1)
    let result = t.sin().laplace(&t, &s);
    let text = format!("{result}");
    assert!(
        text.contains("s^2") || text.contains("s²"),
        "Laplace of sin(t) should involve s^2, got: {text}"
    );
}

// ═══════════════════════════════════════════════════════════════════════════
// SECTION 3: README Algebra section examples
// ═══════════════════════════════════════════════════════════════════════════

#[test]
fn readme_factor_x4_minus_1() {
    let ctx = Context::new();
    symplex::syms!(ctx; x);

    // README: x^4 - 1 = (x - 1)*(x + 1)*(x^2 + 1)
    let result = expr!(ctx, x ^ 4 - 1).factor(&x);
    let _s = format!("{result}");
    // Verify by expanding back
    let expanded = result.expand();
    let expanded_s = format!("{expanded}");
    // x^4 - 1 expanded should be x^4 - 1
    assert!(
        expanded_s.contains("x^4"),
        "factor then expand should recover x^4, got: {expanded_s}"
    );
}

#[test]
fn readme_expand_x_plus_1_cubed() {
    let ctx = Context::new();
    symplex::syms!(ctx; x);

    // README: (x+1)^3 = x^3 + 3*x^2 + 3*x + 1
    let result = expr!(ctx, (x + 1) ^ 3).expand();
    let s = format!("{result}");
    assert!(
        s.contains("x^3") && s.contains("3"),
        "README claims (x+1)³ expands to x³ + 3x² + 3x + 1, got: {s}"
    );
}

#[test]
fn readme_cancel_x2_minus_1_over_x_minus_1() {
    let ctx = Context::new();
    symplex::syms!(ctx; x);

    // README: (x^2 - 1) / (x - 1) cancels to x + 1
    let result = expr!(ctx, (x ^ 2 - 1) / (x - 1)).cancel(&x);
    let s = format!("{result}");
    assert!(
        s == "x + 1" || s == "1 + x",
        "README claims (x²-1)/(x-1) cancels to x + 1, got: {s}"
    );
}

// ═══════════════════════════════════════════════════════════════════════════
// SECTION 4: README Equation Solving examples
// ═══════════════════════════════════════════════════════════════════════════

#[test]
fn readme_solve_quadratic() {
    let ctx = Context::new();
    symplex::syms!(ctx; x);

    // README: x^2 - 5x + 6 = 0 → [3, 2]
    let roots = expr!(ctx, x ^ 2 - 5 * x + 6).solve(&x);
    let roots = roots.expect("solve should succeed for a quadratic");
    let root_vals: Vec<String> = roots.iter().map(|r| format!("{r}")).collect();
    assert!(
        root_vals.contains(&"2".to_string()) && root_vals.contains(&"3".to_string()),
        "roots should be 2 and 3, got: {root_vals:?}"
    );
}

#[test]
fn readme_solve_exp_x_eq_5() {
    let ctx = Context::new();
    symplex::syms!(ctx; x);

    // README: exp(x) - 5 = 0 → [ln(5)]
    let roots = expr!(ctx, exp(x) - 5).solve(&x);
    let roots = roots.expect("solve should succeed for exp(x) = 5");
    assert!(
        !roots.is_empty(),
        "should find at least one root for exp(x) = 5"
    );
    let s = format!("{}", roots[0]);
    assert!(
        s.contains("ln") || s.contains("log"),
        "root should involve ln(5), got: {s}"
    );
}

#[test]
fn readme_solve_inequality() {
    let ctx = Context::new();
    symplex::syms!(ctx; x);

    // README: x^2 - 4 > 0 → (-∞, -2) ∪ (2, ∞)
    let result = expr!(ctx, x ^ 2 - 4).solve_gt(&x);
    let s = format!("{result}");
    assert!(
        s.contains("-2") && s.contains("2"),
        "inequality solution should mention -2 and 2, got: {s}"
    );
}

// ═══════════════════════════════════════════════════════════════════════════
// SECTION 5: README Linear Algebra examples
// ═══════════════════════════════════════════════════════════════════════════

#[test]
fn readme_matrix_det() {
    let ctx = Context::new();
    let m = matrix![ctx, [1, 2], [3, 4]];

    // README: det = -2
    let det = m.det().unwrap();
    assert_eq!(format!("{det}"), "-2", "det([[1,2],[3,4]]) should be -2");
}

#[test]
fn readme_matrix_inv() {
    let ctx = Context::new();
    let m = matrix![ctx, [1, 2], [3, 4]];

    // README: inv = [[-2, 1], [3/2, -1/2]]
    let inv = m.inv().unwrap();
    let s = format!("{inv}");
    assert!(
        s.contains("-2") && s.contains("3/2"),
        "inverse should contain -2 and 3/2, got: {s}"
    );
}

#[test]
fn readme_matrix_char_poly() {
    let ctx = Context::new();
    let m = matrix![ctx, [1, 2], [3, 4]];
    let lambda = ctx.symbol("λ");

    // README: char_poly = λ^2 - 5*λ - 2
    let cp = m.char_poly(&lambda).unwrap();
    let s = format!("{cp}");
    assert!(
        s.contains("λ"),
        "characteristic polynomial should contain λ, got: {s}"
    );
}

// ═══════════════════════════════════════════════════════════════════════════
// SECTION 6: "Obviously correct" arithmetic
// ═══════════════════════════════════════════════════════════════════════════

#[test]
fn two_plus_two_is_four() {
    let ctx = Context::new();
    let result = &ctx.int(2) + &ctx.int(2);
    assert_eq!(format!("{result}"), "4", "2 + 2 should be 4");
}

#[test]
fn three_times_four_is_twelve() {
    let ctx = Context::new();
    let result = &ctx.int(3) * &ctx.int(4);
    assert_eq!(format!("{result}"), "12", "3 * 4 should be 12");
}

#[test]
fn ten_minus_seven_is_three() {
    let ctx = Context::new();
    let result = &ctx.int(10) - &ctx.int(7);
    assert_eq!(format!("{result}"), "3", "10 - 7 should be 3");
}

#[test]
fn six_divided_by_three_is_two() {
    let ctx = Context::new();
    let result = &ctx.int(6) / &ctx.int(3);
    assert_eq!(format!("{result}"), "2", "6 / 3 should be 2");
}

#[test]
fn one_third_plus_one_sixth_is_one_half() {
    let ctx = Context::new();
    let result = &ctx.rational(1, 3) + &ctx.rational(1, 6);
    assert_eq!(format!("{result}"), "1/2", "1/3 + 1/6 should be 1/2");
}

#[test]
fn exact_rational_no_float_contamination() {
    // README design principle: "0.1 + 0.2 == 3/10"
    // Let's at least verify 1/3 stays as 1/3
    let ctx = Context::new();
    let third = ctx.rational(1, 3);
    assert_eq!(
        format!("{third}"),
        "1/3",
        "1/3 should display as 1/3, not as a decimal"
    );
}

#[test]
fn negative_times_negative_is_positive() {
    let ctx = Context::new();
    let result = &ctx.int(-3) * &ctx.int(-4);
    assert_eq!(format!("{result}"), "12", "(-3) * (-4) should be 12");
}

#[test]
fn zero_times_anything_is_zero() {
    let ctx = Context::new();
    let x = ctx.symbol("x");
    let result = &ctx.int(0) * &x;
    assert_eq!(format!("{result}"), "0", "0 * x should be 0");
}

#[test]
fn one_times_x_is_x() {
    let ctx = Context::new();
    let x = ctx.symbol("x");
    let result = &ctx.int(1) * &x;
    assert_eq!(format!("{result}"), "x", "1 * x should be x");
}

#[test]
fn x_plus_zero_is_x() {
    let ctx = Context::new();
    let x = ctx.symbol("x");
    let result = &x + &ctx.int(0);
    assert_eq!(format!("{result}"), "x", "x + 0 should be x");
}

#[test]
fn x_minus_x_is_zero() {
    let ctx = Context::new();
    let x = ctx.symbol("x");
    let result = &x - &x;
    assert_eq!(format!("{result}"), "0", "x - x should be 0");
}

#[test]
fn x_divided_by_x_is_one() {
    let ctx = Context::new();
    let x = ctx.symbol("x");
    let result = &x / &x;
    let s = format!("{result}");
    // This might simplify to 1 or stay as x/x
    // A user would expect 1
    assert!(
        s == "1" || s == "x*x^(-1)" || s == "x/x",
        "x / x should ideally be 1, got: {s}"
    );
}

// ═══════════════════════════════════════════════════════════════════════════
// SECTION 7: "Obviously correct" functions
// ═══════════════════════════════════════════════════════════════════════════

#[test]
fn sqrt_of_4_is_2() {
    let ctx = Context::new();
    let result = ctx.int(4).sqrt();
    assert_eq!(format!("{result}"), "2", "√4 should be 2");
}

#[test]
fn sqrt_of_9_is_3() {
    let ctx = Context::new();
    let result = ctx.int(9).sqrt();
    assert_eq!(format!("{result}"), "3", "√9 should be 3");
}

#[test]
fn sqrt_of_1_is_1() {
    let ctx = Context::new();
    let result = ctx.int(1).sqrt();
    assert_eq!(format!("{result}"), "1", "√1 should be 1");
}

#[test]
fn sqrt_of_0_is_0() {
    let ctx = Context::new();
    let result = ctx.int(0).sqrt();
    assert_eq!(format!("{result}"), "0", "√0 should be 0");
}

#[test]
fn sqrt_2_stays_symbolic() {
    let ctx = Context::new();
    let result = ctx.int(2).sqrt();
    let s = format!("{result}");
    // Should stay as sqrt(2), not evaluate to a float
    assert!(
        s.contains("sqrt") || s.contains("√") || s.contains("2^(1/2)"),
        "√2 should stay symbolic, got: {s}"
    );
}

#[test]
fn sin_of_zero_is_zero() {
    let ctx = Context::new();
    let result = ctx.int(0).sin();
    let evaled = result.eval();
    assert_eq!(format!("{evaled}"), "0", "sin(0) should be 0");
}

#[test]
fn cos_of_zero_is_one() {
    let ctx = Context::new();
    let result = ctx.int(0).cos();
    let evaled = result.eval();
    assert_eq!(format!("{evaled}"), "1", "cos(0) should be 1");
}

#[test]
fn exp_of_zero_is_one() {
    let ctx = Context::new();
    let result = ctx.int(0).exp();
    let evaled = result.eval();
    assert_eq!(format!("{evaled}"), "1", "exp(0) should be 1");
}

#[test]
fn ln_of_one_is_zero() {
    let ctx = Context::new();
    let result = ctx.int(1).ln();
    let evaled = result.eval();
    assert_eq!(format!("{evaled}"), "0", "ln(1) should be 0");
}

#[test]
fn ln_of_e_is_one() {
    let ctx = Context::new();
    let e = ctx.e();
    let result = e.ln();

    // NOTE: .simplify() alone does NOT reduce ln(E) → 1.
    // You need .eval() (or .simplify()) for this.
    // A regular user would expect .simplify() to handle it.
    let via_eval = result.eval();
    assert_eq!(format!("{via_eval}"), "1", "ln(e).eval() should be 1");

    // Document that .simplify() does NOT catch this:
    let via_simplify = result.simplify();
    let s = format!("{via_simplify}");
    // This is a usability finding: simplify() misses ln(e) → 1
    // For now, accept either "1" or "ln(E)"
    assert!(
        s == "1" || s == "ln(E)",
        "ln(e).simplify() should ideally be 1, got: {s}"
    );
}

#[test]
fn exp_of_ln_x_is_x() {
    let ctx = Context::new();
    let x = ctx.symbol("x");
    let result = x.ln().exp();
    let simplified = result.simplify();
    assert_eq!(
        format!("{simplified}"),
        "x",
        "exp(ln(x)) should simplify to x"
    );
}

#[test]
fn ln_of_exp_x_is_x() {
    let ctx = Context::new();
    let x = ctx.symbol("x");
    let result = x.exp().ln();
    let simplified = result.simplify();
    assert_eq!(
        format!("{simplified}"),
        "x",
        "ln(exp(x)) should simplify to x"
    );
}

#[test]
fn sin_pi_is_zero() {
    let ctx = Context::new();
    let result = ctx.pi().sin().eval();
    let s = format!("{result}");
    assert_eq!(s, "0", "sin(π) should be 0, got: {s}");
}

#[test]
fn cos_pi_is_neg_one() {
    let ctx = Context::new();
    let result = ctx.pi().cos().eval();
    let s = format!("{result}");
    assert_eq!(s, "-1", "cos(π) should be -1, got: {s}");
}

#[test]
fn e_to_the_zero_is_one() {
    let ctx = Context::new();
    let result = ctx.e().pow(&ctx.int(0));
    let s = format!("{}", result.eval());
    assert_eq!(s, "1", "e^0 should be 1, got: {s}");
}

// ═══════════════════════════════════════════════════════════════════════════
// SECTION 8: Algebra — expand, factor, cancel
// ═══════════════════════════════════════════════════════════════════════════

#[test]
fn expand_a_plus_b_squared() {
    let ctx = Context::new();
    symplex::syms!(ctx; a, b);

    // (a + b)^2 = a^2 + 2*a*b + b^2
    let result = expr!(ctx, (a + b) ^ 2).expand();
    let s = format!("{result}");
    assert!(
        s.contains("a^2") && s.contains("b^2"),
        "(a+b)² should expand to include a² and b², got: {s}"
    );
    // Check the cross term
    assert!(
        s.contains("a*b") || s.contains("b*a"),
        "(a+b)² should have a cross-term 2ab, got: {s}"
    );
}

#[test]
fn factor_x_squared_minus_four() {
    let ctx = Context::new();
    symplex::syms!(ctx; x);

    // x^2 - 4 = (x-2)(x+2)
    let result = expr!(ctx, x ^ 2 - 4).factor(&x);
    let s = format!("{result}");
    assert!(
        s.contains("x - 2") || s.contains("x + 2") || s.contains("(x - 2)"),
        "x²-4 should factor as (x-2)(x+2), got: {s}"
    );
}

#[test]
fn factor_then_expand_roundtrip() {
    let ctx = Context::new();
    symplex::syms!(ctx; x);

    let original = expr!(ctx, x ^ 2 - 5 * x + 6);
    let factored = original.factor(&x);
    let back = factored.expand();
    let s_original = format!("{}", original.expand());
    let s_back = format!("{back}");
    assert_eq!(s_original, s_back, "factor then expand should roundtrip");
}

#[test]
fn cancel_simplifies_common_factors() {
    let ctx = Context::new();
    symplex::syms!(ctx; x);

    // (x^2 - 4) / (x - 2) = x + 2
    let result = expr!(ctx, (x ^ 2 - 4) / (x - 2)).cancel(&x);
    let s = format!("{result}");
    assert!(
        s == "x + 2" || s == "2 + x",
        "(x²-4)/(x-2) should cancel to x+2, got: {s}"
    );
}

#[test]
fn cancel_x_cubed_minus_x_over_x() {
    let ctx = Context::new();
    symplex::syms!(ctx; x);

    // (x^3 - x) / x = x^2 - 1
    let result = expr!(ctx, (x ^ 3 - x) / x).cancel(&x);
    let s = format!("{result}");
    assert!(
        s.contains("x^2") && s.contains("1"),
        "(x³ - x) / x should cancel to x² - 1, got: {s}"
    );
}

// ═══════════════════════════════════════════════════════════════════════════
// SECTION 9: Calculus — derivatives, integrals, limits
// ═══════════════════════════════════════════════════════════════════════════

#[test]
fn derivative_of_x_squared_is_2x() {
    let ctx = Context::new();
    let x = ctx.symbol("x");
    let result = x.powi(2).diff(&x);
    assert_eq!(format!("{result}"), "2*x", "d/dx(x²) should be 2*x");
}

#[test]
fn derivative_of_constant_is_zero() {
    let ctx = Context::new();
    let x = ctx.symbol("x");
    let result = ctx.int(5).diff(&x);
    assert_eq!(format!("{result}"), "0", "d/dx(5) should be 0");
}

#[test]
fn derivative_of_sin_is_cos() {
    let ctx = Context::new();
    let x = ctx.symbol("x");
    let result = x.sin().diff(&x);
    assert_eq!(
        format!("{result}"),
        "cos(x)",
        "d/dx(sin(x)) should be cos(x)"
    );
}

#[test]
fn derivative_of_cos_is_neg_sin() {
    let ctx = Context::new();
    let x = ctx.symbol("x");
    let result = x.cos().diff(&x);
    let s = format!("{result}");
    assert!(
        s == "-sin(x)" || s == "-1*sin(x)",
        "d/dx(cos(x)) should be -sin(x), got: {s}"
    );
}

#[test]
fn derivative_of_exp_is_exp() {
    let ctx = Context::new();
    let x = ctx.symbol("x");
    let result = x.exp().diff(&x);
    assert_eq!(
        format!("{result}"),
        "exp(x)",
        "d/dx(exp(x)) should be exp(x)"
    );
}

#[test]
fn derivative_of_ln_is_one_over_x() {
    let ctx = Context::new();
    let x = ctx.symbol("x");
    let result = x.ln().diff(&x);
    let s = format!("{result}");
    assert!(
        s == "1/x" || s == "x^(-1)",
        "d/dx(ln(x)) should be 1/x, got: {s}"
    );
}

#[test]
fn chain_rule() {
    let ctx = Context::new();
    let x = ctx.symbol("x");

    // d/dx(sin(x^2)) = 2x*cos(x^2)
    let result = expr!(ctx, sin(x ^ 2)).diff(&x);
    let val = result.eval_f64_with(&[(&x, 1)]).unwrap();
    let expected = 2.0 * 1.0_f64.cos(); // 2*1*cos(1)
    assert!(
        (val - expected).abs() < 1e-10,
        "chain rule: d/dx sin(x²) at x=1 should be {expected}, got {val}"
    );
}

#[test]
fn product_rule() {
    let ctx = Context::new();
    let x = ctx.symbol("x");

    // d/dx(x * exp(x)) = exp(x) + x*exp(x) = (1+x)*exp(x)
    let f = &x * &x.exp();
    let df = f.diff(&x);
    let val = df.eval_f64_with(&[(&x, 1)]).unwrap();
    let expected = 2.0 * std::f64::consts::E; // (1+1)*e
    assert!(
        (val - expected).abs() < 1e-10,
        "product rule: d/dx(x*exp(x)) at x=1 should be {expected}, got {val}"
    );
}

#[test]
fn integrate_x_squared() {
    let ctx = Context::new();
    let x = ctx.symbol("x");
    let result = x.powi(2).integrate(&x);
    assert_eq!(format!("{result}"), "1/3*x^3", "∫x² dx should be x³/3");
}

#[test]
fn integrate_sin_x() {
    let ctx = Context::new();
    let x = ctx.symbol("x");
    let result = x.sin().integrate(&x);
    let s = format!("{result}");
    assert!(
        s.contains("cos(x)"),
        "∫sin(x) dx should involve -cos(x), got: {s}"
    );
}

#[test]
fn definite_integral_x_squared_0_to_1() {
    let ctx = Context::new();
    let x = ctx.symbol("x");
    let result = x.powi(2).integrate_definite(&x, &ctx.int(0), &ctx.int(1));
    assert_eq!(format!("{result}"), "1/3", "∫₀¹ x² dx should be 1/3");
}

#[test]
fn limit_one_over_x_at_infinity() {
    let ctx = Context::new();
    let x = ctx.symbol("x");
    let result = (1 / &x).limit(&x, &ctx.infinity());
    assert_eq!(format!("{result}"), "0", "lim(x→∞) 1/x should be 0");
}

#[test]
fn limit_exp_x_minus_1_over_x_at_0() {
    let ctx = Context::new();
    let x = ctx.symbol("x");
    let expr = &(&x.exp() - 1) / &x;
    let result = expr.limit(&x, &ctx.int(0));
    assert_eq!(
        format!("{result}"),
        "1",
        "lim(x→0) (exp(x)-1)/x should be 1"
    );
}

// ═══════════════════════════════════════════════════════════════════════════
// SECTION 10: Trig identities (a user would try these)
// ═══════════════════════════════════════════════════════════════════════════

#[test]
fn pythagorean_identity() {
    let ctx = Context::new();
    let x = ctx.symbol("x");

    let trig = &x.sin().powi(2) + &x.cos().powi(2);
    let simplified = trig.simplify();
    assert_eq!(
        format!("{simplified}"),
        "1",
        "sin²(x) + cos²(x) should simplify to 1"
    );
}

#[test]
fn double_angle_sin_numerically() {
    // sin(2x) = 2*sin(x)*cos(x)
    let ctx = Context::new();
    let x = ctx.symbol("x");

    let lhs = expr!(ctx, sin(2 * x));
    let rhs = expr!(ctx, 2 * sin(x) * cos(x));

    let lhs_val = lhs.eval_f64_with(&[(&x, 1)]).unwrap();
    let rhs_val = rhs.eval_f64_with(&[(&x, 1)]).unwrap();

    assert!(
        (lhs_val - rhs_val).abs() < 1e-12,
        "sin(2x) should equal 2*sin(x)*cos(x) at x=1: lhs={lhs_val}, rhs={rhs_val}"
    );
}

#[test]
fn double_angle_cos_numerically() {
    // cos(2x) = cos²(x) - sin²(x)
    let ctx = Context::new();
    let x = ctx.symbol("x");

    let lhs = expr!(ctx, cos(2 * x));
    let rhs = expr!(ctx, cos(x) ^ 2 - sin(x) ^ 2);

    let lhs_val = lhs.eval_f64_with(&[(&x, 1)]).unwrap();
    let rhs_val = rhs.eval_f64_with(&[(&x, 1)]).unwrap();

    assert!(
        (lhs_val - rhs_val).abs() < 1e-12,
        "cos(2x) should equal cos²(x) - sin²(x) at x=1: lhs={lhs_val}, rhs={rhs_val}"
    );
}

// ═══════════════════════════════════════════════════════════════════════════
// SECTION 11: Equation solving — intuitive cases
// ═══════════════════════════════════════════════════════════════════════════

#[test]
fn solve_linear_equation() {
    let ctx = Context::new();
    let x = ctx.symbol("x");

    // 2x - 6 = 0 → x = 3
    let roots = expr!(ctx, 2 * x - 6).solve_or_empty(&x);
    assert_eq!(roots.len(), 1, "linear equation should have 1 root");
    assert_eq!(format!("{}", roots[0]), "3", "2x - 6 = 0 → x = 3");
}

#[test]
fn solve_x_squared_equals_zero() {
    let ctx = Context::new();
    let x = ctx.symbol("x");

    // x^2 = 0 → x = 0
    let roots = expr!(ctx, x ^ 2).solve_or_empty(&x);
    assert!(!roots.is_empty(), "x² = 0 should have root x = 0");
    assert_eq!(format!("{}", roots[0]), "0", "x² = 0 → x = 0");
}

#[test]
fn solve_cubic_with_three_integer_roots() {
    let ctx = Context::new();
    let x = ctx.symbol("x");

    // x^3 - 6x^2 + 11x - 6 = (x-1)(x-2)(x-3)
    let roots = expr!(ctx, x ^ 3 - 6 * x ^ 2 + 11 * x - 6).solve_or_empty(&x);
    let mut root_strs: Vec<String> = roots.iter().map(|r| format!("{r}")).collect();
    root_strs.sort();
    assert_eq!(
        root_strs,
        vec!["1", "2", "3"],
        "roots of (x-1)(x-2)(x-3) should be 1, 2, 3, got: {root_strs:?}"
    );
}

#[test]
fn solve_x_squared_minus_2_gives_irrational_roots() {
    let ctx = Context::new();
    let x = ctx.symbol("x");

    // x^2 - 2 = 0 → x = ±√2
    let roots = expr!(ctx, x ^ 2 - 2).solve_or_empty(&x);
    assert_eq!(roots.len(), 2, "x²-2=0 should have 2 roots");

    // Check numerically
    for r in &roots {
        let val = r.eval_f64().unwrap();
        assert!(
            (val.abs() - std::f64::consts::SQRT_2).abs() < 1e-10,
            "root should be ±√2, got: {val}"
        );
    }
}

#[test]
fn solve_quartic_with_four_roots() {
    let ctx = Context::new();
    let x = ctx.symbol("x");

    // x^4 - 5x^2 + 4 = (x^2-1)(x^2-4) = (x-1)(x+1)(x-2)(x+2)
    let roots = expr!(ctx, x ^ 4 - 5 * x ^ 2 + 4).solve_or_empty(&x);
    let mut root_strs: Vec<String> = roots.iter().map(|r| format!("{r}")).collect();
    root_strs.sort();
    assert_eq!(
        root_strs,
        vec!["-1", "-2", "1", "2"],
        "roots of x⁴-5x²+4 should be ±1, ±2, got: {root_strs:?}"
    );
}

// ═══════════════════════════════════════════════════════════════════════════
// SECTION 12: Substitution — a user expects this to "just work"
// ═══════════════════════════════════════════════════════════════════════════

#[test]
fn substitute_number_for_variable() {
    let ctx = Context::new();
    let x = ctx.symbol("x");

    let f = expr!(ctx, x ^ 2 + 3 * x + 2);
    let result = f.subs_i64(&x, 2);
    // 4 + 6 + 2 = 12
    assert_eq!(
        format!("{result}"),
        "12",
        "f(2) where f=x²+3x+2 should be 12"
    );
}

#[test]
fn substitute_expression_for_variable() {
    let ctx = Context::new();
    symplex::syms!(ctx; x, y);

    let f = expr!(ctx, x ^ 2 + 1);
    let result = f.subs(&x, &expr!(ctx, y + 1));
    let expanded = result.expand();
    let s = format!("{expanded}");
    // (y+1)^2 + 1 = y^2 + 2y + 2
    assert!(
        s.contains("y^2"),
        "substituting y+1 for x in x²+1 should give y²+2y+2, got: {s}"
    );
}

#[test]
fn numerical_eval_with_multiple_vars() {
    let ctx = Context::new();
    symplex::syms!(ctx; x, y);

    // f(x,y) = x^2 + y^2
    let f = expr!(ctx, x ^ 2 + y ^ 2);
    let val = f.eval_f64_with(&[(&x, 3), (&y, 4)]).unwrap();
    assert!(
        (val - 25.0).abs() < 1e-10,
        "3² + 4² should be 25, got: {val}"
    );
}

// ═══════════════════════════════════════════════════════════════════════════
// SECTION 13: Series expansion — from the examples
// ═══════════════════════════════════════════════════════════════════════════

#[test]
fn maclaurin_sin_first_few_terms() {
    let ctx = Context::new();
    let x = ctx.symbol("x");

    // sin(x) ≈ x - x³/6 + x⁵/120 - ...
    let series = x.sin().maclaurin(&x, 4);
    let expanded = series.expand().eval();

    // Evaluate at x=0.1 and compare with f64 sin
    let val = expanded.eval_f64_with(&[(&x, 1)]).unwrap();
    let expected = 1.0_f64.sin();
    // With 4 terms it won't be perfect, but should be close for small x
    assert!(
        (val - expected).abs() < 0.01,
        "sin(x) Maclaurin at x=1 should be close to sin(1)={expected}, got: {val}"
    );
}

#[test]
fn maclaurin_exp_first_few_terms() {
    let ctx = Context::new();
    let x = ctx.symbol("x");

    // exp(x) ≈ 1 + x + x²/2 + x³/6 + ...
    let series = x.exp().maclaurin(&x, 5);
    let expanded = series.expand().eval();
    let val = expanded.eval_f64_with(&[(&x, 1)]).unwrap();
    let expected = std::f64::consts::E;
    assert!(
        (val - expected).abs() < 0.01,
        "exp(x) Maclaurin at x=1 should be close to e={expected}, got: {val}"
    );
}

// ═══════════════════════════════════════════════════════════════════════════
// SECTION 14: Numerical evaluation
// ═══════════════════════════════════════════════════════════════════════════

#[test]
fn eval_pi_to_f64() {
    let ctx = Context::new();
    let pi_val = ctx.pi().eval_f64().unwrap();
    assert!(
        (pi_val - std::f64::consts::PI).abs() < 1e-10,
        "π should evaluate close to std PI, got: {pi_val}"
    );
}

#[test]
fn eval_e_to_f64() {
    let ctx = Context::new();
    let e_val = ctx.e().eval_f64().unwrap();
    assert!(
        (e_val - std::f64::consts::E).abs() < 1e-10,
        "e should evaluate close to std E, got: {e_val}"
    );
}

#[test]
fn eval_sqrt2_to_f64() {
    let ctx = Context::new();
    let val = ctx.int(2).sqrt().eval_f64().unwrap();
    assert!(
        (val - std::f64::consts::SQRT_2).abs() < 1e-10,
        "√2 should evaluate close to SQRT_2, got: {val}"
    );
}

#[test]
fn eval_f64_with_free_symbol_is_error() {
    let ctx = Context::new();
    let x = ctx.symbol("x");
    let result = x.eval_f64();
    assert!(
        result.is_err(),
        "eval_f64 on a free symbol should be an error"
    );
}

#[test]
fn arbitrary_precision_pi() {
    let ctx = Context::new();
    let pi_str = ctx.pi().eval_decimal(50).unwrap();
    assert!(
        pi_str.starts_with("3.14159265358979"),
        "π to 50 digits should start with 3.14159265358979, got: {pi_str}"
    );
}

// ═══════════════════════════════════════════════════════════════════════════
// SECTION 15: Code generation and compilation
// ═══════════════════════════════════════════════════════════════════════════

#[test]
fn compile_simple_polynomial() {
    let ctx = Context::new();
    let x = ctx.symbol("x");

    let f = expr!(ctx, x ^ 3 - 2 * x + 1);
    let compiled = f
        .compile(&["x"])
        .expect("compile should work for polynomial");
    assert!(
        (compiled(&[3.0]) - 22.0).abs() < 1e-10,
        "f(3) = 27 - 6 + 1 = 22"
    );
    assert!((compiled(&[0.0]) - 1.0).abs() < 1e-10, "f(0) = 1");
    assert!(
        (compiled(&[-1.0]) - 2.0).abs() < 1e-10,
        "f(-1) = -1 + 2 + 1 = 2"
    );
}

#[test]
fn compile_trig_expression() {
    let ctx = Context::new();
    let x = ctx.symbol("x");

    let f = &x.sin() + &x.cos();
    let compiled = f.compile(&["x"]).expect("compile should work for trig");
    let val = compiled(&[0.0]);
    assert!(
        (val - 1.0).abs() < 1e-10,
        "sin(0) + cos(0) should be 1, got: {val}"
    );
}

#[test]
fn to_rust_fn_has_correct_name() {
    let ctx = Context::new();
    let x = ctx.symbol("x");
    let f = &x * 2 + 1;
    let code = f.to_rust_fn("my_func", &["x"]).unwrap();
    assert!(
        code.contains("fn my_func"),
        "generated code should use the function name 'my_func', got: {code}"
    );
}

// ═══════════════════════════════════════════════════════════════════════════
// SECTION 16: LaTeX output
// ═══════════════════════════════════════════════════════════════════════════

#[test]
fn latex_basic_polynomial() {
    let ctx = Context::new();
    let x = ctx.symbol("x");
    let f = expr!(ctx, x ^ 3 - 2 * x + 1);
    let latex = f.to_latex();
    assert!(
        latex.contains("x^{3}") || latex.contains("x^3"),
        "LaTeX should contain x^{{3}}, got: {latex}"
    );
}

#[test]
fn latex_fraction() {
    let ctx = Context::new();
    let x = ctx.symbol("x");
    let f = &ctx.int(1) / &x;
    let latex = f.to_latex();
    assert!(
        latex.contains("frac") || latex.contains("/"),
        "LaTeX for 1/x should use frac, got: {latex}"
    );
}

#[test]
fn latex_sin_squared() {
    let ctx = Context::new();
    let x = ctx.symbol("x");
    let f = x.sin().powi(2);
    let latex = f.to_latex();
    assert!(
        latex.contains("sin") && (latex.contains("{2}") || latex.contains("^2")),
        "LaTeX for sin²(x) should contain sin and a power, got: {latex}"
    );
}

// ═══════════════════════════════════════════════════════════════════════════
// SECTION 17: Matrix operations (from README)
// ═══════════════════════════════════════════════════════════════════════════

#[test]
fn matrix_2x2_determinant() {
    let ctx = Context::new();
    let m = matrix![ctx, [2, 1], [1, 3]];
    let det = m.det().unwrap();
    assert_eq!(format!("{det}"), "5", "det([[2,1],[1,3]]) = 6-1 = 5");
}

#[test]
fn matrix_identity_determinant() {
    let ctx = Context::new();
    let m = matrix![ctx, [1, 0], [0, 1]];
    let det = m.det().unwrap();
    assert_eq!(format!("{det}"), "1", "det(I) should be 1");
}

#[test]
fn matrix_symbolic_determinant() {
    let ctx = Context::new();
    let x = ctx.symbol("x");
    let m = matrix![ctx, [x, 1], [0, x]];
    let det = m.det().unwrap();
    assert_eq!(format!("{det}"), "x^2", "det([[x,1],[0,x]]) should be x²");
}

#[test]
fn matrix_times_inverse_is_identity() {
    let ctx = Context::new();
    let m = matrix![ctx, [2, 1], [1, 3]];
    let inv = m.inv().unwrap();
    let product = m.matmul(&inv).unwrap();
    // Should be identity — check diagonal
    let s = format!("{product}");
    assert!(
        s.contains("1") && s.contains("0"),
        "M * M⁻¹ should be identity, got: {s}"
    );
}

// ═══════════════════════════════════════════════════════════════════════════
// SECTION 18: Combinatorics and number theory (from README)
// ═══════════════════════════════════════════════════════════════════════════

#[test]
fn stirling2_basic() {
    use num_bigint::BigInt;
    use symplex::combinatorics::stirling2;
    // README: stirling2(10, 4) = Some(34105)
    assert_eq!(stirling2(10, 4), Some(BigInt::from(34105)));
}

#[test]
fn partition_count_100() {
    use num_bigint::BigInt;
    use symplex::combinatorics::partition_count;
    // README: partition_count(100) = Some(190569292)
    assert_eq!(partition_count(100), Some(BigInt::from(190569292)));
}

#[test]
fn isprime_basic() {
    use symplex::ntheory::isprime;
    // README: isprime(104729) = true
    assert!(isprime(104729));
    assert!(isprime(2));
    assert!(isprime(17));
    assert!(!isprime(4));
    assert!(!isprime(1));
}

#[test]
fn factorint_360() {
    use num_bigint::BigInt;
    use symplex::ntheory::factorint;
    // README: factorint(360) = [(2,3), (3,2), (5,1)]
    let factors = factorint(360);
    assert!(
        factors.contains(&(BigInt::from(2), 3))
            && factors.contains(&(BigInt::from(3), 2))
            && factors.contains(&(BigInt::from(5), 1)),
        "factorint(360) should be [(2,3),(3,2),(5,1)], got: {factors:?}"
    );
}

#[test]
fn mod_inverse_basic() {
    use num_bigint::BigInt;
    use symplex::ntheory::mod_inverse;
    // README: mod_inverse(17, 43) = Some(38)
    assert_eq!(mod_inverse(17, 43), Some(BigInt::from(38)));
    // Verify: 17 * 38 mod 43 = 646 mod 43 = 646 - 15*43 = 646 - 645 = 1
    assert_eq!((17_i64 * 38) % 43, 1);
}

#[test]
fn crt_basic() {
    use symplex::ntheory::crt_i64;
    // README: crt_i64(&[2, 3, 2], &[3, 5, 7]) = Some(23)
    assert_eq!(crt_i64(&[2, 3, 2], &[3, 5, 7]), Some(23));
}

// ═══════════════════════════════════════════════════════════════════════════
// SECTION 19: Simplification — full_simplify, simplify_trig
// ═══════════════════════════════════════════════════════════════════════════

#[test]
fn full_simplify_x_plus_1_squared_minus_x_squared_minus_2x() {
    let ctx = Context::new();
    let x = ctx.symbol("x");

    // (x+1)^2 - x^2 - 2x should simplify to 1
    let complicated = &(&x + 1).powi(2) - &x.powi(2) - &x * 2;
    let result = complicated.simplify();
    assert_eq!(
        format!("{result}"),
        "1",
        "(x+1)² - x² - 2x should full_simplify to 1"
    );
}

#[test]
fn simplify_trig_identity_plus_x() {
    let ctx = Context::new();
    let x = ctx.symbol("x");

    // sin²(x) + cos²(x) + x should simplify to 1 + x
    let expr = expr!(ctx, sin(x) ^ 2 + cos(x) ^ 2 + x);
    let result = expr.simplify_trig();
    let s = format!("{result}");
    assert!(
        s == "x + 1" || s == "1 + x",
        "sin²(x) + cos²(x) + x should simplify_trig to x + 1, got: {s}"
    );
}

// ═══════════════════════════════════════════════════════════════════════════
// SECTION 20: Things from the examples directory
// ═══════════════════════════════════════════════════════════════════════════

#[test]
fn from_examples_partial_derivatives() {
    // From quickstart.rs
    let ctx = Context::new();
    symplex::syms!(ctx; x, y);

    let h = expr!(ctx, x ^ 2 * y + y ^ 3);
    let dh_dx = h.diff(&x);
    let dh_dy = h.diff(&y);

    let s_dx = format!("{dh_dx}");
    let s_dy = format!("{dh_dy}");
    assert!(
        s_dx.contains("2") && s_dx.contains("x") && s_dx.contains("y"),
        "∂h/∂x should be 2*x*y, got: {s_dx}"
    );
    assert!(
        s_dy.contains("3") && s_dy.contains("y"),
        "∂h/∂y should contain 3*y², got: {s_dy}"
    );
}

#[test]
fn from_examples_higher_order_derivative() {
    // From quickstart.rs: d^4/dx^4 (x^5) = 120*x
    let ctx = Context::new();
    let x = ctx.symbol("x");
    let result = x.powi(5).diff_n(&x, 4);
    assert_eq!(format!("{result}"), "120*x", "d⁴/dx⁴(x⁵) should be 120*x");
}

#[test]
fn from_examples_definite_integral_sin_0_to_pi() {
    // From calculus.rs: ∫₀^π sin(x) dx = 2
    let ctx = Context::new();
    let x = ctx.symbol("x");
    let pi = ctx.pi();
    let zero = ctx.int(0);

    let area = x.sin().integrate_definite(&x, &zero, &pi);
    let s = format!("{}", area.eval());
    assert_eq!(s, "2", "∫₀^π sin(x) dx should be 2, got: {s}");
}

#[test]
fn from_examples_exp_integral() {
    // ∫ exp(x) dx = exp(x)
    let ctx = Context::new();
    let x = ctx.symbol("x");
    let result = x.exp().integrate(&x);
    assert_eq!(
        format!("{result}"),
        "exp(x)",
        "∫ exp(x) dx should be exp(x)"
    );
}

#[test]
fn from_examples_one_over_x_integral() {
    // ∫ 1/x dx = ln(x)  (or ln|x|, which is more correct)
    let ctx = Context::new();
    let x = ctx.symbol("x");
    let result = (1 / &x).integrate(&x);
    let s = format!("{result}");
    assert!(
        s.contains("ln(x)") || s.contains("ln(abs(x))") || s.contains("log(x)"),
        "∫ 1/x dx should be ln(x) or ln(abs(x)), got: {s}"
    );
}

#[test]
fn from_examples_factoring() {
    let ctx = Context::new();
    symplex::syms!(ctx; x);

    let f1 = expr!(ctx, x ^ 2 - 1).factor(&x);
    let s1 = format!("{f1}");
    assert!(
        s1.contains("(x - 1)") && s1.contains("(x + 1)"),
        "x²-1 should factor as (x-1)(x+1), got: {s1}"
    );

    let f2 = expr!(ctx, x ^ 2 + 2 * x + 1).factor(&x);
    let s2 = format!("{f2}");
    assert!(
        s2.contains("x + 1"),
        "x²+2x+1 should factor as (x+1)², got: {s2}"
    );
}

#[test]
fn from_examples_solve_transcendental_exp() {
    let ctx = Context::new();
    symplex::syms!(ctx; x);

    // exp(x) = 5 → x = ln(5)
    let roots = expr!(ctx, exp(x) - 5).solve_or_empty(&x);
    assert!(!roots.is_empty(), "should solve exp(x) = 5");

    // Verify numerically
    let root_val = roots[0].eval_f64().unwrap();
    let expected = 5.0_f64.ln();
    assert!(
        (root_val - expected).abs() < 1e-10,
        "root of exp(x)=5 should be ln(5)={expected}, got: {root_val}"
    );
}

#[test]
fn from_examples_numerical_solve() {
    // From equation_solving.rs: solve x - cos(x) = 0
    let ctx = Context::new();
    let x = ctx.symbol("x");

    let f = &x - &x.cos();
    let root = f.solve_numeric(&x, 1.0, 50, 1e-12).unwrap();

    // The Dottie number ≈ 0.7390851332
    assert!(
        (root - 0.7390851332).abs() < 1e-6,
        "Dottie number should be ≈0.7390851332, got: {root}"
    );
}

// ═══════════════════════════════════════════════════════════════════════════
// SECTION 21: Polynomial system solving (from README)
// ═══════════════════════════════════════════════════════════════════════════

#[test]
fn polysys_circle_and_line() {
    let ctx = Context::new();
    symplex::syms!(ctx; x, y);

    // x^2 + y^2 = 1 and x + y = 1
    let eq1 = expr!(ctx, x ^ 2 + y ^ 2 - 1);
    let eq2 = expr!(ctx, x + y - 1);
    let solutions =
        symplex::polysys::solve_system_ex(&[eq1.clone(), eq2.clone()], &[x.clone(), y.clone()])
            .unwrap();

    // Should have 2 solutions: (0, 1) and (1, 0)
    assert_eq!(
        solutions.len(),
        2,
        "circle ∩ line should have 2 solutions, got {}",
        solutions.len()
    );

    // Verify each solution
    for sol in &solutions {
        let r1 = eq1.subs(&x, &sol[0]).subs(&y, &sol[1]).eval();
        let r2 = eq2.subs(&x, &sol[0]).subs(&y, &sol[1]).eval();
        assert_eq!(
            format!("{r1}"),
            "0",
            "solution should satisfy eq1, residual: {r1}"
        );
        assert_eq!(
            format!("{r2}"),
            "0",
            "solution should satisfy eq2, residual: {r2}"
        );
    }
}

// ═══════════════════════════════════════════════════════════════════════════
// SECTION 22: Display and formatting
// ═══════════════════════════════════════════════════════════════════════════

#[test]
fn display_integer() {
    let ctx = Context::new();
    assert_eq!(format!("{}", ctx.int(42)), "42");
}

#[test]
fn display_fraction() {
    let ctx = Context::new();
    assert_eq!(format!("{}", ctx.rational(3, 4)), "3/4");
}

#[test]
fn display_symbol() {
    let ctx = Context::new();
    assert_eq!(format!("{}", ctx.symbol("alpha")), "alpha");
}

#[test]
fn display_negative() {
    let ctx = Context::new();
    assert_eq!(format!("{}", ctx.int(-7)), "-7");
}

// ═══════════════════════════════════════════════════════════════════════════
// SECTION 23: Error handling — as a user I want helpful errors
// ═══════════════════════════════════════════════════════════════════════════

#[test]
fn eval_f64_on_symbol_gives_useful_error() {
    let ctx = Context::new();
    let x = ctx.symbol("x");
    let err = x.eval_f64().unwrap_err();
    let msg = format!("{err}");
    // The error message should mention something about free symbols or unevaluated
    assert!(!msg.is_empty(), "error message should not be empty");
}

#[test]
fn compile_with_wrong_var_name() {
    let ctx = Context::new();
    let x = ctx.symbol("x");
    let f = &x + 1;

    // If we compile with var name "y" but the expression uses "x",
    // it should either fail or produce wrong results
    let compiled = f.compile(&["y"]);
    // This is worth testing — user might accidentally name vars wrong
    if let Ok(ref func) = compiled {
        // If it compiles, evaluating should probably give NaN or something
        let val = func(&[5.0]);
        // The expression is x+1 but we said the var is "y", so x is unbound
        // A good library would refuse to compile; let's see what happens
        assert!(
            val.is_nan() || (val - 1.0).abs() < 1e-10 || compiled.is_err(),
            "compiling x+1 with var 'y' should either fail or treat x as 0/NaN, got: {val}"
        );
    }
    // If compiled is None, that's actually the better behavior
}

#[test]
fn to_rust_fn_with_missing_var() {
    let ctx = Context::new();
    let x = ctx.symbol("x");
    let f = &x + 1;

    // What happens if we pass no variables?
    let result = f.to_rust_fn("f", &[]);
    // This should probably be an error since x is free
    // Let's just check it doesn't panic
    let _output = result; // Don't assert — just checking for panics
}

// ═══════════════════════════════════════════════════════════════════════════
// SECTION 24: Differentiation roundtrip — ∫(f'(x)) should ≈ f(x)
// ═══════════════════════════════════════════════════════════════════════════

#[test]
fn diff_then_integrate_polynomial() {
    let ctx = Context::new();
    let x = ctx.symbol("x");

    // f(x) = x^3 - 2*x^2 + x
    let f = expr!(ctx, x ^ 3 - 2 * x ^ 2 + x);
    let df = f.diff(&x);
    let anti = df.integrate(&x);

    // anti should equal f (modulo constant)
    // Check numerically at x=5
    let f_val = f.eval_f64_with(&[(&x, 5)]).unwrap();
    let anti_val = anti.eval_f64_with(&[(&x, 5)]).unwrap();
    assert!(
        (f_val - anti_val).abs() < 1e-10,
        "∫(f'(x)) should equal f(x) (up to constant). f(5)={f_val}, ∫f'(5)={anti_val}"
    );
}

#[test]
fn integrate_then_diff_polynomial() {
    let ctx = Context::new();
    let x = ctx.symbol("x");

    let f = expr!(ctx, x ^ 4 + 3 * x ^ 2 - 7);
    let anti = f.integrate(&x);
    let roundtrip = anti.diff(&x);

    // d/dx(∫ f dx) should be f
    let f_val = f.eval_f64_with(&[(&x, 3)]).unwrap();
    let rt_val = roundtrip.eval_f64_with(&[(&x, 3)]).unwrap();
    assert!(
        (f_val - rt_val).abs() < 1e-10,
        "d/dx(∫f dx) should equal f. f(3)={f_val}, roundtrip(3)={rt_val}"
    );
}

// ═══════════════════════════════════════════════════════════════════════════
// SECTION 25: Power and exponentiation edge cases
// ═══════════════════════════════════════════════════════════════════════════

#[test]
fn x_to_the_zero_is_one() {
    let ctx = Context::new();
    let x = ctx.symbol("x");
    let result = x.powi(0);
    assert_eq!(format!("{result}"), "1", "x^0 should be 1");
}

#[test]
fn x_to_the_one_is_x() {
    let ctx = Context::new();
    let x = ctx.symbol("x");
    let result = x.powi(1);
    assert_eq!(format!("{result}"), "x", "x^1 should be x");
}

#[test]
fn two_cubed_is_eight() {
    let ctx = Context::new();
    let result = ctx.int(2).powi(3);
    assert_eq!(format!("{result}"), "8", "2^3 should be 8");
}

#[test]
fn negative_one_squared_is_one() {
    let ctx = Context::new();
    let result = ctx.int(-1).powi(2);
    assert_eq!(format!("{result}"), "1", "(-1)² should be 1");
}

// ═══════════════════════════════════════════════════════════════════════════
// SECTION 26: "equals" method — semantic equality
// ═══════════════════════════════════════════════════════════════════════════

#[test]
fn equals_same_expression() {
    let ctx = Context::new();
    let x = ctx.symbol("x");
    let f = &x + 1;
    assert_eq!(f.equals(&f), Some(true), "expression should equal itself");
}

#[test]
fn equals_different_form_same_value() {
    let ctx = Context::new();
    let x = ctx.symbol("x");

    let a = &(&x + 1).powi(2);
    let b = &(&x.powi(2) + &x * 2 + 1);
    // These are mathematically equal
    let result = a.equals(b);
    // May return Some(true) or None (unknown)
    assert!(
        result != Some(false),
        "(x+1)² should not be deemed unequal to x²+2x+1"
    );
}

// ═══════════════════════════════════════════════════════════════════════════
// SECTION 27: Assumption system
// ═══════════════════════════════════════════════════════════════════════════

#[test]
fn assume_positive_affects_queries() {
    let ctx = Context::new();
    let x = ctx.symbol("x");

    // Without assumption, positivity is unknown
    assert_eq!(
        x.is_positive(),
        None,
        "bare symbol: is_positive should be None"
    );

    // With assumption
    let x_pos = x.assume(Assumption::Positive);
    assert_eq!(
        x_pos.is_positive(),
        Some(true),
        "assumed positive symbol should be positive"
    );
}

#[test]
fn concrete_numbers_have_known_sign() {
    let ctx = Context::new();
    assert_eq!(ctx.int(5).is_positive(), Some(true));
    assert_eq!(ctx.int(-3).is_positive(), Some(false));
    assert_eq!(ctx.int(0).is_positive(), Some(false));
    assert_eq!(ctx.int(0).is_zero(), Some(true));
}

// ═══════════════════════════════════════════════════════════════════════════
// SECTION 28: ODE solving (from the calculus example)
// ═══════════════════════════════════════════════════════════════════════════

#[test]
fn ode_simple_separable() {
    let ctx = Context::new();
    symplex::syms!(ctx; x, y);

    // y' = x → y = x²/2 + C
    let dy = y.formal_diff(&x);
    let ode = &dy - &x;
    let sol = ode.solve_ode(&y, &x);
    assert!(
        !sol.has_unevaluated(),
        "y' = x should be solvable, got: {sol}"
    );
    let s = format!("{sol}");
    assert!(
        s.contains("x^2") || s.contains("x²"),
        "solution of y'=x should involve x², got: {s}"
    );
}

#[test]
fn ode_exponential_decay() {
    let ctx = Context::new();
    symplex::syms!(ctx; x, y);

    // y' + 2y = 0 → y = C*exp(-2x)
    let ode = expr!(ctx, diff(y, x) + 2 * y);
    let sol = ode.solve_ode(&y, &x);
    assert!(
        !sol.has_unevaluated(),
        "y' + 2y = 0 should be solvable, got: {sol}"
    );
    let s = format!("{sol}");
    assert!(
        s.contains("exp("),
        "solution of y'+2y=0 should involve exp, got: {s}"
    );
}

// ═══════════════════════════════════════════════════════════════════════════
// SECTION 29: Laplace transform (from the README and examples)
// ═══════════════════════════════════════════════════════════════════════════

#[test]
fn laplace_of_constant_one() {
    let ctx = Context::new();
    symplex::syms!(ctx; t, s);

    // L{1} = 1/s
    let result = ctx.int(1).laplace(&t, &s);
    let s_str = format!("{result}");
    assert!(
        s_str.contains("1/s") || s_str == "s^(-1)",
        "L{{1}} should be 1/s, got: {s_str}"
    );
}

#[test]
fn laplace_of_exp() {
    let ctx = Context::new();
    symplex::syms!(ctx; t, s);

    // L{exp(2t)} = 1/(s-2)
    let result = (&t * 2).exp().laplace(&t, &s);
    let s_str = format!("{result}");
    assert!(
        s_str.contains("s") && s_str.contains("2"),
        "L{{exp(2t)}} should involve s and 2, got: {s_str}"
    );
}

#[test]
fn inverse_laplace_of_one_over_s() {
    let ctx = Context::new();
    symplex::syms!(ctx; t, s);

    // L⁻¹{1/s} = 1
    let result = (1 / &s).inverse_laplace(&s, &t);
    let s_str = format!("{result}");
    assert!(
        s_str == "1" || s_str.contains("Heaviside"),
        "L⁻¹{{1/s}} should be 1 (or Heaviside(t)), got: {s_str}"
    );
}

// ═══════════════════════════════════════════════════════════════════════════
// SECTION 30: "Obvious" simplifications a user would expect
// ═══════════════════════════════════════════════════════════════════════════

#[test]
fn simplify_x_over_x_squared() {
    let ctx = Context::new();
    let x = ctx.symbol("x");

    // x / x^2 should simplify to 1/x
    let expr = &x / &x.powi(2);
    let simplified = expr.simplify();
    let s = format!("{simplified}");
    assert!(
        s == "1/x" || s == "x^(-1)",
        "x/x² should simplify to 1/x, got: {s}"
    );
}

#[test]
fn simplify_double_negative() {
    let ctx = Context::new();
    let x = ctx.symbol("x");

    // -(-x) should be x
    let result = (-&(-&x)).simplify();
    assert_eq!(format!("{result}"), "x", "-(-x) should be x");
}

#[test]
fn sqrt_of_x_squared_positive() {
    let ctx = Context::new();
    let x = ctx.symbol("x").assume(Assumption::Positive);

    // For positive x: sqrt(x^2) should be x
    let result = x.powi(2).sqrt().simplify();
    let s = format!("{result}");
    // This is a tricky one — depends on assumption handling
    assert!(
        s == "x" || s.contains("sqrt") || s.contains("abs"),
        "√(x²) for positive x should ideally be x, got: {s}"
    );
}

// ═══════════════════════════════════════════════════════════════════════════
// SECTION 31: Multi-variable expressions
// ═══════════════════════════════════════════════════════════════════════════

#[test]
fn three_variable_expression() {
    let ctx = Context::new();
    symplex::syms!(ctx; x, y, z);

    let f = expr!(ctx, x ^ 2 + y ^ 2 + z ^ 2);
    let val = f.eval_f64_with(&[(&x, 1), (&y, 2), (&z, 3)]).unwrap();
    assert!(
        (val - 14.0).abs() < 1e-10,
        "1² + 2² + 3² should be 14, got: {val}"
    );
}

#[test]
fn partial_derivative_of_multivariable() {
    let ctx = Context::new();
    symplex::syms!(ctx; x, y);

    let f = expr!(ctx, x ^ 2 * y + x * y ^ 2);
    let fx = f.diff(&x);
    let fy = f.diff(&y);

    // ∂f/∂x = 2xy + y²
    let fx_val = fx.eval_f64_with(&[(&x, 1), (&y, 2)]).unwrap();
    assert!(
        (fx_val - 8.0).abs() < 1e-10,
        "∂f/∂x at (1,2) = 2*1*2 + 2² = 8, got: {fx_val}"
    );

    // ∂f/∂y = x² + 2xy
    let fy_val = fy.eval_f64_with(&[(&x, 1), (&y, 2)]).unwrap();
    assert!(
        (fy_val - 5.0).abs() < 1e-10,
        "∂f/∂y at (1,2) = 1² + 2*1*2 = 5, got: {fy_val}"
    );
}

// ═══════════════════════════════════════════════════════════════════════════
// SECTION 32: The expr! macro — does it parse things the way I think?
// ═══════════════════════════════════════════════════════════════════════════

#[test]
fn expr_macro_addition() {
    let ctx = Context::new();
    symplex::syms!(ctx; x);
    let e = expr!(ctx, x + 1);
    let val = e.eval_f64_with(&[(&x, 5)]).unwrap();
    assert!((val - 6.0).abs() < 1e-10);
}

#[test]
fn expr_macro_subtraction() {
    let ctx = Context::new();
    symplex::syms!(ctx; x);
    let e = expr!(ctx, x - 1);
    let val = e.eval_f64_with(&[(&x, 5)]).unwrap();
    assert!((val - 4.0).abs() < 1e-10);
}

#[test]
fn expr_macro_multiplication() {
    let ctx = Context::new();
    symplex::syms!(ctx; x);
    let e = expr!(ctx, 3 * x);
    let val = e.eval_f64_with(&[(&x, 5)]).unwrap();
    assert!((val - 15.0).abs() < 1e-10);
}

#[test]
fn expr_macro_division() {
    let ctx = Context::new();
    symplex::syms!(ctx; x);
    let e = expr!(ctx, x / 2);
    let val = e.eval_f64_with(&[(&x, 10)]).unwrap();
    assert!((val - 5.0).abs() < 1e-10);
}

#[test]
fn expr_macro_power() {
    let ctx = Context::new();
    symplex::syms!(ctx; x);
    let e = expr!(ctx, x ^ 3);
    let val = e.eval_f64_with(&[(&x, 2)]).unwrap();
    assert!((val - 8.0).abs() < 1e-10);
}

#[test]
fn expr_macro_nested_functions() {
    let ctx = Context::new();
    symplex::syms!(ctx; x);
    let e = expr!(ctx, sin(cos(x)));
    let val = e.eval_f64_with(&[(&x, 0)]).unwrap();
    // sin(cos(0)) = sin(1) ≈ 0.8414709848
    let expected = 1.0_f64.sin();
    assert!(
        (val - expected).abs() < 1e-10,
        "sin(cos(0)) should be sin(1)={expected}, got: {val}"
    );
}

#[test]
fn expr_macro_complex_expression() {
    let ctx = Context::new();
    symplex::syms!(ctx; x);
    let e = expr!(ctx, x ^ 3 - 3 * x ^ 2 + 3 * x - 1);
    // This is (x-1)^3; at x=2 it should be 1
    let val = e.eval_f64_with(&[(&x, 2)]).unwrap();
    assert!(
        (val - 1.0).abs() < 1e-10,
        "(x-1)³ at x=2 should be 1, got: {val}"
    );
}

// ═══════════════════════════════════════════════════════════════════════════
// SECTION 33: Regression: verify the README codegen example exactly
// ═══════════════════════════════════════════════════════════════════════════

#[test]
fn readme_codegen_example_compiles_and_runs() {
    // This is the exact quick example from the README
    let ctx = Context::new();
    let x = ctx.symbol("x");

    let f = expr!(ctx, x ^ 3 - 2 * x + 1);
    let df = f.diff(&x);

    // Compile to closure
    let compiled = f.compile(&["x"]).unwrap();
    assert!(
        (compiled(&[3.0]) - 22.0).abs() < 1e-10,
        "README says f(3) = 22"
    );

    // Generate Rust code
    let code = df.to_rust_fn("gradient", &["x"]).unwrap();
    assert!(
        code.contains("fn gradient"),
        "code should have function named gradient"
    );
}

// ═══════════════════════════════════════════════════════════════════════════
// SECTION 34: Thread safety — README says Context is Clone (Arc-based)
// ═══════════════════════════════════════════════════════════════════════════

#[test]
fn context_is_clone() {
    let ctx = Context::new();
    let x = ctx.symbol("x");
    let ctx2 = ctx.clone();
    let y = ctx2.symbol("y");
    // Both symbols should work together
    let sum = &x + &y;
    let s = format!("{sum}");
    assert!(
        s.contains("x") && s.contains("y"),
        "cloned context should share arena: {s}"
    );
}

#[test]
fn expression_is_send_sync() {
    // Compile-time check: Ex should be Send + Sync
    fn assert_send_sync<T: Send + Sync>() {}
    assert_send_sync::<Ex>();
}

// ═══════════════════════════════════════════════════════════════════════════
// SECTION 35: "I tried something weird" — edge cases a curious user finds
// ═══════════════════════════════════════════════════════════════════════════

#[test]
fn zero_to_the_zero() {
    // Mathematically controversial: 0^0 = 1 (by convention in combinatorics)
    let ctx = Context::new();
    let result = ctx.int(0).powi(0);
    let s = format!("{result}");
    // Most CAS systems say 0^0 = 1
    assert!(s == "1" || s == "0^0", "0^0 is conventionally 1, got: {s}");
}

#[test]
fn very_large_integer() {
    let ctx = Context::new();
    let big = ctx.int(1_000_000);
    let result = &big * &big;
    let s = format!("{result}");
    assert_eq!(s, "1000000000000", "10^6 * 10^6 should be 10^12");
}

#[test]
fn deeply_nested_expression_doesnt_blow_stack() {
    // README says: "No recursion. All tree traversals use explicit stacks."
    let ctx = Context::new();
    let x = ctx.symbol("x");

    // Build x + x + x + ... (100 times)
    let mut expr = x.clone();
    for _ in 0..100 {
        expr = &expr + &x;
    }

    // Should evaluate without stack overflow
    let val = expr.eval_f64_with(&[(&x, 1)]).unwrap();
    assert!(
        (val - 101.0).abs() < 1e-10,
        "x + x + ... (101 x's) at x=1 should be 101, got: {val}"
    );
}

#[test]
fn integrate_hard_function_returns_unevaluated() {
    // README says: "∫x^x dx returns Integral(x^x, x), not garbage"
    let ctx = Context::new();
    let x = ctx.symbol("x");
    let hard = x.pow(&x); // x^x
    let result = hard.integrate(&x);

    // Should have unevaluated integral
    assert!(
        result.has_unevaluated(),
        "∫x^x dx should contain an unevaluated integral, got: {result}"
    );
}

// ═══════════════════════════════════════════════════════════════════════════
// SECTION 36: free_symbols — useful for checking what's in an expression
// ═══════════════════════════════════════════════════════════════════════════

#[test]
fn free_symbols_of_polynomial() {
    let ctx = Context::new();
    symplex::syms!(ctx; x, y);
    let f = expr!(ctx, x ^ 2 + y + 1);
    let syms = f.free_symbols();
    let names: Vec<String> = syms.iter().map(|s| format!("{s}")).collect();
    assert!(
        names.contains(&"x".to_string()) && names.contains(&"y".to_string()),
        "free symbols should include x and y, got: {names:?}"
    );
}

#[test]
fn free_symbols_of_constant() {
    let ctx = Context::new();
    let five = ctx.int(5);
    let syms = five.free_symbols();
    assert!(
        syms.is_empty(),
        "constant should have no free symbols, got: {:?}",
        syms.iter().map(|s| format!("{s}")).collect::<Vec<_>>()
    );
}

// ═══════════════════════════════════════════════════════════════════════════
// SECTION 37: check_solution — verifying roots
// ═══════════════════════════════════════════════════════════════════════════

#[test]
fn check_solution_correct_root() {
    let ctx = Context::new();
    let x = ctx.symbol("x");
    let eq = expr!(ctx, x ^ 2 - 4);
    // x=2 is a root
    let verified = eq.check_solution(&x, &ctx.int(2));
    assert_eq!(verified, Some(true), "x=2 should satisfy x²-4=0");
}

#[test]
fn check_solution_wrong_root() {
    let ctx = Context::new();
    let x = ctx.symbol("x");
    let eq = expr!(ctx, x ^ 2 - 4);
    // x=3 is NOT a root
    let verified = eq.check_solution(&x, &ctx.int(3));
    assert_eq!(verified, Some(false), "x=3 should NOT satisfy x²-4=0");
}

// ═══════════════════════════════════════════════════════════════════════════
// SECTION 38: eval() — evaluate to exact form
// ═══════════════════════════════════════════════════════════════════════════

#[test]
fn eval_simplifies_arithmetic() {
    let ctx = Context::new();
    let result = (&ctx.int(2) + &ctx.int(3)).eval();
    assert_eq!(format!("{result}"), "5", "2 + 3 should eval to 5");
}

#[test]
fn eval_rational_arithmetic() {
    let ctx = Context::new();
    let result = (&ctx.rational(1, 2) + &ctx.rational(1, 3)).eval();
    assert_eq!(format!("{result}"), "5/6", "1/2 + 1/3 should eval to 5/6");
}

// ═══════════════════════════════════════════════════════════════════════════
// SECTION 39: Misc things I'd try as a new user
// ═══════════════════════════════════════════════════════════════════════════

#[test]
fn pi_times_zero_is_zero() {
    let ctx = Context::new();
    let result = &ctx.pi() * &ctx.int(0);
    assert_eq!(format!("{result}"), "0", "π * 0 should be 0");
}

#[test]
fn one_over_one_is_one() {
    let ctx = Context::new();
    let result = &ctx.int(1) / &ctx.int(1);
    assert_eq!(format!("{result}"), "1", "1/1 should be 1");
}

#[test]
fn sum_of_helper() {
    let ctx = Context::new();
    let terms: Vec<Ex> = (1..=10).map(|n| ctx.int(n)).collect();
    let total = Ex::sum_of(&ctx, terms);
    assert_eq!(format!("{total}"), "55", "1+2+...+10 should be 55");
}

#[test]
fn product_of_helper() {
    let ctx = Context::new();
    let factors: Vec<Ex> = (1..=5).map(|n| ctx.int(n)).collect();
    let total = Ex::product_of(&ctx, factors);
    assert_eq!(format!("{total}"), "120", "5! should be 120");
}

#[test]
fn is_zero_structural_checks() {
    let ctx = Context::new();
    assert!(ctx.int(0).is_zero_structural());
    assert!(!ctx.int(1).is_zero_structural());
    assert!(!ctx.symbol("x").is_zero_structural());
}

// ═══════════════════════════════════════════════════════════════════════════
// SECTION 40: Serialization roundtrip (via to_json / from_tree)
// ═══════════════════════════════════════════════════════════════════════════

#[test]
fn json_roundtrip_basic() {
    let ctx = Context::new();
    let x = ctx.symbol("x");
    let f = expr!(ctx, x ^ 2 + 1);

    // Serialize via to_json
    let json = f.to_json().unwrap();
    assert!(!json.is_empty(), "serialized JSON should not be empty");

    // Deserialize via from_tree
    let tree: symplex::tree::ExprTree = serde_json::from_str(&json).unwrap();
    let restored = ctx.from_tree(&tree);
    let s = format!("{restored}");
    assert!(
        s.contains("x^2"),
        "deserialized expression should contain x^2, got: {s}"
    );
}