symplex 0.57.0

//! Products and quotients of Γ at real arguments, for the functions `evalf`
//! had no kernel for (or refused at the poles of Γ): the rising and falling
//! factorials, the binomial coefficient at the poles, harmonic numbers of a
//! non-integer argument, and polygamma functions of large order; and the
//! rising and falling factorials of non-real arguments ([`eval_complex`],
//! since 0.56).
//!
//! * `rf(x, n) = Γ(x + n)/Γ(x)`, `ff(x, n) = Γ(x + 1)/Γ(x − n + 1)` and
//!   `C(n, k) = Γ(n + 1)/(Γ(k + 1)·Γ(n − k + 1))` are quotients of Γ
//!   ([`gamma_ratio`], mpmath's `gammaprod`, BSD, see THIRD-PARTY-NOTICES.md):
//!   a pole of Γ more in the denominator makes the value 0, one more in the
//!   numerator infinite, and poles that pair up cancel through
//!   `lim Γ(−i + ε)/Γ(−j + ε) = (−1)^{i+j}·Γ(1 + j)/Γ(1 + i)` — so
//!   `C(−7, 2500) = (−1)^{2500}·C(2506, 6)`.  The quotient is taken in log
//!   space, `Σ ± ln|Γ|` with the signs apart, at enough bits for the
//!   logarithms' absolute error (a huge `n` makes `ln Γ(x + n)` large and
//!   the quotient moderate: `rf(1/3, 10⁹)`).
//! * `H(z) = ψ(z + 1) + γ` ([`harmonic`]), with the cancellation near
//!   `z = 0` measured and made up.
//! * `ψ⁽ⁿ⁾(x) = (−1)^{n+1}·n!·ζ(n + 1, x)` for `n > 10⁴` ([`polygamma_large`]):
//!   `ln n!` plus the logarithm of the Hurwitz zeta function, scaled by its
//!   largest term `d^{−s}` (`d` the distance from `x` to the nearest pole
//!   among the terms, `s = n + 1`), summed directly while the terms matter
//!   (`(d/(x + k))^s` falls off at once for a large `s`) with a rigorous
//!   bound on the rest, `u^{−s}·(1 + u/(s − 1))`, or with the
//!   Euler–Maclaurin tail once `x + k ≥ max(s/π, wp)`.  Left of 0 an odd
//!   `s` gives terms of both signs; they are summed in pairs `(f + j,
//!   1 − f + j)` (`f` the fractional part), which cancel exactly at a
//!   half-integer `x`.

use super::ball::Cache;
use super::wp::{RM, Wp};
use crate::base::bigcomplex::{Complex, c_add, c_div, c_mul, c_sub};
use crate::base::bigfloat::is_strictly_positive;
use crate::base::libfn::LibFn;
use crate::base::node::ExprId;
use astro_float::BigFloat;

use crate::base::errors::SymplexError;

fn unevaluable(reason: impl Into<String>) -> SymplexError {
    SymplexError::Unevaluable {
        reason: reason.into(),
    }
}

/// `x > 0` strictly (astro-float's `is_positive` is true for `+0`).
fn pos(x: &BigFloat) -> bool {
    is_strictly_positive(x)
}

/// Is `x` a pole of Γ, a non-positive integer?
fn is_npint(x: &BigFloat) -> bool {
    x.is_zero() || (x.is_int() && x.is_negative())
}

/// Is the integer `x` even?
fn is_even(x: &BigFloat) -> bool {
    if x.is_zero() {
        return true;
    }
    let mut half = x.clone();
    match x.exponent() {
        Some(e) => {
            half.set_exponent(e - 1);
            half.is_int()
        }
        None => true,
    }
}

/// The sign of `Γ(x)` at a real `x` off the poles: `+` right of 0,
/// `(−1)^{⌈−x⌉}` left of it (negative on `(−1, 0)`, positive on `(−2, −1)`).
fn gamma_sign(x: &BigFloat) -> i32 {
    if pos(x) || is_even(&x.floor()) { 1 } else { -1 }
}

/// `a + b` exactly (the precision spans both operands' bits), or to
/// `2¹⁸ + prec` bits when they are further apart than that.
pub(super) fn exact_sum(a: &BigFloat, b: &BigFloat, prec: usize) -> BigFloat {
    let span = |x: &BigFloat| -> Option<(i64, i64)> {
        let e = i64::from(x.exponent()?);
        let m = x.mantissa_max_bit_len().unwrap_or(64) as i64;
        (!x.is_zero()).then_some((e, e - m))
    };
    let p = match (span(a), span(b)) {
        (Some((ea, la)), Some((eb, lb))) => {
            let bits = ea.max(eb) + 2 - la.min(lb);
            usize::try_from(bits).unwrap_or(prec).min(prec + (1 << 18))
        }
        (Some((ea, la)), None) | (None, Some((ea, la))) => {
            usize::try_from(ea - la + 2).unwrap_or(prec)
        }
        (None, None) => prec,
    }
    .max(64);
    a.add(b, p, RM)
}

fn int(n: i32) -> BigFloat {
    BigFloat::from_i32(n, 64)
}

