symplex 0.22.1

//! 0.22 — `integrate` regressions from the 0.21 differential audit.
//!
//! Four confirmed wrong antiderivatives: `∫ sqrt(x²) dx` treated as
//! `∫ x dx`; the Weierstrass substitution applying a numeric factor twice;
//! the same substitution producing bogus `atan(3·tan(x/2))` forms for
//! `∫ cos x·R(sin x) dx` (rational integrals with an irreducible quartic
//! denominator); and `∫ 1/((x²−4)²+1) dx` emitting an `re(…)`/`im(…)`
//! closed form whose derivative is not the integrand.
//!
//! The right assertion for an antiderivative is `d/dx F = f`, so every
//! test differentiates the result and compares with the integrand at
//! four rational points (relative tolerance 1e-10).  The reference value
//! at each point is `f` itself, quoted from sympy 1.14
//! (`symplex/.venv/bin/python`) as `N(f.subs(x, p), 17)`; the sympy
//! antiderivative `integrate(f, x)` is cited alongside for the record.
//! Where the simplifier can do it, `(F′ − f).simplify()` is also checked
//! to be zero.  Nothing below was computed by hand.

// The reference values are quoted with sympy's 17 significant digits.
#![allow(clippy::excessive_precision)]

use symplex::prelude::*;

/// The rational sample points `p = n/d` shared by every table below:
/// `Rational(1,3), Rational(-5,7), Rational(13,11), Rational(-17,5)`.
const POINTS: [(i64, i64); 4] = [(1, 3), (-5, 7), (13, 11), (-17, 5)];

/// Relative closeness, absolute when the reference is (essentially) zero.
fn close(actual: f64, expected: f64, rel: f64, label: &str) {
    if expected.abs() < 1e-15 && actual.abs() < 1e-15 {
        return;
    }
    let scale = expected.abs().max(1e-300);
    assert!(
        (actual - expected).abs() <= rel * scale,
        "{label}: got {actual:e}, expected {expected:e} (relative error {:e})",
        (actual - expected).abs() / scale
    );
}

/// `F = ∫ f dx` must be a closed form with `F′(p) = f(p)` at every point
/// of [`POINTS`]; `f_ref[i]` is sympy's `N(f.subs(x, p_i), 17)` and is
/// also checked against symplex's own evaluation of `f` (so a typo in the
/// integrand would be caught, not silently agreed with).  Returns `F`.
fn assert_antiderivative(f: &Ex, x: &Ex, f_ref: &[f64; 4], label: &str) -> Ex {
    let big_f = f.integrate(x);
    assert!(
        !big_f.has_unevaluated(),
        "{label}: expected a closed form, got {big_f}"
    );
    let d_big_f = big_f.diff(x);
    let ctx = x.context();
    for (&(n, d), &expected) in POINTS.iter().zip(f_ref) {
        let p = ctx.rational(n, d);
        let f_at = f
            .subs(x, &p)
            .eval_f64()
            .unwrap_or_else(|e| panic!("{label}: f({n}/{d}) did not evaluate: {e}"));
        close(
            f_at,
            expected,
            1e-12,
            &format!("{label}: f({n}/{d}) vs sympy"),
        );
        let df_at = d_big_f
            .subs(x, &p)
            .eval_f64()
            .unwrap_or_else(|e| panic!("{label}: F'({n}/{d}) did not evaluate: {e} (F = {big_f})"));
        close(
            df_at,
            expected,
            1e-10,
            &format!("{label}: F'({n}/{d}) vs f, F = {big_f}"),
        );
    }
    big_f
}

/// `(F′ − f).simplify()` is literally zero.
fn assert_derivative_simplifies_to_f(f: &Ex, big_f: &Ex, x: &Ex, label: &str) {
    let residual = (big_f.diff(x) - f).simplify();
    assert_eq!(
        residual.to_string(),
        "0",
        "{label}: (F' - f).simplify() should be 0, F = {big_f}"
    );
}

