symplex 0.39.0

//! Rothstein-Trager algorithm for the logarithmic part of rational function
//! integration.
//!
//! Given `∫ A(x)/D(x) dx` where `D` is **square-free** and `deg(A) < deg(D)`,
//! computes the logarithmic part:
//!
//! ```text
//!   ∫ A/D dx  =  Σ cᵢ · ln(vᵢ(x))
//! ```
//!
//! # Algorithm
//!
//! 1. Compute `R(t) = res_x(D, A − t·D')` — a polynomial in `t`.
//! 2. Factor `R(t)` over ℤ: irreducible factors `q` with multiplicities `i`.
//! 3. For each linear factor `(t − c)`, `c ∈ ℚ`:
//!    - `v(x) = gcd(D(x), A(x) − c·D'(x))` — the logarithmic argument.
//!    - Emit `c · ln(v)`.
//! 4. For each irreducible factor `q` of degree > 1 and multiplicity `i`
//!    (Lazard–Rioboo–Trager): the log argument at a root `α` of `q` is
//!    `gcd(D, A − α·D')`, of degree `i` in `x`, and equals `S_i(α, x)` up to
//!    a non-zero constant for the member `S_i` of degree `i` of the
//!    subresultant remainder sequence of `D` and `A − t·D'` in `ℤ[t][x]`
//!    (which also yields `R(t)` in step 1).  Emit an `Algebraic` term
//!    `Σ_{q(α)=0} α·ln(S_i(α, x))`, with `S_i` made monic in `x` and its
//!    coefficients reduced modulo `q`.
//!
//! The distinct roots of `R(t)` are exactly the residues of `A/D`, and the
//! splitting field of `R(t)` is the **minimal** algebraic extension needed
//! to express the integral.  A residue of multiplicity `i` belongs to a log
//! argument of degree `i`: before 0.25 every factor got the degree-1 member,
//! which only exists when all multiplicities are 1, so `∫ x³/(x⁸+1) dx`
//! (`R = (64t²+1)⁴`) fell through to a partial-fraction route over the
//! complex roots whose expansion exhausted memory.
//!
//! # References
//!
//! - Rothstein, "A new algorithm for the integration of exponential and
//!   logarithmic functions", 1977
//! - Trager, "Algebraic factoring and rational function integration", 1976
//! - Bronstein, *Symbolic Integration I*, §2.4–2.5
//! - Lazard & Rioboo, "Integration of rational functions: rational computation
//!   of the logarithmic part", 1990
//! - Bronstein, *Symbolic Integration I*, §2.5 (`IntRationalLogPart`); SymPy's
//!   `ratint_logpart` (`sympy/integrals/rationaltools.py`) is the same
//!   algorithm.

use std::collections::BTreeMap;

use num_bigint::BigInt;
use num_rational::Ratio;
use num_traits::Zero;

use super::LogTerm;
use super::log_to_real::poly_to_genpoly_rf;
use crate::base::numeric::Q;
use crate::poly::dense::Poly;
use crate::poly::generic::GenPoly;
use crate::poly::ratfn::RationalFn;

// ═══════════════════════════════════════════════════════════════════════════
// Public types
// ═══════════════════════════════════════════════════════════════════════════

/// Result of the Rothstein-Trager algorithm.
#[derive(Clone, Debug)]
pub struct LogPartResult {
    /// The logarithmic terms `Σ cᵢ ln(vᵢ)`.
    pub terms: Vec<LogTerm>,
}

// ═══════════════════════════════════════════════════════════════════════════
// Main algorithm
// ═══════════════════════════════════════════════════════════════════════════

/// Compute the logarithmic part of `∫ A/D dx` where `D` is square-free
/// and `deg(A) < deg(D)`.
///
/// If `A` is zero, returns an empty result (no logarithmic part).
///
/// `None` if `D` is zero.  The integrator passes the square-free
/// denominator of a Hermite reduction, which is never zero; up to 0.28 this
/// was a runtime `assert!`.
pub fn logarithmic_part(a: &Poly, d: &Poly) -> Option<LogPartResult> {
    (!d.is_zero()).then(|| logarithmic_part_nonzero(a, d))
}