/// `Π Γ(aᵢ) / Π Γ(bⱼ)` at real arguments (mpmath's `gammaprod`): 0 when
/// the denominator has more poles than the numerator, infinite (an error)
/// when the numerator has more, and otherwise the limit, the poles paired
/// by `lim Γ(i + ε)/Γ(j + ε) = (−1)^{i+j}·Γ(1 − j)/Γ(1 − i)`.  Computed as
/// `±exp(Σ ln|Γ(aᵢ)| − Σ ln|Γ(bⱼ)|)` with the logarithms to an absolute
/// error below `2^{−prec−32}` (see the module documentation).
///
/// # Errors
///
/// [`SymplexError::Unevaluable`] for an infinite value, or one beyond the
/// floating-point exponent range.
pub(super) fn gamma_ratio(
    num: &[BigFloat],
    den: &[BigFloat],
    w: &mut Wp<'_>,
) -> Result<BigFloat, SymplexError> {
    let prec = w.bits();
    if num.iter().chain(den).any(|x| x.is_nan() || x.is_inf()) {
        return Err(unevaluable("Gamma quotient of a non-finite argument"));
    }
    let (poles_num, mut regular_num): (Vec<&BigFloat>, Vec<BigFloat>) = split_poles(num);
    let (poles_den, mut regular_den): (Vec<&BigFloat>, Vec<BigFloat>) = split_poles(den);
    if poles_num.len() < poles_den.len() {
        return Ok(BigFloat::new(prec));
    }
    if poles_num.len() > poles_den.len() {
        return Err(unevaluable(
            "the Gamma quotient has a pole here (more poles of Gamma above than below): the value is infinite",
        ));
    }
    let mut sign = 1;
    for (i, j) in poles_num.into_iter().zip(poles_den) {
        if !is_even(&exact_sum(i, j, prec)) {
            sign = -sign;
        }
        regular_num.push(exact_sum(&int(1), &j.neg(), prec));
        regular_den.push(exact_sum(&int(1), &i.neg(), prec));
    }
    // Γ(1) = Γ(2) = 1: such factors drop out.  (Their logarithm is exactly
    // 0, which `arb_log_gamma` cannot certify to any relative precision: up
    // to 0.39 `rf(1/2, 1/2) = 1/√π`, `ff(3/2, 1/2)` and `C(−1, −1) = 1` were
    // refused with `PrecisionExhausted`.)
    let unit = |x: &BigFloat| x.cmp(&int(1)) == Some(0) || x.cmp(&int(2)) == Some(0);
    regular_num.retain(|x| !unit(x));
    regular_den.retain(|x| !unit(x));
    for x in regular_num.iter().chain(&regular_den) {
        sign *= gamma_sign(x);
    }
    // A first pass for the logarithms' size, a second with their bits.
    let mut wp = prec + 64;
    loop {
        let mut total = BigFloat::new(wp);
        let mut size = 0.0f64;
        for (x, plus) in regular_num
            .iter()
            .map(|x| (x, true))
            .chain(regular_den.iter().map(|x| (x, false)))
        {
            let l = super::arb_log_gamma(x, &mut w.at(wp))?.0;
            size += super::f64_of(&l).unwrap_or(f64::INFINITY).abs();
            total = if plus {
                total.add(&l, wp, RM)
            } else {
                total.sub(&l, wp, RM)
            };
        }
        let needed = prec + 32 + (size + 1.0).log2().ceil().max(0.0) as usize;
        if needed <= wp || !size.is_finite() {
            let v = total.exp(wp, RM, w.cc);
            if v.is_inf() || v.is_nan() || (v.is_zero() && !total.is_zero()) {
                return Err(unevaluable(
                    "the Gamma quotient is beyond the floating-point exponent range",
                ));
            }
            let v = super::round_to(v, prec);
            return Ok(if sign < 0 { v.neg() } else { v });
        }
        wp = needed;
    }
}

fn split_poles(xs: &[BigFloat]) -> (Vec<&BigFloat>, Vec<BigFloat>) {
    let mut poles = Vec::new();
    let mut regular = Vec::new();
    for x in xs {
        if is_npint(x) {
            poles.push(x);
        } else {
            regular.push(x.clone());
        }
    }
    (poles, regular)
}

/// `rf(x, n) = Γ(x + n)/Γ(x)` (the Pochhammer symbol) for real `x`, `n`.
pub(super) fn rising_factorial(
    x: &BigFloat,
    n: &BigFloat,
    w: &mut Wp<'_>,
) -> Result<BigFloat, SymplexError> {
    let prec = w.bits();
    if n.is_zero() {
        return Ok(BigFloat::from_i32(1, prec));
    }
    let xn = exact_sum(x, n, prec);
    gamma_ratio(&[xn], std::slice::from_ref(x), w)
}

/// `ff(x, n) = Γ(x + 1)/Γ(x − n + 1)` for real `x`, `n`.
pub(super) fn falling_factorial(
    x: &BigFloat,
    n: &BigFloat,
    w: &mut Wp<'_>,
) -> Result<BigFloat, SymplexError> {
    let prec = w.bits();
    if n.is_zero() {
        return Ok(BigFloat::from_i32(1, prec));
    }
    let x1 = exact_sum(x, &int(1), prec);
    let xn1 = exact_sum(&x1, &n.neg(), prec);
    gamma_ratio(&[x1], &[xn1], w)
}

/// Does `f(args)` go to [`eval_complex`]: `rising_factorial` or
/// `falling_factorial` with a non-real (finite) argument?
pub(super) fn nonreal_args(f: LibFn, args: &[ExprId], cache: Cache<'_>) -> bool {
    if !matches!(f, LibFn::RisingFactorial | LibFn::FallingFactorial) {
        return false;
    }
    let vals: Option<Vec<&Complex>> = args.iter().map(|a| cache.get(a)).collect();
    vals.is_some_and(|vals| {
        vals.iter().any(|v| !v.1.is_zero())
            && vals
                .iter()
                .all(|v| [&v.0, &v.1].iter().all(|x| !x.is_nan() && !x.is_inf()))
    })
}