// ═══════════════════════════════════════════════════════════════════════════
// 1. sqrt(x²) is |x|, not x
// ═══════════════════════════════════════════════════════════════════════════

#[test]
fn sqrt_of_x_squared_integrates_as_abs_x() {
    // sympy: integrate(sqrt(x**2), x) = x*sqrt(x**2)/2
    // f = sqrt(x**2): 0.33333333333333333, 0.71428571428571429, 1.1818181818181818, 3.4000000000000000
    // (0.21 returned x**2/2, whose derivative is -5/7 at the second point.)
    let ctx = Context::new();
    let x = ctx.symbol("x");
    let f = x.pow(&ctx.int(2)).sqrt();
    let big_f = assert_antiderivative(
        &f,
        &x,
        &[
            0.33333333333333333,
            0.71428571428571429,
            1.1818181818181818,
            3.4000000000000000,
        ],
        "∫ sqrt(x^2)",
    );
    assert_ne!(big_f.to_string(), "1/2*x^2", "∫ sqrt(x^2) must not be ∫ x");
}

#[test]
fn neg_sqrt_of_x_squared_integrates_as_neg_abs_x() {
    // sympy: integrate(-sqrt(x**2), x) = -x*sqrt(x**2)/2
    // f = -sqrt(x**2): -0.33333333333333333, -0.71428571428571429, -1.1818181818181818, -3.4000000000000000
    let ctx = Context::new();
    let x = ctx.symbol("x");
    let f = -x.pow(&ctx.int(2)).sqrt();
    assert_antiderivative(
        &f,
        &x,
        &[
            -0.33333333333333333,
            -0.71428571428571429,
            -1.1818181818181818,
            -3.4000000000000000,
        ],
        "∫ -sqrt(x^2)",
    );
}

#[test]
fn x_squared_to_three_halves_is_abs_x_cubed() {
    // sympy: integrate((x**2)**Rational(3,2), x) = x*(x**2)**(3/2)/4
    // f = (x**2)**(3/2): 0.037037037037037037, 0.36443148688046647, 1.6506386175807663, 39.304000000000000
    let ctx = Context::new();
    let x = ctx.symbol("x");
    let f = x.pow(&ctx.int(2)).pow(&ctx.rational(3, 2));
    assert_antiderivative(
        &f,
        &x,
        &[
            0.037037037037037037,
            0.36443148688046647,
            1.6506386175807663,
            39.304000000000000,
        ],
        "∫ (x^2)^(3/2)",
    );
}

#[test]
fn sqrt_of_x_fourth_flattens_to_x_squared() {
    // The flattening (x^m)^n → x^(m·n) stays legitimate when m·n is an even
    // integer: sqrt(x**4) = x**2 for every real x.
    // sympy: integrate(sqrt(x**4), x) = x*sqrt(x**4)/3
    // f = sqrt(x**4): 0.11111111111111111, 0.51020408163265306, 1.3966942148760331, 11.560000000000000
    let ctx = Context::new();
    let x = ctx.symbol("x");
    let f = x.pow(&ctx.int(4)).sqrt();
    let big_f = assert_antiderivative(
        &f,
        &x,
        &[
            0.11111111111111111,
            0.51020408163265306,
            1.3966942148760331,
            11.560000000000000,
        ],
        "∫ sqrt(x^4)",
    );
    assert_eq!(big_f.to_string(), "1/3*x^3");
}

#[test]
fn reciprocal_sqrt_flattening_still_reaches_asinh() {
    // The flattening (x²+1)^(1/2) → (x²+1)^(-1/2) that the standard-form
    // handler relies on is unaffected (odd numerator in m = 1/2).
    // sympy: integrate(1/sqrt(x**2+1), x) = asinh(x)
    // f = 1/sqrt(x**2+1): 0.94868329805051380, 0.81373347120673496, 0.64594224146617384, 0.28216632399155017
    let ctx = Context::new();
    let x = ctx.symbol("x");
    let f = ctx.int(1) / (x.pow(&ctx.int(2)) + 1).sqrt();
    let big_f = assert_antiderivative(
        &f,
        &x,
        &[
            0.94868329805051380,
            0.81373347120673496,
            0.64594224146617384,
            0.28216632399155017,
        ],
        "∫ 1/sqrt(x^2+1)",
    );
    assert_eq!(big_f.to_string(), "asinh(x)");
}