/// [`logarithmic_part`] for a nonzero `d`.
fn logarithmic_part_nonzero(a: &Poly, d: &Poly) -> LogPartResult {
    if a.is_zero() {
        return LogPartResult { terms: vec![] };
    }

    let d_prime = d.derivative();

    // Handle the trivial case: D is linear (degree 1).
    // Then ∫ A/D dx = (A/lc(D)) · ln(D).  Since deg(A) < deg(D) = 1,
    // A is a constant.
    if d.degree() == Some(1) {
        let a_val = a.coeff(0);
        let d_lc = d.coeff(1);
        let coeff = a_val / d_lc;
        if coeff.is_zero() {
            return LogPartResult { terms: vec![] };
        }
        return LogPartResult {
            terms: vec![LogTerm::Rational {
                coeff,
                argument: d.make_monic(),
            }],
        };
    }

    // R(t) = res_x(D, A − t·D') and the subresultant sequence of D and
    // A − t·D' in one computation, as in Bronstein's IntRationalLogPart.
    // (`Poly::resultant_poly` is the evaluation–interpolation route; its
    // scalar resultants run Euclid over ℚ, 10 s in a debug build for the
    // degree-36 denominator of `∫ atan(√x − x⁹) dx`.)
    let (r_poly, prs) = match lrt_subresultants(a, d, &d_prime) {
        Some((r, prs)) => (r, Some(prs)),
        None => (Poly::resultant_poly(d, a, &d_prime), None),
    };

    if r_poly.is_zero() {
        // Degenerate case: resultant is identically zero.
        // This means D and A - t·D' share a common factor for all t,
        // which shouldn't happen for a proper fraction with square-free D.
        return LogPartResult { terms: vec![] };
    }

    // Factor R(t) over ℤ.  A factor's multiplicity is the degree in x of
    // the log argument belonging to its roots.
    let (_content, r_factors) = r_poly.factor_over_z();

    let mut terms: Vec<LogTerm> = Vec::new();

    for (factor, multiplicity) in &r_factors {
        let deg = match factor.degree() {
            Some(d) => d,
            None => continue, // zero polynomial — skip
        };

        if deg == 0 {
            // Constant factor — no roots, skip.
            continue;
        }

        if deg == 1 {
            // Linear factor: (t - c) where c = -factor[0] / factor[1].
            let c = extract_linear_root(factor);

            // v(x) = gcd(D(x), A(x) − c · D'(x))
            let a_minus_c_dprime = a - &d_prime.scale(&c);
            let v = Poly::gcd(d, &a_minus_c_dprime);

            if v.degree().unwrap_or(0) == 0 {
                // GCD is constant — this root doesn't contribute a log term.
                continue;
            }

            if !c.is_zero() {
                terms.push(LogTerm::Rational {
                    coeff: c,
                    argument: v.make_monic(),
                });
            }
        } else {
            // Irreducible factor of degree > 1: roots are algebraic numbers
            // not in ℚ, with log argument S_i(α, x), i = multiplicity.
            let i = *multiplicity as usize;
            let log_arg = if d.degree() == Some(i) {
                Some(poly_to_genpoly_rf(&d.make_monic()))
            } else {
                prs.as_ref()
                    .and_then(|p| p.get(&i))
                    .and_then(|s| reduce_log_arg(s, factor))
            };
            if log_arg.is_none() {
                tracing::debug!(
                    min_poly_degree = deg,
                    multiplicity = i,
                    "logarithmic_part: no usable remainder of degree i"
                );
            }
            terms.push(LogTerm::Algebraic {
                min_poly: factor.clone(),
                log_arg,
            });
        }
    }

    LogPartResult { terms }
}