/// `rf(x, n) = Γ(x + n)/Γ(x)` and `ff(x, n) = Γ(x + 1)/Γ(x − n + 1)` with a
/// non-real argument (refused before 0.56: "rising_factorial of complex
/// argument not yet supported"), each part to its own precision (up to
/// `sensitivity::gap_cap` bits below the modulus).  An integer `|n| ≤ 4096`
/// is the product `x(x + 1)⋯(x + n − 1)` (`1/((x − 1)⋯(x − |n|))` for `n <
/// 0`; `ff` with `x − k`), whose relative error is at most `4|n|` units of the
/// working precision; otherwise `e^{ln Γ(u) − ln Γ(d)}` with the principal
/// `loggamma` of a non-real argument and `ln|Γ| + iπ⌊x⌋` of a real one (the
/// branch does not matter in the exponential), at the bits of the
/// logarithms' magnitude (their absolute error is the value's relative
/// one).  A pole of `Γ(d)` (a real `d`) makes the value 0, one of `Γ(u)`
/// infinite (refused).  mpmath: `rf(1 + i, 3) = 10i`, `rf(2, i) =
/// 0.6529654964 + 0.3430658398i`, `ff(1 + i, 1/2) = 1.198230682 +
/// 0.4125896629i`.
pub(super) fn eval_complex(
    f: LibFn,
    args: &[ExprId],
    cache: Cache<'_>,
    w: &mut Wp<'_>,
) -> Result<Complex, SymplexError> {
    let prec = w.bits();
    let (Some(x), Some(n)) = (
        args.first().and_then(|a| cache.get(a)),
        args.get(1).and_then(|a| cache.get(a)),
    ) else {
        return Err(unevaluable("rising/falling factorial lost its arguments"));
    };
    let rising = f == LibFn::RisingFactorial;
    let cap = super::sensitivity::gap_cap(prec);
    let mut extra = 32usize;
    for _ in 0..4 {
        let wp = prec + extra;
        let (value, lost) = complex_quotient(rising, x, n, &mut w.at(wp))?;
        let gap = part_gap(&value).min(cap);
        if lost + gap + 16 <= extra {
            return Ok((
                super::round_to(value.0, prec),
                super::round_to(value.1, prec),
            ));
        }
        if lost > 4 * prec + 512 {
            break;
        }
        extra = lost + gap + 32;
    }
    Err(unevaluable(format!(
        "{}: the value cancels beyond the working precision",
        f.name()
    )))
}

/// Is the value of [`eval_complex`] at `prec` bits exact: exact arguments and
/// an integer `n` whose product (`n = 0`: the empty product 1) rounds
/// nowhere at `prec` bits?  (`rising_factorial(I, 0) = 1` must stay an
/// exactly real argument: `rising_factorial(rising_factorial(I, 0), y)` was
/// refused.)
pub(super) fn exact_product(f: LibFn, args: &[ExprId], cache: Cache<'_>, prec: usize) -> bool {
    let (Some(x), Some(n)) = (
        args.first().and_then(|a| cache.get(a)),
        args.get(1).and_then(|a| cache.get(a)),
    ) else {
        return false;
    };
    if !args
        .iter()
        .all(|a| cache.bound(a).is_some_and(|b| b.is_exact()))
    {
        return false;
    }
    let Some(k) = integer_steps(n) else {
        return false;
    };
    if k == 0 {
        return true;
    }
    if k < 0 || [&x.0, &x.1].iter().any(|p| p.inexact()) {
        return false;
    }
    let rising = f == LibFn::RisingFactorial;
    let mut acc = (BigFloat::from_i32(1, prec), BigFloat::new(prec));
    for j in 0..k.unsigned_abs() {
        let offset = BigFloat::from_u64(j, prec);
        let re = if rising {
            x.0.add(&offset, prec, RM)
        } else {
            x.0.sub(&offset, prec, RM)
        };
        acc = c_mul(&acc, &(re, x.1.clone()), prec, RM);
        if [&acc.0, &acc.1].iter().any(|p| p.inexact()) {
            return false;
        }
    }
    true
}

/// An exact integer `n` with `|n| ≤ 4096` (the product route).
fn integer_steps(n: &Complex) -> Option<i64> {
    (n.1.is_zero() && (n.0.is_zero() || n.0.is_int()))
        .then(|| crate::base::bigfloat::to_f64(&n.0).and_then(crate::base::cast::f64_to_i64))
        .flatten()
        .filter(|k| k.unsigned_abs() <= 4096)
}

/// How many bits below the modulus of `v` its smaller nonzero part lies.
fn part_gap(v: &Complex) -> usize {
    let e = |x: &BigFloat| {
        (!x.is_zero())
            .then(|| x.exponent().map(i64::from))
            .flatten()
    };
    match (e(&v.0), e(&v.1)) {
        (Some(a), Some(b)) => usize::try_from((a - b).abs()).unwrap_or(usize::MAX / 8),
        _ => 0,
    }
}