// ═══════════════════════════════════════════════════════════════════════════
// 2. Weierstrass: the constant factor is applied exactly once
// ═══════════════════════════════════════════════════════════════════════════

#[test]
fn neg_of_product_sin_plus_one_times_cos() {
    // sympy: integrate(-(sin(x)+1)*cos(x), x) = -sin(x)**2/2 - sin(x)
    // f = -(sin(x)+1)*cos(x): -1.2541418478496062, -0.26060980851463465, -0.73015561900524185, 1.2138548681487652
    // (0.21 returned an F with F' = -f.)
    let ctx = Context::new();
    let x = ctx.symbol("x");
    let f = -((x.sin() + 1) * x.cos());
    let big_f = assert_antiderivative(
        &f,
        &x,
        &[
            -1.2541418478496062,
            -0.26060980851463465,
            -0.73015561900524185,
            1.2138548681487652,
        ],
        "∫ -(sin x + 1) cos x",
    );
    assert_derivative_simplifies_to_f(&f, &big_f, &x, "∫ -(sin x + 1) cos x");
}

#[test]
fn three_times_cos_minus_four_times_sin() {
    // sympy: integrate(3*(cos(x)-4)*sin(x), x) = 3*sin(x)**2/2 + 12*cos(x)
    // f = 3*(cos(x)-4)*sin(x): -2.9987816569492214, 6.3760801515840370, -10.050827316564941, -3.8076632510298883
    // (0.21 returned an F with F' = 3·f.)
    let ctx = Context::new();
    let x = ctx.symbol("x");
    let f = (x.cos() - 4) * x.sin() * 3;
    let big_f = assert_antiderivative(
        &f,
        &x,
        &[
            -2.9987816569492214,
            6.3760801515840370,
            -10.050827316564941,
            -3.8076632510298883,
        ],
        "∫ 3 (cos x - 4) sin x",
    );
    assert_derivative_simplifies_to_f(&f, &big_f, &x, "∫ 3 (cos x - 4) sin x");
}

#[test]
fn minus_two_times_cos_minus_four_times_sin() {
    // sympy: integrate(-2*(cos(x)-4)*sin(x), x) = -sin(x)**2 - 8*cos(x)
    // f = -2*(cos(x)-4)*sin(x): 1.9991877712994809, -4.2507201010560247, 6.7005515443766276, 2.5384421673532589
    // (0.21 returned an F with F' = -2·f.)
    let ctx = Context::new();
    let x = ctx.symbol("x");
    let f = (x.cos() - 4) * x.sin() * -2;
    let big_f = assert_antiderivative(
        &f,
        &x,
        &[
            1.9991877712994809,
            -4.2507201010560247,
            6.7005515443766276,
            2.5384421673532589,
        ],
        "∫ -2 (cos x - 4) sin x",
    );
    assert_derivative_simplifies_to_f(&f, &big_f, &x, "∫ -2 (cos x - 4) sin x");
}

#[test]
fn three_cos_over_sin_squared_plus_one() {
    // sympy: integrate(3*cos(x)/(sin(x)**2+1), x) = 3*atan(sin(x))
    // f = 3*cos(x)/(sin(x)**2+1): 2.5607285380951193, 1.5860619515432745, 0.61294300386480045, -2.7226050514833956
    // (0.21 returned 9*atan(3*tan(x/2)): both the double constant and the
    // wrong quartic integration at once.)
    let ctx = Context::new();
    let x = ctx.symbol("x");
    let f = x.cos() * 3 / (x.sin().pow(&ctx.int(2)) + 1);
    let big_f = assert_antiderivative(
        &f,
        &x,
        &[
            2.5607285380951193,
            1.5860619515432745,
            0.61294300386480045,
            -2.7226050514833956,
        ],
        "∫ 3 cos x/(sin²x + 1)",
    );
    assert_eq!(big_f.to_string(), "3*atan(sin(x))");
}