/// `R(t) = res_x(D, A − t·D')` (up to a non-zero constant factor) and the
/// subresultant remainder sequence of `D` and `A − t·D'`, keyed by degree
/// in `x`, both computed in `ℤ[t][x]` after scaling `D` and `A − t·D'` to
/// integer coefficients; each member as an element of `ℚ(t)[x]` with
/// polynomial coefficients.
///
/// The member of degree `i` is a multiple of the subresultant `S_i` by a
/// factor that does not vanish at a root `α` of a factor of multiplicity
/// `i` of the resultant, so it specialises there to a non-zero multiple of
/// `gcd(D, A − α·D')`, which is all [`reduce_log_arg`] needs (Lazard–Rioboo;
/// Bronstein, *Symbolic Integration I*, §2.5, `IntRationalLogPart`, uses
/// the same subresultant PRS for both).  Up to 0.28.0 the resultant came from
/// evaluation–interpolation and the sequence was Euclid's over `ℚ(t)` with
/// every member made primitive: each coefficient operation there
/// normalises a rational function by a gcd in `ℚ[t]`, and for the
/// degree-36 denominator of `∫ atan(√x − x⁹) dx` the sequence took 24 s of
/// a 28 s integral (release).  `None` if the sequence fails (an inexact
/// division, which theory rules out).
fn lrt_subresultants(
    a: &Poly,
    d: &Poly,
    d_prime: &Poly,
) -> Option<(Poly, BTreeMap<usize, GenPoly<RationalFn>>)> {
    use crate::poly::zpoly::{denominator_lcm, integer_scaled, ztx_subresultant_prs};
    let constant = |c: BigInt| if c.is_zero() { Vec::new() } else { vec![c] };
    let d_z: Vec<Vec<BigInt>> = integer_scaled(d.coeffs())
        .into_iter()
        .map(constant)
        .collect();
    // B = A − t·D' over the common denominator of A and D'.
    let den = denominator_lcm(a.coeffs().iter().chain(d_prime.coeffs()));
    let len = a.coeffs().len().max(d_prime.coeffs().len());
    let b_z: Vec<Vec<BigInt>> = (0..len)
        .map(|k| {
            let mut c = vec![
                (a.coeff(k) * &den).to_integer(),
                -(d_prime.coeff(k) * &den).to_integer(),
            ];
            crate::poly::zpoly::z_normalize(&mut c);
            c
        })
        .collect();
    let chain = ztx_subresultant_prs(&d_z, &b_z)?;
    let to_poly =
        |c: Vec<BigInt>| Poly::from_coeffs(c.into_iter().map(Ratio::from_integer).collect());
    let members = chain
        .members
        .into_iter()
        .map(|(deg, member)| {
            let coeffs = member
                .into_iter()
                .map(|c| RationalFn::from_poly(to_poly(c)))
                .collect();
            (deg, GenPoly::from_coeffs(coeffs))
        })
        .collect();
    Some((to_poly(chain.resultant), members))
}

/// `p ∈ ℚ(t)[x]` scaled by a non-zero element of ℚ(t) so that its
/// coefficients are polynomials in `t` with no common factor: clear the
/// denominators (their lcm in `ℚ[t]`) and divide by the gcd of the
/// numerators.  The normalisation of the Euclidean sequence `lrt_prs` used
/// up to 0.28.0, kept as the reference the subresultant sequence is tested
/// against.
#[cfg(test)]
fn primitive_over_q_t(p: GenPoly<RationalFn>) -> GenPoly<RationalFn> {
    let coeffs = &p.coeffs;
    let mut lcm = Poly::from_int(1);
    for c in coeffs {
        let d = c.denom();
        if !d.is_constant() {
            let g = Poly::gcd(&lcm, d);
            lcm = (&lcm * d).div(&g);
        }
    }
    let numers: Vec<Poly> = coeffs
        .iter()
        .map(|c| (c.numer() * &lcm).div(c.denom()))
        .collect();
    let content = numers
        .iter()
        .filter(|n| !n.is_zero())
        .fold(Poly::zero(), |g, n| Poly::gcd(&g, n));
    if content.is_zero() {
        return p;
    }
    GenPoly::from_coeffs(
        numers
            .iter()
            .map(|n| RationalFn::from_poly(n.div(&content)))
            .collect(),
    )
}