/// The quotient of [`eval_complex`] at the working precision, and the bits
/// its route loses (relative to the modulus).
fn complex_quotient(
    rising: bool,
    x: &Complex,
    n: &Complex,
    w: &mut Wp<'_>,
) -> Result<(Complex, usize), SymplexError> {
    let wp = w.bits();
    let one = (BigFloat::from_i32(1, wp), BigFloat::new(wp));
    if let Some(k) = integer_steps(n) {
        // x(x + 1)⋯ or x(x − 1)⋯, k factors; 1/((x ∓ 1)⋯) for k < 0.
        let up = rising == (k > 0);
        let mut acc = one.clone();
        for j in 0..k.unsigned_abs() {
            let off = if k > 0 { j } else { j + 1 };
            let offset = BigFloat::from_u64(off, wp);
            let factor = if up {
                (x.0.add(&offset, wp, RM), x.1.clone())
            } else {
                (x.0.sub(&offset, wp, RM), x.1.clone())
            };
            acc = c_mul(&acc, &factor, wp, RM);
        }
        let lost = usize::try_from(super::accuracy::ceil_log2(k.unsigned_abs() as usize + 1))
            .unwrap_or(0)
            + 3;
        if k >= 0 {
            return Ok((acc, lost));
        }
        if acc.0.is_zero() && acc.1.is_zero() {
            return Err(unevaluable("the value is infinite (a pole of Gamma)"));
        }
        return Ok((c_div(&one, &acc, wp, RM), lost));
    }
    // Γ(u)/Γ(d).
    let (u, d) = if rising {
        (c_add(x, n, wp, RM), x.clone())
    } else {
        let x1 = c_add(x, &one, wp, RM);
        (x1.clone(), c_sub(&x1, n, wp, RM))
    };
    let pole = |z: &Complex| z.1.is_zero() && is_npint(&z.0);
    if pole(&d) && !pole(&u) {
        return Ok(((BigFloat::new(wp), BigFloat::new(wp)), 0));
    }
    if pole(&u) {
        return Err(unevaluable("the value is infinite (a pole of Gamma)"));
    }
    // A real argument below 2¹⁰: the logarithm of `Γ` itself (`arb_log_gamma`
    // measures its own cancellation and exhausts the precision at the zeros
    // 1, 2 of `ln Γ`); its absolute error is `Γ`'s relative one.
    let log_gamma = |z: &Complex, w: &mut Wp<'_>| -> Result<Complex, SymplexError> {
        if z.1.is_zero() && z.0.exponent().is_none_or(|e| e > 10) {
            super::arb_log_gamma(&z.0, w)
        } else if z.1.is_zero() {
            let p = w.bits();
            let g = (super::arb_gamma_real(&z.0, w)?, BigFloat::new(p));
            Ok(super::c_ln(&g, w))
        } else {
            super::loggamma::loggamma_complex(z, w)
        }
    };
    // The logarithms' magnitude first: their absolute error is the value's
    // relative one.
    let rough = |z: &Complex, w: &mut Wp<'_>| -> Result<i64, SymplexError> {
        Ok(super::accuracy::mag(&log_gamma(z, &mut w.at(64))?)
            .unwrap_or(0)
            .max(0))
    };
    let big = rough(&u, w)?.max(rough(&d, w)?) + 2;
    let p = wp + usize::try_from(big).unwrap_or(0);
    let lu = log_gamma(&u, &mut w.at(p))?;
    let ld = log_gamma(&d, &mut w.at(p))?;
    let diff = c_sub(&lu, &ld, p, RM);
    let limit = f64::from(astro_float::EXPONENT_MAX) * std::f64::consts::LN_2;
    if crate::base::bigfloat::to_f64(&diff.0).is_none_or(|r| r.abs() > limit) {
        return Err(unevaluable("Gamma quotient beyond the exponent range"));
    }
    Ok((super::c_exp(&diff, &mut w.at(wp)), 4))
}

/// `C(n, k) = Γ(n + 1)/(Γ(k + 1)·Γ(n − k + 1))` through [`gamma_ratio`]:
/// the limit where poles of Γ meet (`C(−7, 2500) = C(2506, 6)`).  For a
/// negative integer `k` it is 0 for every `n` (Graham–Knuth–Patashnik, as
/// `eval` and `compile`): at a negative integer `n` too, where the quotient
/// of Γ has no limit (mpmath's, along `n − k` fixed, gave `C(−2, −4) = 3`
/// up to 0.52).
pub(super) fn binomial(
    n: &BigFloat,
    k: &BigFloat,
    w: &mut Wp<'_>,
) -> Result<BigFloat, SymplexError> {
    let prec = w.bits();
    if k.is_int() && k.is_negative() && !k.is_zero() {
        return Ok(BigFloat::new(prec));
    }
    let n1 = exact_sum(n, &int(1), prec);
    let k1 = exact_sum(k, &int(1), prec);
    let nk1 = exact_sum(&n1, &k.neg(), prec);
    gamma_ratio(&[n1], &[k1, nk1], w)
}

/// Does `C(n, k)`'s Γ formula meet a pole: is `n + 1`, `k + 1` or
/// `n − k + 1` a non-positive integer?
pub(super) fn binomial_meets_pole(n: &BigFloat, k: &BigFloat, prec: usize) -> bool {
    let n1 = exact_sum(n, &int(1), prec);
    let k1 = exact_sum(k, &int(1), prec);
    let nk1 = exact_sum(&n1, &k.neg(), prec);
    is_npint(&n1) || is_npint(&k1) || is_npint(&nk1)
}