// ═══════════════════════════════════════════════════════════════════════════
// 3. cos x·R(sin x), sin x·R(cos x): u-substitution wins; a·sin² + b
// ═══════════════════════════════════════════════════════════════════════════

#[test]
fn cos_over_sin_squared_plus_one_is_atan_sin() {
    // sympy: integrate(cos(x)/(sin(x)**2+1), x) = atan(sin(x))
    // f = cos(x)/(sin(x)**2+1): 0.85357617936503976, 0.52868731718109150, 0.20431433462160015, -0.90753501716113186
    // (0.21 returned atan(3*tan(x/2)).)
    let ctx = Context::new();
    let x = ctx.symbol("x");
    let f = x.cos() / (x.sin().pow(&ctx.int(2)) + 1);
    let big_f = assert_antiderivative(
        &f,
        &x,
        &[
            0.85357617936503976,
            0.52868731718109150,
            0.20431433462160015,
            -0.90753501716113186,
        ],
        "∫ cos x/(sin²x + 1)",
    );
    assert_eq!(big_f.to_string(), "atan(sin(x))");
    assert_derivative_simplifies_to_f(&f, &big_f, &x, "∫ cos x/(sin²x + 1)");
}

#[test]
fn minus_two_cos_over_sin_squared_plus_four() {
    // sympy: integrate(-2*cos(x)/(sin(x)**2+4), x) = -atan(sin(x)/2)
    // f = -2*cos(x)/(sin(x)**2+4): -0.46016263779896368, -0.34117844800441081, -0.15619005457761024, 0.47563421846387038
    // (0.21 returned -2*atan(-3*tan(x/2)); at x = 1, f = -0.2295 but F' = 1.0568.)
    let ctx = Context::new();
    let x = ctx.symbol("x");
    let f = x.cos() * -2 / (x.sin().pow(&ctx.int(2)) + 4);
    let big_f = assert_antiderivative(
        &f,
        &x,
        &[
            -0.46016263779896368,
            -0.34117844800441081,
            -0.15619005457761024,
            0.47563421846387038,
        ],
        "∫ -2 cos x/(sin²x + 4)",
    );
    assert_eq!(big_f.to_string(), "-atan(1/2*sin(x))");
}

#[test]
fn three_sin_over_scaled_cos_squared_plus_four() {
    // sympy: integrate(3*sin(x)/(Rational(9,4)*cos(x)**2+4), x) = -atan(3*cos(x)/4)
    // f: 0.16334897170264714, -0.37188892601949924, 0.64203134933413098, 0.12561268880521519
    // (0.21 returned an F with F' = 3·f.)
    let ctx = Context::new();
    let x = ctx.symbol("x");
    let f = x.sin() * 3 / (x.cos().pow(&ctx.int(2)) * ctx.rational(9, 4) + 4);
    let big_f = assert_antiderivative(
        &f,
        &x,
        &[
            0.16334897170264714,
            -0.37188892601949924,
            0.64203134933413098,
            0.12561268880521519,
        ],
        "∫ 3 sin x/(9/4 cos²x + 4)",
    );
    assert_eq!(big_f.to_string(), "-atan(3/4*cos(x))");
}