/// `S(t, x)` made monic in `x`, each coefficient reduced to a polynomial in
/// `t` of degree below `deg q` — an element of `ℚ[t]/(q)`, the field of
/// the roots of the irreducible `q`.  `None` if the leading coefficient
/// vanishes at the roots of `q` (for a primitive remainder of the right
/// degree it cannot: the specialisation is then a non-zero multiple of
/// the gcd, which has that degree).
fn reduce_log_arg(s: &GenPoly<RationalFn>, q: &Poly) -> Option<GenPoly<RationalFn>> {
    let lc = s.leading_coeff()?.clone();
    let mut coeffs = Vec::with_capacity(s.coeffs.len());
    for c in &s.coeffs {
        // c / lc = (c.numer · lc.denom) / (c.denom · lc.numer)
        let num = (c.numer() * lc.denom()).rem(q);
        let den = (c.denom() * lc.numer()).rem(q);
        let inv = inverse_mod(&den, q)?;
        coeffs.push(RationalFn::from_poly((&num * &inv).rem(q)));
    }
    Some(GenPoly::from_coeffs(coeffs))
}

/// The inverse of `a` modulo `q`, if `gcd(a, q) = 1`.
fn inverse_mod(a: &Poly, q: &Poly) -> Option<Poly> {
    if a.is_zero() {
        return None;
    }
    let eg = Poly::extended_gcd(a, q);
    if eg.gcd.degree() != Some(0) {
        return None;
    }
    // eg.gcd is monic, i.e. 1: eg.x · a ≡ 1 (mod q).
    Some(eg.x.rem(q))
}

// ═══════════════════════════════════════════════════════════════════════════
// Helpers
// ═══════════════════════════════════════════════════════════════════════════

/// Extract the rational root of a monic-ish linear polynomial `a·t + b`.
/// Returns `c = -b/a`.
fn extract_linear_root(factor: &Poly) -> Q {
    let a = factor.coeff(1);
    let b = factor.coeff(0);
    if a.is_zero() { Ratio::zero() } else { -b / a }
}

// ═══════════════════════════════════════════════════════════════════════════
// Tests
// ═══════════════════════════════════════════════════════════════════════════

#[cfg(test)]
mod tests {
    use super::*;
    use num_bigint::BigInt;
    use num_traits::Signed;

    fn rat(n: i64, d: i64) -> Q {
        Ratio::new(BigInt::from(n), BigInt::from(d))
    }

    /// Verify a logarithmic part result by differentiating and comparing
    /// to the original integrand A/D.
    ///
    /// For Σ cᵢ ln(vᵢ), the derivative is Σ cᵢ · vᵢ'/vᵢ.
    /// This should equal A/D (as rational functions).
    fn verify_log_part(a: &Poly, d: &Poly, result: &LogPartResult) {
        // Compute Σ cᵢ · vᵢ' / vᵢ  as a single fraction.
        let mut sum_numer = Poly::zero();
        let mut sum_denom = Poly::from_int(1);

        for term in &result.terms {
            match term {
                LogTerm::Rational { coeff, argument } => {
                    let v = argument;
                    let v_prime = v.derivative();
                    // Add c · v'/v to the running sum.
                    // sum = sum_numer/sum_denom + c·v'/(v)
                    //     = (sum_numer · v + c · v' · sum_denom) / (sum_denom · v)
                    let c_poly = Poly::constant(coeff.clone());
                    sum_numer = &(&sum_numer * v) + &(&(&c_poly * &v_prime) * &sum_denom);
                    sum_denom = &sum_denom * v;
                }
                LogTerm::Algebraic { .. } => {
                    // Can't verify algebraic terms numerically at the Poly level.
                    // Skip verification for these.
                    return;
                }
            }
        }

        // Now check: sum_numer/sum_denom == A/D
        // i.e., sum_numer * D == A * sum_denom
        let lhs = &sum_numer * d;
        let rhs = &(a * &sum_denom);
        let diff = &lhs - rhs;

        assert!(
            diff.is_zero(),
            "Rothstein-Trager verification failed:\n  d/dx(Σ cᵢ ln(vᵢ)) ≠ A/D\n  diff = {diff}"
        );
    }

    #[test]
    fn log_part_one_over_x() {
        // ∫ 1/x dx = ln(x)
        let a = Poly::from_int(1);
        let d = Poly::x(); // x

        let result = logarithmic_part(&a, &d).unwrap();
        assert_eq!(result.terms.len(), 1, "should have one log term");
        match &result.terms[0] {
            LogTerm::Rational { coeff, argument } => {
                assert_eq!(*coeff, rat(1, 1), "coefficient should be 1");
                assert_eq!(argument.degree(), Some(1), "argument should be linear (x)");
            }
            _ => panic!("expected rational log term"),
        }
        verify_log_part(&a, &d, &result);
    }