/// The Catalan function `C(x) = Γ(2x + 1)/(Γ(x + 1)·Γ(x + 2))` for real
/// `x` (SymPy's `catalan`, which `evalf`s it this way: `catalan(1/2) =
/// 8/(3π) = 0.848826…`).  At the negative integers the quotient is
/// `0/0`-like with unequal rates (`Γ(2x + 1)` meets its pole twice as fast)
/// and continuous: the limits `C(−1) = −1/2`, `C(n) = 0` for `n ≤ −2`;
/// the negative half-integers are poles.  (Before 0.40 `evalf` had no
/// routine for `catalan` at all.)
pub(super) fn catalan(x: &BigFloat, w: &mut Wp<'_>) -> Result<BigFloat, SymplexError> {
    let prec = w.bits();
    if x.is_nan() || x.is_inf() {
        return Err(unevaluable("catalan of a non-finite argument"));
    }
    if is_npint(x) && !x.is_zero() {
        return Ok(if x.cmp(&int(-1)) == Some(0) {
            int(-1).div(&int(2), prec, RM)
        } else {
            BigFloat::new(prec)
        });
    }
    let numerator = exact_sum(&exact_sum(x, x, prec), &int(1), prec);
    if is_npint(&numerator) {
        return Err(unevaluable(
            "catalan has a pole at a negative half-integer (the value is infinite)",
        ));
    }
    let x1 = exact_sum(x, &int(1), prec);
    let x2 = exact_sum(x, &int(2), prec);
    gamma_ratio(&[numerator], &[x1, x2], w)
}

/// `cos(πx)` to full relative precision: `x = k + f` with `k` the nearest
/// integer and `|f| ≤ 1/2` (both exact), `cos(πx) = (−1)ᵏ·sin(π(1/2 − |f|))`,
/// exactly 0 at the half-integers.  (The `cos` of a rounded `πx` has an
/// absolute error of `2^{(bits of x) − wp}`, all of the value next to a
/// half-integer.)
fn cos_pi(x: &BigFloat, w: &mut Wp<'_>) -> BigFloat {
    let wp = w.bits();
    let half = BigFloat::from_f64(0.5, 64);
    let mut k = x.floor();
    let mut f = exact_sum(x, &k.neg(), wp);
    if f.cmp(&half).is_some_and(|c| c > 0) {
        f = exact_sum(&f, &int(-1), wp);
        k = exact_sum(&k, &int(1), wp);
    }
    let g = exact_sum(&half, &f.abs().neg(), wp);
    let s = if g.is_zero() {
        BigFloat::new(wp)
    } else {
        w.cc.pi(wp, RM).clone().mul(&g, wp, RM).sin(wp, RM, w.cc)
    };
    if is_even(&k) { s } else { s.neg() }
}

/// The Fibonacci (`lucas = false`) or Lucas function of a real `x`, as
/// SymPy extends them off the integers: `F(x) = (φˣ − cos(πx)·φ⁻ˣ)/√5`,
/// `L(x) = φˣ + cos(πx)·φ⁻ˣ` (`fibonacci(1/2) = 0.568864…`; before 0.40
/// `evalf` had no routine for either).  The two terms cancel near the
/// zeros (`F(0) = 0`); the loss is measured and the sum redone with as
/// many more bits, up to the cancellation cap.  `x·ln φ` is formed at the
/// bits of `x`'s integer part on top, and `cos(πx)` by [`cos_pi`] (with a
/// rounded `πx`, `fibonacci(−2097153/2)` — the tiny `φˣ/√5`, the cosine
/// being 0 — came out near `±10²¹⁸⁰⁰⁰`).
pub(super) fn fibonacci(
    x: &BigFloat,
    lucas: bool,
    w: &mut Wp<'_>,
) -> Result<BigFloat, SymplexError> {
    let prec = w.bits();
    if x.is_nan() || x.is_inf() {
        return Err(unevaluable("fibonacci of a non-finite argument"));
    }
    let int_bits = x
        .exponent()
        .map_or(0, |e| usize::try_from(e).unwrap_or(0))
        .min(1 << 20);
    let cap = super::cancellation_cap(prec);
    let mut extra = 32 + int_bits;
    loop {
        let wp = prec + extra;
        let sqrt5 = int(5).sqrt(wp, RM);
        let phi = int(1).add(&sqrt5, wp, RM).div(&int(2), wp, RM);
        let t = x.mul(&phi.ln(wp, RM, w.cc), wp, RM);
        let up = t.exp(wp, RM, w.cc);
        let down = t.neg().exp(wp, RM, w.cc);
        if up.is_inf() || up.is_nan() || down.is_inf() || down.is_nan() {
            return Err(unevaluable(
                "fibonacci of this argument is beyond the floating-point exponent range",
            ));
        }
        let cosine = cos_pi(x, &mut w.at(wp));
        let tail = cosine.mul(&down, wp, RM);
        let sum = if lucas {
            up.add(&tail, wp, RM)
        } else {
            up.sub(&tail, wp, RM)
        };
        // A zero tail (the cosine at a half-integer) has no exponent to
        // compare: the sum is `up` exactly.
        let top = [&up, &tail]
            .into_iter()
            .filter(|v| !v.is_zero())
            .filter_map(|v| v.exponent())
            .max()
            .unwrap_or(0);
        let lost = match sum.exponent() {
            Some(e) if !sum.is_zero() => {
                usize::try_from(i64::from(top) - i64::from(e)).unwrap_or(0)
            }
            _ => wp,
        };
        if lost + 16 <= extra {
            let value = if lucas { sum } else { sum.div(&sqrt5, wp, RM) };
            return Ok(super::round_to(value, prec));
        }
        if extra >= cap {
            return Err(super::special_exhausted(prec));
        }
        extra = (lost + 32).max(2 * extra).min(cap);
    }
}