#[test]
fn one_over_sin_squared_plus_four() {
    // Pure Weierstrass case: t = tan(x/2) gives (1+t²)/(2t⁴+6t²+2), a
    // biquadratic whose roots in t² are real and negative.
    // sympy: integrate(1/(sin(x)**2+4), x) =
    //   521*sqrt(10)*sqrt(3 - sqrt(5))*(atan(sqrt(2)*tan(x/2)/sqrt(3 - sqrt(5))) + pi*floor((x/2 - pi/2)/pi))/(2880*sqrt(5) + 6440) + … (four such terms)
    // f = 1/(sin(x)**2+4): 0.24348338810226434, 0.22577812476394662, 0.20592339995226665, 0.24598423027398144
    // (0.21 returned sqrt(5)/10*atan(3*sqrt(5)/5*tan(x/2)); at x = 1, f = 0.2124 but F' = 0.1267.)
    let ctx = Context::new();
    let x = ctx.symbol("x");
    let f = ctx.int(1) / (x.sin().pow(&ctx.int(2)) + 4);
    assert_antiderivative(
        &f,
        &x,
        &[
            0.24348338810226434,
            0.22577812476394662,
            0.20592339995226665,
            0.24598423027398144,
        ],
        "∫ 1/(sin²x + 4)",
    );
}

#[test]
fn sin_over_cos_squared_plus_one_unchanged() {
    // Was already right in 0.21 (pins the u = cos x route).
    // sympy: integrate(sin(x)/(cos(x)**2+1), x) = -atan(cos(x))
    // f = sin(x)/(cos(x)**2+1): 0.17284967790034225, -0.41701520021118986, 0.80894963134186558, 0.13208314868872689
    let ctx = Context::new();
    let x = ctx.symbol("x");
    let f = x.sin() / (x.cos().pow(&ctx.int(2)) + 1);
    let big_f = assert_antiderivative(
        &f,
        &x,
        &[
            0.17284967790034225,
            -0.41701520021118986,
            0.80894963134186558,
            0.13208314868872689,
        ],
        "∫ sin x/(cos²x + 1)",
    );
    assert_eq!(big_f.to_string(), "-atan(cos(x))");
}

// ═══════════════════════════════════════════════════════════════════════════
// 4. Irreducible quartic denominators with nested-radical roots
// ═══════════════════════════════════════════════════════════════════════════

#[test]
fn one_over_x_squared_minus_four_squared_plus_one() {
    // sympy: integrate(1/((x**2-4)**2+1), x) — an atan/log combination in
    //   sqrt(-1/136 + sqrt(17)/544), sqrt(33 - 8*sqrt(17)), … (too long to quote).
    // f = 1/((x**2-4)**2+1): 0.062021439509954058, 0.075880159281967006, 0.12858096358877979, 0.017195839982391460
    // (0.21 returned a form in re(…)/im(…) with F'(-2) = 0.98879 for f(-2) = 1,
    //  F'(1) = 0.016535 for f(1) = 0.1.)
    let ctx = Context::new();
    let x = ctx.symbol("x");
    let f = ctx.int(1) / ((x.pow(&ctx.int(2)) - 4).pow(&ctx.int(2)) + 1);
    let big_f = assert_antiderivative(
        &f,
        &x,
        &[
            0.062021439509954058,
            0.075880159281967006,
            0.12858096358877979,
            0.017195839982391460,
        ],
        "∫ 1/((x²-4)²+1)",
    );
    let s = big_f.to_string();
    assert!(
        !s.contains("re(") && !s.contains("im("),
        "closed form must not contain opaque re/im constants: {s}"
    );
}

#[test]
fn ln_of_x_squared_minus_four_squared_plus_one() {
    // By parts, then the same quartic: ∫ ln((x²−4)²+1) = x·ln(…) − ∫ 4x²(x²−4)/((x²−4)²+1).
    // sympy: integrate(log((x**2-4)**2+1), x) = x*log((x**2 - 4)**2 + 1) - 4*x
    //   - sqrt(2 + sqrt(17)/2)*log(x**2 - sqrt(2)*x*sqrt(4 + sqrt(17)) - sqrt(8*sqrt(17) + 33) + 4 + 2*sqrt(17)) + … (atan terms)
    // f = log((x**2-4)**2+1): 2.7802751551639376, 2.5786000347877560, 2.0511965062168838, 4.0630877859048051
    // (0.21: F'(-2) = -1.0076 for f(-2) = 0.)
    let ctx = Context::new();
    let x = ctx.symbol("x");
    let f = ((x.pow(&ctx.int(2)) - 4).pow(&ctx.int(2)) + 1).ln();
    let big_f = assert_antiderivative(
        &f,
        &x,
        &[
            2.7802751551639376,
            2.5786000347877560,
            2.0511965062168838,
            4.0630877859048051,
        ],
        "∫ ln((x²-4)²+1)",
    );
    let s = big_f.to_string();
    assert!(
        !s.contains("re(") && !s.contains("im("),
        "closed form must not contain opaque re/im constants: {s}"
    );
}