    #[test]
    fn log_part_one_over_x_minus_1() {
        // ∫ 1/(x-1) dx = ln(x-1)
        let a = Poly::from_int(1);
        let d = Poly::from_coeffs(vec![rat(-1, 1), rat(1, 1)]); // x - 1

        let result = logarithmic_part(&a, &d).unwrap();
        assert_eq!(result.terms.len(), 1);
        match &result.terms[0] {
            LogTerm::Rational { coeff, argument } => {
                assert_eq!(*coeff, rat(1, 1));
                // argument should be (x - 1) or equivalent monic form
                assert_eq!(argument.degree(), Some(1));
            }
            _ => panic!("expected rational log term"),
        }
        verify_log_part(&a, &d, &result);
    }

    #[test]
    fn log_part_one_over_x_squared_minus_1() {
        // ∫ 1/(x^2 - 1) dx = 1/2 · ln(x - 1) - 1/2 · ln(x + 1)
        // D = x^2 - 1 = (x-1)(x+1), already square-free.
        let a = Poly::from_int(1);
        let d = Poly::from_coeffs(vec![rat(-1, 1), rat(0, 1), rat(1, 1)]); // x^2 - 1

        let result = logarithmic_part(&a, &d).unwrap();
        // Should have 2 rational log terms with coefficients ±1/2.
        let rational_count = result
            .terms
            .iter()
            .filter(|t| matches!(t, LogTerm::Rational { .. }))
            .count();
        assert!(
            rational_count >= 1,
            "should have rational log terms for 1/(x^2-1), got {} rational terms out of {} total",
            rational_count,
            result.terms.len()
        );
        verify_log_part(&a, &d, &result);
    }

    #[test]
    fn log_part_2x_over_x_squared_plus_1() {
        // ∫ 2x/(x^2 + 1) dx = ln(x^2 + 1)
        // A = 2x, D = x^2 + 1 (square-free)
        let a = Poly::from_coeffs(vec![rat(0, 1), rat(2, 1)]); // 2x
        let d = Poly::from_coeffs(vec![rat(1, 1), rat(0, 1), rat(1, 1)]); // x^2 + 1

        let result = logarithmic_part(&a, &d).unwrap();
        // Should produce 1 · ln(x^2 + 1)
        assert!(
            !result.terms.is_empty(),
            "should have at least one log term for 2x/(x^2+1)"
        );
        verify_log_part(&a, &d, &result);
    }

    #[test]
    fn log_part_one_over_x_squared_plus_1() {
        // ∫ 1/(x^2 + 1) dx = arctan(x)
        // But in the Risch framework, this is expressed as:
        //   (i/2)·ln(x - i) - (i/2)·ln(x + i)
        // The resultant R(t) = 4t^2 + 1, which is irreducible over ℤ.
        // So we should get an Algebraic term.
        let a = Poly::from_int(1);
        let d = Poly::from_coeffs(vec![rat(1, 1), rat(0, 1), rat(1, 1)]); // x^2 + 1

        let result = logarithmic_part(&a, &d).unwrap();
        // Should have an algebraic term (roots are ±i/2).
        let has_algebraic = result
            .terms
            .iter()
            .any(|t| matches!(t, LogTerm::Algebraic { .. }));
        assert!(
            has_algebraic || result.terms.is_empty(),
            "1/(x^2+1) should produce algebraic log terms or be handled as arctan"
        );
    }

    #[test]
    fn log_part_zero_numerator() {
        let a = Poly::zero();
        let d = Poly::from_coeffs(vec![rat(-1, 1), rat(0, 1), rat(1, 1)]); // x^2 - 1
        let result = logarithmic_part(&a, &d).unwrap();
        assert!(
            result.terms.is_empty(),
            "zero numerator should give empty result"
        );
    }