/// `H(z) = ψ(z + 1) + γ` for real `z` (poles at the negative integers).
/// The terms cancel near `z = 0` (`H(z) ≈ ζ(2)·z`): the loss is measured
/// and the sum repeated with that many more bits, up to the cancellation
/// cap of the other special functions.
/// The double factorial of a real non-integer `x`:
/// `x!! = 2^(x/2)·(π/2)^((cos πx − 1)/4)·Γ(x/2 + 1)`, the continuation
/// that takes the values of `n!!` at every integer where those are defined
/// (`(−3)!! = −1`, `(−5)!! = 1/3`, as `eval` folds them) — mpmath's
/// `fac2`, Mathematica's `Factorial2`.  Before 0.54 refused (`factorial2(1/2)`;
/// SymPy raises there).  Integers are `eval`'s; the poles `x = −2, −4, …`
/// are integers too.
pub(super) fn factorial2(x: &BigFloat, w: &mut Wp<'_>) -> Result<BigFloat, SymplexError> {
    let prec = w.bits();
    if x.is_zero() || x.is_int() || x.is_nan() || x.is_inf() {
        return Err(unevaluable(
            "evalf has no numerical routine for 'factorial2' at an integer (eval folds it exactly)",
        ));
    }
    // |x|·π is rounded relative to itself before the cosine: the bits of
    // its integer part come on top.
    let wp = prec
        + 32
        + x.exponent()
            .map_or(0, |e| usize::try_from(e.max(0)).unwrap_or(0));
    let two = BigFloat::from_i32(2, wp);
    let half_x = x.div(&two, wp, RM);
    let pi = w.cc.pi(wp, RM).clone();
    let c = x.mul(&pi, wp, RM).cos(wp, RM, w.cc);
    let e = c
        .sub(&BigFloat::from_i32(1, wp), wp, RM)
        .div(&BigFloat::from_i32(4, wp), wp, RM);
    let ln_half_pi = pi.div(&two, wp, RM).ln(wp, RM, w.cc);
    let ln_two = two.ln(wp, RM, w.cc);
    let log = half_x
        .mul(&ln_two, wp, RM)
        .add(&e.mul(&ln_half_pi, wp, RM), wp, RM);
    let g = super::arb_gamma_real(
        &half_x.add(&BigFloat::from_i32(1, wp), wp, RM),
        &mut w.at(wp),
    )?;
    let r = log.exp(wp, RM, w.cc).mul(&g, wp, RM);
    if r.is_inf() || r.is_nan() {
        return Err(unevaluable(
            "factorial2: the value leaves the exponent range",
        ));
    }
    Ok(super::round_to(r, prec))
}

pub(super) fn harmonic(z: &BigFloat, w: &mut Wp<'_>) -> Result<BigFloat, SymplexError> {
    let prec = w.bits();
    if z.is_nan() || z.is_inf() {
        return Err(unevaluable("harmonic of a non-finite argument"));
    }
    if z.is_zero() {
        return Ok(BigFloat::new(prec));
    }
    if z.is_int() && z.is_negative() {
        return Err(unevaluable(
            "harmonic has a pole at a negative integer (the value is infinite)",
        ));
    }
    let z1 = exact_sum(z, &int(1), prec);
    let cap = super::cancellation_cap(prec);
    let mut extra = 32;
    loop {
        let wp = prec + extra;
        let psi = super::arb_digamma(&z1, &mut w.at(wp))?;
        let gamma = super::arb_euler_gamma(&mut w.at(wp))?;
        let h = psi.add(&gamma, wp, RM);
        let top = psi
            .exponent()
            .unwrap_or(0)
            .max(gamma.exponent().unwrap_or(0));
        let lost = match h.exponent() {
            Some(e) if !h.is_zero() => usize::try_from(i64::from(top) - i64::from(e)).unwrap_or(0),
            _ => wp,
        };
        if lost + 16 <= extra {
            return Ok(super::round_to(h, prec));
        }
        if extra >= cap {
            return Err(super::special_exhausted(prec));
        }
        extra = (lost + 32).max(2 * extra).min(cap);
    }
}

/// The largest order `polygamma_large` accepts (`s = n + 1` stays exact in
/// an `f64`).
const MAX_ORDER: f64 = 9.0e15;