#[test]
fn biquadratic_with_real_negative_roots_in_x_squared() {
    // x⁴ + 6x² + 1 = (x² + 3 − 2√2)(x² + 3 + 2√2).
    // sympy: integrate(1/(x**4+6*x**2+1), x)
    //   = 2*(sqrt(2)/16 + 1/8)*atan(x/(-1 + sqrt(2))) - 2*(1/8 - sqrt(2)/16)*atan(x/(1 + sqrt(2)))
    // f = 1/(x**4+6*x**2+1): 0.59558823529411765, 0.23139938319198150, 0.088254086897815499, 0.0049021145761435653
    let ctx = Context::new();
    let x = ctx.symbol("x");
    let f = ctx.int(1) / (x.pow(&ctx.int(4)) + x.pow(&ctx.int(2)) * 6 + 1);
    assert_antiderivative(
        &f,
        &x,
        &[
            0.59558823529411765,
            0.23139938319198150,
            0.088254086897815499,
            0.0049021145761435653,
        ],
        "∫ 1/(x⁴+6x²+1)",
    );
}

#[test]
fn biquadratic_with_cubic_numerator() {
    // sympy: integrate((x**3+1)/(x**4+3*x**2+1), x)
    //   = (1/4 - 3*sqrt(5)/20)*log(x**2 - 19737*sqrt(5)/980 - …) + (1/4 + 3*sqrt(5)/20)*log(…) - 2*sqrt(…)*atan(…) - …
    // f = (x**3+1)/(x**4+3*x**2+1): 0.77064220183486239, 0.22772720489479182, 0.37119436819099178, -0.22623108834730347
    let ctx = Context::new();
    let x = ctx.symbol("x");
    let f = (x.pow(&ctx.int(3)) + 1) / (x.pow(&ctx.int(4)) + x.pow(&ctx.int(2)) * 3 + 1);
    assert_antiderivative(
        &f,
        &x,
        &[
            0.77064220183486239,
            0.22772720489479182,
            0.37119436819099178,
            -0.22623108834730347,
        ],
        "∫ (x³+1)/(x⁴+3x²+1)",
    );
}

#[test]
fn one_over_x_fourth_plus_one_unchanged() {
    // Was already right in 0.21 (the Rothstein–Trager route stays in charge
    // when its answer verifies).
    // sympy: integrate(1/(x**4+1), x) = -sqrt(2)*log(x**2 - sqrt(2)*x + 1)/8
    //   + sqrt(2)*log(x**2 + sqrt(2)*x + 1)/8 + sqrt(2)*atan(sqrt(2)*x - 1)/4 + sqrt(2)*atan(sqrt(2)*x + 1)/4
    // f = 1/(x**4+1): 0.98780487804878049, 0.79345670852610707, 0.33889634739132448, 0.0074275663727331067
    let ctx = Context::new();
    let x = ctx.symbol("x");
    let f = ctx.int(1) / (x.pow(&ctx.int(4)) + 1);
    let big_f = assert_antiderivative(
        &f,
        &x,
        &[
            0.98780487804878049,
            0.79345670852610707,
            0.33889634739132448,
            0.0074275663727331067,
        ],
        "∫ 1/(x⁴+1)",
    );
    let s = big_f.to_string();
    assert!(
        s.contains("atan") && s.contains("ln"),
        "expected atan/ln form: {s}"
    );
}