    #[test]
    fn log_part_partial_fraction_style() {
        // ∫ (3x + 5) / ((x + 1)(x + 2)) dx = ∫ (3x+5)/(x^2+3x+2) dx
        // = 2·ln(x+1) + 1·ln(x+2)    [partial fractions: 2/(x+1) + 1/(x+2)]
        let a = Poly::from_coeffs(vec![rat(5, 1), rat(3, 1)]); // 3x + 5
        let d = Poly::from_coeffs(vec![rat(2, 1), rat(3, 1), rat(1, 1)]); // x^2 + 3x + 2

        let result = logarithmic_part(&a, &d).unwrap();
        verify_log_part(&a, &d, &result);

        // Should have exactly 2 rational log terms.
        let rational_terms: Vec<_> = result
            .terms
            .iter()
            .filter_map(|t| match t {
                LogTerm::Rational { coeff, argument } => Some((coeff.clone(), argument.clone())),
                _ => None,
            })
            .collect();
        assert_eq!(
            rational_terms.len(),
            2,
            "should have 2 rational log terms, got {}: {:?}",
            rational_terms.len(),
            rational_terms
        );

        // Check coefficients sum to 3 (= leading coeff of A divided by leading coeff of D)
        // Actually, the coefficients should be 2 and 1.
        let mut coeffs: Vec<Q> = rational_terms.iter().map(|(c, _)| c.clone()).collect();
        coeffs.sort_by(|a, b| a.partial_cmp(b).unwrap());
        assert_eq!(
            coeffs,
            vec![rat(1, 1), rat(2, 1)],
            "coefficients should be [1, 2], got {coeffs:?}"
        );
    }

    #[test]
    fn log_part_one_over_x_cubed_minus_1() {
        // ∫ 1/(x^3 - 1) dx
        // D = x^3 - 1 = (x - 1)(x^2 + x + 1), square-free.
        // The (x-1) factor gives a rational log term.
        // The (x^2+x+1) factor gives algebraic log terms (cube roots of unity).
        let a = Poly::from_int(1);
        let d = Poly::from_coeffs(vec![rat(-1, 1), rat(0, 1), rat(0, 1), rat(1, 1)]); // x^3 - 1

        let result = logarithmic_part(&a, &d).unwrap();
        // Should have at least one rational term (from the x-1 factor)
        // and possibly algebraic terms.
        assert!(!result.terms.is_empty(), "1/(x^3-1) should have log terms");
        // Only verify if all terms are rational (can't verify algebraic).
        let all_rational = result
            .terms
            .iter()
            .all(|t| matches!(t, LogTerm::Rational { .. }));
        if all_rational {
            verify_log_part(&a, &d, &result);
        }
    }

    #[test]
    fn log_part_x_over_x_squared_minus_1() {
        // ∫ x/(x^2 - 1) dx = 1/2 · ln(x^2 - 1) = 1/2·ln(x-1) + 1/2·ln(x+1)
        let a = Poly::x(); // x
        let d = Poly::from_coeffs(vec![rat(-1, 1), rat(0, 1), rat(1, 1)]); // x^2 - 1

        let result = logarithmic_part(&a, &d).unwrap();
        verify_log_part(&a, &d, &result);

        // Both log terms should have coefficient 1/2.
        let rational_terms: Vec<_> = result
            .terms
            .iter()
            .filter_map(|t| match t {
                LogTerm::Rational { coeff, .. } => Some(coeff.clone()),
                _ => None,
            })
            .collect();
        for c in &rational_terms {
            assert_eq!(c.abs(), rat(1, 2), "coefficients should be ±1/2, got {c}");
        }
    }

    /// `S(t, x)` divides `D(x)` in `(ℚ[t]/q)[x]`: every root `α` of `q` makes
    /// `S(α, x)` a factor of `D`.
    fn assert_log_arg_divides(d: &Poly, q: &Poly, s: &GenPoly<RationalFn>) {
        let d_gp = poly_to_genpoly_rf(d);
        let (_, rem) = d_gp.div_rem(s);
        for c in &rem.coeffs {
            assert!(c.denom().is_constant(), "remainder coefficient {c:?}");
            assert!(
                c.numer().rem(q).is_zero(),
                "S does not divide D modulo q: remainder coefficient {c:?}"
            );
        }
    }