/// `ψ⁽ⁿ⁾(x)` for an integer order `n` (beyond the recurrence and
/// asymptotic series of `arb_polygamma`, which handles `n ≤ 10⁴`) and a real
/// `x` off the poles: `(−1)^{n+1}·exp(ln Γ(n + 1) − s·ln d + ln Z)·sign Z`,
/// `ζ(s, x) = d^{−s}·Z` (see the module documentation).
///
/// # Errors
///
/// [`SymplexError::Unevaluable`] at a pole or beyond the floating-point
/// exponent range; [`SymplexError::PrecisionExhausted`] when the pairs of a
/// negative `x` cancel beyond the cap.
pub(super) fn polygamma_large(
    n: &BigFloat,
    x: &BigFloat,
    w: &mut Wp<'_>,
) -> Result<BigFloat, SymplexError> {
    let prec = w.bits();
    if x.is_nan() || x.is_inf() {
        return Err(unevaluable("polygamma of a non-finite argument"));
    }
    if is_npint(x) {
        return Err(unevaluable("polygamma at non-positive integer pole"));
    }
    let s = exact_sum(n, &int(1), prec);
    let sf = super::f64_of(&s)?;
    if !(2.0..=MAX_ORDER).contains(&sf) {
        return Err(SymplexError::NotImplemented(format!(
            "polygamma of order {} (at most {MAX_ORDER:e})",
            sf - 1.0
        )));
    }
    let odd = !is_even(&s);
    let cap = super::cancellation_cap(prec);
    let mut extra = 32;
    loop {
        let wp = prec + extra;
        let (scale, z) = scaled_hurwitz(&s, sf, odd, x, &mut w.at(wp))?;
        // `Z` is a sum of terms of magnitude at most 1: its loss is `−log₂|Z|`.
        let lost = match z.exponent() {
            Some(e) if !z.is_zero() => usize::try_from(-i64::from(e)).unwrap_or(0),
            _ => wp,
        };
        if lost + 16 > extra {
            if extra >= cap {
                return Err(super::special_exhausted(prec));
            }
            extra = (lost + 32).max(2 * extra).min(cap);
            continue;
        }
        // ln|ψ| = ln Γ(s) − s·ln d + ln|Z|, to an absolute error below 2^{−prec−32}.
        let lg = super::arb_log_gamma(&s, &mut w.at(wp))?.0;
        let ln_d = scale.ln(wp, RM, w.cc);
        let size = super::f64_of(&lg)?.abs() + sf * super::f64_of(&ln_d)?.abs();
        let lp = prec + 32 + (size + 1.0).log2().ceil().max(0.0) as usize;
        let lg = super::arb_log_gamma(&s, &mut w.at(lp))?.0;
        let s_ln_d = s.mul(&scale.ln(lp, RM, w.cc), lp, RM);
        let ln_z = z.abs().ln(lp, RM, w.cc);
        let total = lg.sub(&s_ln_d, lp, RM).add(&ln_z, lp, RM);
        let v = total.exp(lp, RM, w.cc);
        if v.is_inf() || v.is_nan() || v.is_zero() {
            return Err(unevaluable(
                "polygamma of this order is beyond the floating-point exponent range",
            ));
        }
        // (−1)^{n+1} = (−1)^s.
        let negative = odd != z.is_negative();
        let v = super::round_to(v, prec);
        return Ok(if negative { v.neg() } else { v });
    }
}

/// `ζ(s, x) = d^{−s}·Z`: the scale `d` (the smallest `|x + k|` contributing)
/// and `Z` at `wp` bits, whose largest term is 1.
fn scaled_hurwitz(
    s: &BigFloat,
    sf: f64,
    odd: bool,
    x: &BigFloat,
    w: &mut Wp<'_>,
) -> Result<(BigFloat, BigFloat), SymplexError> {
    let wp = w.bits();
    if pos(x) {
        return Ok((x.clone(), positive_part(s, sf, x, x, w)?));
    }
    // x = f − m, f ∈ (0, 1), m ≥ 1: the terms left of 0 are (x + j)^{−s} =
    // (−1)^s·(i − f)^{−s}, i = m − j = 1 … m.
    let m = x.floor().neg();
    let f = exact_sum(x, &m, wp);
    let one_minus_f = exact_sum(&int(1), &f.neg(), wp);
    let mf = super::f64_of(&m)?;
    if odd && f == one_minus_f {
        // A half-integer x: the pairs (f + j, 1 − f + j), j < m, cancel exactly.
        let start = exact_sum(&f, &m, wp);
        return Ok((start.clone(), positive_part(s, sf, &start, &start, w)?));
    }
    let g = if f < one_minus_f {
        f.clone()
    } else {
        one_minus_f.clone()
    };
    let negligible = |t: &BigFloat, sum: &BigFloat| -> bool {
        match (t.exponent(), sum.exponent()) {
            (_, _) if t.is_zero() => true,
            (Some(te), Some(se)) if !sum.is_zero() => {
                i64::from(se) - i64::from(te) > wp as i64 + 32
            }
            _ => false,
        }
    };
    let mut z = BigFloat::new(wp);
    let mut j = 0.0f64;
    let mut u = f.clone(); // f + j
    let mut v = one_minus_f.clone(); // 1 − f + j
    while j < mf {
        let a = scaled_term(s, &u, &g, w);
        let b = scaled_term(s, &v, &g, w);
        let pair = if odd {
            a.sub(&b, wp, RM)
        } else {
            a.add(&b, wp, RM)
        };
        z = z.add(&pair, wp, RM);
        if negligible(&a, &z) && negligible(&b, &z) {
            // The rest of both sums: below `(1 + u/(s − 1))` times these
            // terms each (the bound of `positive_part`), under `2^{−wp−8}`
            // of the sum for `|x| ≤ 10⁶`, `s > 10⁴`.
            return Ok((g, z));
        }
        j += 1.0;
        u = exact_sum(&u, &int(1), wp);
        v = exact_sum(&v, &int(1), wp);
        if j > 1.0e6 {
            return Err(SymplexError::NotImplemented(
                "polygamma of large order far left of 0 (more than a million terms)".into(),
            ));
        }
    }
    // The terms right of the last pair: f + m, f + m + 1, …
    let rest = positive_part(s, sf, &u, &g, w)?;
    Ok((g, z.add(&rest, wp, RM)))
}