#[test]
fn x_cubed_over_shifted_biquadratic() {
    // sympy: integrate(x**3/(x**4-8*x**2+17), x) = log(x**4 - 8*x**2 + 17)/4 + 2*atan(x**2 - 4)
    // f = x**3/(x**4-8*x**2+17): 0.0022970903522205207, -0.027653119271853865, 0.21224070398538633, -0.67586529466791394
    let ctx = Context::new();
    let x = ctx.symbol("x");
    let f = x.pow(&ctx.int(3)) / (x.pow(&ctx.int(4)) - x.pow(&ctx.int(2)) * 8 + 17);
    assert_antiderivative(
        &f,
        &x,
        &[
            0.0022970903522205207,
            -0.027653119271853865,
            0.21224070398538633,
            -0.67586529466791394,
        ],
        "∫ x³/(x⁴-8x²+17)",
    );
}

#[test]
fn general_quartic_is_never_silently_wrong() {
    // x⁴ + x + 1 is not biquadratic; sympy answers with a RootSum:
    //   integrate(1/(x**4+x+1), x)
    //   = RootSum(229*_t**4 + 18*_t**2 + 8*_t + 1, Lambda(_t, _t*log(2061*_t**3/64 - 687*_t**2/64 + 391*_t/64 + x + 27/64)))
    // f = 1/(x**4+x+1): 0.74311926605504587, 1.8314263920671243, 0.24198000165275597, 0.0076199997561600078
    // The contract: either an unevaluated Integral, or a form whose derivative is f.
    let ctx = Context::new();
    let x = ctx.symbol("x");
    let f = ctx.int(1) / (x.pow(&ctx.int(4)) + &x + 1);
    let big_f = f.integrate(&x);
    if big_f.has_unevaluated() {
        return;
    }
    let d_big_f = big_f.diff(&x);
    let f_ref = [
        0.74311926605504587,
        1.8314263920671243,
        0.24198000165275597,
        0.0076199997561600078,
    ];
    for (&(n, d), &expected) in POINTS.iter().zip(&f_ref) {
        let p = ctx.rational(n, d);
        let df_at = d_big_f
            .subs(&x, &p)
            .eval_f64()
            .unwrap_or_else(|e| panic!("F'({n}/{d}) did not evaluate: {e} (F = {big_f})"));
        close(
            df_at,
            expected,
            1e-10,
            &format!("∫ 1/(x⁴+x+1): F'({n}/{d}), F = {big_f}"),
        );
    }
}

// ═══════════════════════════════════════════════════════════════════════════
// u-substitution must not absorb a stray `x`
// ═══════════════════════════════════════════════════════════════════════════

#[test]
fn u_substitution_rejects_factor_with_bare_var() {
    // With the wider candidate set (u = sin x inside a sum) the factor
    // x + sin²x must not be turned into x + x²: either unevaluated or right.
    // sympy: integrate((x+sin(x)**2)*cos(x), x) = x*sin(x) + sin(x)**3/3 + cos(x)
    // f = (x+sin(x)**2)*cos(x): 0.41614930888322871, -0.21545486337458458, 0.77289471651060822, 3.2239807196321021
    let ctx = Context::new();
    let x = ctx.symbol("x");
    let f = (&x + x.sin().pow(&ctx.int(2))) * x.cos();
    let big_f = f.integrate(&x);
    if big_f.has_unevaluated() {
        return;
    }
    let d_big_f = big_f.diff(&x);
    let f_ref = [
        0.41614930888322871,
        -0.21545486337458458,
        0.77289471651060822,
        3.2239807196321021,
    ];
    for (&(n, d), &expected) in POINTS.iter().zip(&f_ref) {
        let p = ctx.rational(n, d);
        let df_at = d_big_f
            .subs(&x, &p)
            .eval_f64()
            .unwrap_or_else(|e| panic!("F'({n}/{d}) did not evaluate: {e} (F = {big_f})"));
        close(
            df_at,
            expected,
            1e-10,
            &format!("∫ (x + sin²x) cos x: F'({n}/{d}), F = {big_f}"),
        );
    }
}