    fn algebraic_terms(result: &LogPartResult) -> Vec<(Poly, GenPoly<RationalFn>)> {
        result
            .terms
            .iter()
            .filter_map(|t| match t {
                LogTerm::Algebraic { min_poly, log_arg } => {
                    Some((min_poly.clone(), log_arg.clone().expect("log argument")))
                }
                LogTerm::Rational { .. } => None,
            })
            .collect()
    }

    fn poly(coeffs: &[i64]) -> Poly {
        Poly::from_coeffs(coeffs.iter().map(|&c| rat(c, 1)).collect())
    }

    #[test]
    fn log_argument_has_the_degree_of_the_multiplicity() {
        // ∫ x³/(x⁸ + 1): R(t) = (64t² + 1)⁴ (SymPy 1.14:
        // factor_list(resultant(x**8+1, x**3 - t*diff(x**8+1, x), x), t)),
        // so each residue ±i/8 has a log argument of degree 4:
        // gcd(x⁸ + 1, x³ − α·8x⁷) = x⁴ + 8α.
        let a = poly(&[0, 0, 0, 1]);
        let d = poly(&[1, 0, 0, 0, 0, 0, 0, 0, 1]);
        let result = logarithmic_part(&a, &d).unwrap();
        let alg = algebraic_terms(&result);
        assert_eq!(alg.len(), 1, "{:?}", result.terms);
        let (q, s) = &alg[0];
        assert_eq!(q.primitive_part(), poly(&[1, 0, 64]));
        assert_eq!(s.degree(), Some(4));
        assert_eq!(s.coeff(0).numer(), &poly(&[0, 8]));
        for k in 1..4 {
            assert!(s.coeff(k).numer().is_zero());
        }
        assert_log_arg_divides(&d, q, s);
    }

    #[test]
    fn log_arguments_of_mixed_factors() {
        // ∫ x/(x⁶ + 1): R(t) = (6t − 1)²·(36t² + 6t + 1)² (SymPy 1.14, as
        // above): a rational residue 1/6 with log argument x² + 1, and a
        // quadratic factor with log arguments of degree 2.
        let a = poly(&[0, 1]);
        let d = poly(&[1, 0, 0, 0, 0, 0, 1]);
        let result = logarithmic_part(&a, &d).unwrap();
        let rational: Vec<_> = result
            .terms
            .iter()
            .filter_map(|t| match t {
                LogTerm::Rational { coeff, argument } => Some((coeff.clone(), argument.clone())),
                LogTerm::Algebraic { .. } => None,
            })
            .collect();
        assert_eq!(rational, vec![(rat(1, 6), poly(&[1, 0, 1]))]);
        let alg = algebraic_terms(&result);
        assert_eq!(alg.len(), 1);
        let (q, s) = &alg[0];
        assert_eq!(q.primitive_part(), poly(&[1, 6, 36]));
        assert_eq!(s.degree(), Some(2));
        assert_log_arg_divides(&d, q, s);
    }

    #[test]
    fn log_argument_of_a_simple_quadratic_factor() {
        // ∫ 1/(x² + x + 1): q = 3t² + 1, S = x + 3t/2 + 1/2 (the example in
        // SymPy's `ratint_logpart` docstring).
        let a = poly(&[1]);
        let d = poly(&[1, 1, 1]);
        let alg = algebraic_terms(&logarithmic_part(&a, &d).unwrap());
        assert_eq!(alg.len(), 1);
        let (q, s) = &alg[0];
        assert_eq!(q.primitive_part(), poly(&[1, 0, 3]));
        assert_eq!(s.degree(), Some(1));
        assert_eq!(
            s.coeff(0).numer(),
            &Poly::from_coeffs(vec![rat(1, 2), rat(3, 2)])
        );
        assert_log_arg_divides(&d, q, s);
    }

    #[test]
    fn extract_linear_root_basic() {
        // t - 3  →  root = 3
        let p = Poly::from_coeffs(vec![rat(-3, 1), rat(1, 1)]);
        assert_eq!(extract_linear_root(&p), rat(3, 1));
    }

    #[test]
    fn extract_linear_root_scaled() {
        // 2t + 6  →  root = -3
        let p = Poly::from_coeffs(vec![rat(6, 1), rat(2, 1)]);
        assert_eq!(extract_linear_root(&p), rat(-3, 1));
    }