/// `(g/u)^s = exp(−s·ln(u/g))` at `wp` bits (`u ≥ g > 0`).
fn scaled_term(s: &BigFloat, u: &BigFloat, g: &BigFloat, w: &mut Wp<'_>) -> BigFloat {
    let wp = w.bits();
    // `s·ln(u/g)` to an absolute `2^{−wp−8}` where the term matters
    // (`s·ln(u/g) ≲ wp`): `log₂ s` bits for the product, `2·log₂ wp` more.
    let s_bits = s.exponent().map_or(0, |e| usize::try_from(e).unwrap_or(0));
    let p = wp + s_bits + 2 * (usize::BITS - wp.leading_zeros()) as usize + 16;
    let l = u.div(g, p, RM).ln(p, RM, w.cc);
    s.mul(&l, p, RM).neg().exp(wp, RM, w.cc)
}

/// `Σ_{j≥0} (g/(f + j))^s` for `f ≥ g > 0` and a large `s`: the terms while
/// they matter, with the rigorous bound `(g/u)^s·(1 + u/(s − 1))` on the sum
/// from `u` on, or the Euler–Maclaurin tail once `f + j ≥ max(s/π, wp)`.
fn positive_part(
    s: &BigFloat,
    sf: f64,
    f: &BigFloat,
    g: &BigFloat,
    w: &mut Wp<'_>,
) -> Result<BigFloat, SymplexError> {
    let wp = w.bits();
    let switch = (sf / std::f64::consts::PI).max(wp as f64) + 1.0;
    let mut sum = BigFloat::new(wp);
    let mut u = f.clone();
    for _ in 0..1_000_000 {
        let uf = super::f64_of(&u)?;
        if uf >= switch {
            let tail = euler_maclaurin_tail(s, &u, g, w);
            return Ok(sum.add(&tail, wp, RM));
        }
        let t = scaled_term(s, &u, g, w);
        sum = sum.add(&t, wp, RM);
        let u1 = exact_sum(&u, &int(1), wp);
        let next = scaled_term(s, &u1, g, w);
        // Σ_{i≥0} (u1 + i)^{−s} ≤ u1^{−s}·(1 + u1/(s − 1)).
        let factor = 1.0 + (uf + 1.0) / (sf - 1.0);
        let bound_exp =
            next.exponent().map_or(i64::MIN / 2, i64::from) + factor.log2().ceil() as i64 + 1;
        let sum_exp = sum.exponent().map_or(i64::MIN / 2, i64::from);
        if next.is_zero() || sum_exp - bound_exp > wp as i64 + 8 {
            return Ok(sum);
        }
        u = u1;
    }
    Err(SymplexError::NotImplemented(
        "polygamma of large order: the Hurwitz zeta sum did not settle in a million terms".into(),
    ))
}

/// `Σ_{i≥0} (g/(a + i))^s` by Euler–Maclaurin at `a ≥ max(s/π, wp)`:
/// `(g/a)^s·[a/(s − 1) + 1/2 + Σ_{j≥1} B₂ⱼ/(2j)!·(s)₂ⱼ₋₁/a^{2j−1}]`, the
/// series stopped once a term is below `2^{−wp−8}` of the bracket (for the
/// completely monotone `t^{−s}` the remainder is below the first omitted
/// term).
fn euler_maclaurin_tail(s: &BigFloat, a: &BigFloat, g: &BigFloat, w: &mut Wp<'_>) -> BigFloat {
    let wp = w.bits();
    let one = BigFloat::from_i32(1, wp);
    let s_minus_1 = s.sub(&one, wp, RM);
    let mut bracket = a
        .div(&s_minus_1, wp, RM)
        .add(&BigFloat::from_f64(0.5, wp), wp, RM);
    let a2 = a.mul(a, wp, RM);
    // (s)_{2j−1}/((2j)!·a^{2j−1}), j = 1: s/(2a).
    let mut coef = s.div(&a.mul(&BigFloat::from_i32(2, wp), wp, RM), wp, RM);
    let mut prev: Option<i64> = None;
    for j in 1usize..=(wp / 2 + 64) {
        if j > 1 {
            // × (s + 2j − 3)(s + 2j − 2)/((2j − 1)(2j)·a²)
            let jj = j as i64;
            let p1 = s.add(&BigFloat::from_i64(2 * jj - 3, 64), wp, RM);
            let p2 = s.add(&BigFloat::from_i64(2 * jj - 2, 64), wp, RM);
            let q = BigFloat::from_i64((2 * jj - 1) * (2 * jj), 64);
            coef = coef
                .mul(&p1, wp, RM)
                .mul(&p2, wp, RM)
                .div(&q, wp, RM)
                .div(&a2, wp, RM);
        }
        let b = super::ratio_to_bigfloat(&super::bernoulli::even(j), wp);
        let term = b.mul(&coef, wp, RM);
        let te = term.exponent().map_or(i64::MIN / 2, i64::from);
        if prev.is_some_and(|p| te > p) {
            break; // the asymptotic series turned (not reached for a ≥ s/π)
        }
        prev = Some(te);
        bracket = bracket.add(&term, wp, RM);
        let be = bracket.exponent().map_or(0, i64::from);
        if term.is_zero() || be - te > wp as i64 + 8 {
            break;
        }
    }
    scaled_term(s, a, g, w).mul(&bracket, wp, RM)
}