    #[test]
    fn extract_linear_root_fractional() {
        // 3t - 1  →  root = 1/3
        let p = Poly::from_coeffs(vec![rat(-1, 1), rat(3, 1)]);
        assert_eq!(extract_linear_root(&p), rat(1, 3));
    }

    /// The subresultant chain gives exactly the resultant factors and the
    /// reduced log arguments of the route used up to 0.28.0: the resultant by
    /// evaluation–interpolation and Euclid's sequence over `ℚ(t)` with
    /// primitive members.  Covers multiplicities 1–4, a cubic factor, and a
    /// few seeded random denominators.
    #[test]
    fn subresultant_log_arguments_match_the_euclidean_route() {
        use super::super::log_to_real::poly_to_genpoly_rf_times_t;
        use crate::base::rng::SplitMix64;
        let p = |c: &[i64]| Poly::from_coeffs(c.iter().map(|&k| rat(k, 1)).collect());
        let as_pairs = |s: &GenPoly<RationalFn>| -> Vec<(Poly, Poly)> {
            s.coeffs
                .iter()
                .map(|c| (c.numer().clone(), c.denom().clone()))
                .collect()
        };
        let mut cases: Vec<(Poly, Poly)> = vec![
            (p(&[1]), p(&[-2, 0, 0, 1])),                        // 1/(x³ − 2)
            (p(&[0, 0, 0, 1]), p(&[1, 0, 0, 0, 0, 0, 0, 0, 1])), // x³/(x⁸ + 1)
            (p(&[1]), p(&[1, 0, 0, 0, 1])),                      // 1/(x⁴ + 1)
            (p(&[1, 0, 0, 0, 1]), p(&[1, 0, 0, 0, 0, 0, 1])),    // (x⁴ + 1)/(x⁶ + 1)
            (p(&[0, 1]), p(&[1, 0, 0, 0, 1, 0, 0, 0, 1])),       // x/(x⁸ + x⁴ + 1)
            (p(&[0, 0, 1]), p(&[1, 0, 1, 0, 0, -2, 0, 0, 1])),   // x²/((x − x⁴)² + 1)
        ];
        let mut rng = SplitMix64::new(20260924);
        while cases.len() < 16 {
            let mut rand = |deg: usize| {
                let mut c: Vec<i64> = (0..deg).map(|_| (rng.next_u64() % 7) as i64 - 3).collect();
                c.push(1);
                p(&c)
            };
            let d = rand(3 + cases.len() % 5);
            let a = rand(d.degree().unwrap_or(1) - 1);
            if Poly::gcd(&d, &d.derivative()).degree() == Some(0)
                && Poly::gcd(&a, &d).degree() == Some(0)
            {
                cases.push((a, d));
            }
        }
        for (a, d) in &cases {
            let d_prime = d.derivative();
            let (r_new, prs_new) = lrt_subresultants(a, d, &d_prime).unwrap();
            let r_old = Poly::resultant_poly(d, a, &d_prime);
            let (_, factors) = r_old.factor_over_z();
            assert_eq!(r_new.factor_over_z().1, factors, "R(t) for {a:?} / {d:?}");
            let d_gp = poly_to_genpoly_rf(d);
            let b_gp = &poly_to_genpoly_rf(a) - &poly_to_genpoly_rf_times_t(&d_prime);
            let prs_old = GenPoly::<RationalFn>::euclidean_prs_normalized(
                &d_gp,
                &b_gp,
                None,
                primitive_over_q_t,
            );
            for (q, i) in factors.iter().filter(|(q, _)| q.degree() >= Some(2)) {
                let i = *i as usize;
                if d.degree() == Some(i) {
                    continue;
                }
                let new = prs_new.get(&i).and_then(|s| reduce_log_arg(s, q));
                let old = prs_old.get(&i).and_then(|s| reduce_log_arg(s, q));
                assert!(
                    new.is_some(),
                    "no log argument of degree {i} for {a:?} / {d:?}"
                );
                assert_eq!(
                    new.as_ref().map(as_pairs),
                    old.as_ref().map(as_pairs),
                    "log argument of degree {i} for {a:?} / {d:?}"
                );
            }
        }
    }
}