symplex 0.55.0

//! The even Bernoulli numbers `B₂ₖ` of the asymptotic series (`ln Γ`, `ψ`,
//! `ψ⁽ⁿ⁾` by Stirling's series), from the tangent numbers.
//!
//! The series need `B₂ₖ` up to `k ≈ 0.1·p` at a working precision of `p`
//! bits, and a value that is zero to the precision is pursued to about
//! 3,200 bits (`ZERO_SEARCH_BITS`): some 350 numbers.  The exact
//! rational recurrence of `base::bernoulli` (`B_m = −Σ C(m+1, k)·B_k/(m+1)`,
//! a sum of rationals per number) took 5–8 s to reach them.  The tangent
//! numbers `T₁ = 1, T₂ = 2, T₃ = 16, …` need only integer updates by small
//! factors, `O(n²)` of them for `T₁ … Tₙ` (Brent and Harvey, "Fast
//! computation of Bernoulli, Tangent and Secant numbers", 2011,
//! Algorithm TangentNumbers — the published algorithm, not code), and
//!
//! ```text
//! B₂ₖ = (−1)^(k−1) · 2k · Tₖ / (2^(2k) · (2^(2k) − 1)).
//! ```
//!
//! The exact evaluator's `bernoulli(n)` (`eval::exact_bernoulli`) takes its
//! even numbers from here too.  Since 0.40 the algorithm lives in
//! `base::bernoulli`, which every layer reaches (`ntheory::bernoulli`,
//! series and summation used the rational recurrence until then): one
//! implementation and one cache.

use super::wp::{RM, Wp};
use astro_float::BigFloat;

use crate::base::errors::SymplexError;
use crate::base::numeric::Q;

/// `B₂ₖ` for `k ≥ 1` (`B₀ = 1` for `k = 0`), from the tangent numbers
/// ([`crate::base::bernoulli::even`]).
pub(crate) fn even(k: usize) -> Q {
    crate::base::bernoulli::even(k)
}

/// The Bernoulli function `B(s) = −s·ζ(1 − s)` of a real `s` (SymPy's
/// `bernoulli(s)`, which `evalf`s it so: `bernoulli(1/2) = 0.730177…`,
/// `bernoulli(−1) = ζ(2)`).  An entire function, it is the Bernoulli
/// number `Bₙ` at every integer `n ≥ 0` in this crate's convention,
/// `B₁ = +1/2` (`Bₙ = Bₙ(1)`, SymPy's since 1.12).  Up to 0.49 the crate's
/// `B₁` was `−1/2`, returned at `s = 1` exactly: `bernoulli(1)` was `−1/2`
/// and `bernoulli(1 + 10⁻³⁰)` `+0.5`.  (Before 0.40 `evalf` had no routine
/// for `bernoulli`.)
pub(super) fn bernoulli_function(s: &BigFloat, w: &mut Wp<'_>) -> Result<BigFloat, SymplexError> {
    let prec = w.bits();
    if s.is_nan() || s.is_inf() {
        return Err(SymplexError::Unevaluable {
            reason: "bernoulli of a non-finite argument".into(),
        });
    }
    if s.is_zero() {
        return Ok(BigFloat::from_i32(1, prec));
    }
    if *s == BigFloat::from_i32(1, 64) {
        return Ok(BigFloat::from_i32(1, prec).div(&BigFloat::from_i32(2, prec), prec, RM));
    }
    // `1 − s` exactly: next to `s = 0` the pole of `ζ` at 1 is cancelled
    // by the factor `s`, both relative to `s`.
    let wp = prec + 32;
    let one_minus_s = BigFloat::from_i32(1, 64).sub(s, super::exact_bits(s, wp), RM);
    let zeta = super::arb_zeta(&one_minus_s, &mut w.at(wp))?;
    let r = s.mul(&zeta, wp, RM).neg();
    Ok(super::round_to(r, prec))
}

/// Largest `|s|` the Euler function is evaluated at (`Γ(s + 1)` and the
/// Hurwitz sums of the order `−s` beyond).
const MAX_EULER_ORDER: f64 = 1.0e6;

/// The Euler function `E(s)` of a real `s` (SymPy's `euler(s)`, the Euler
/// numbers `Eₙ` at the non-negative integers), as SymPy's `euler._eval_evalf`
/// defines it (`sympy/functions/combinatorial/numbers.py`, BSD):
///
/// ```text
/// E(s) = 2^(s+1)·(ζ(−s, ½) − 2^(s+1)·ζ(−s, ¾)),   E(−1) = π/2.
/// ```
///
/// For `s ≥ ½` the Hurwitz functions of the negative order `−s` cancel about
/// `1.8·s` bits; Hurwitz's formula `ζ(−s, a) = 2Γ(s+1)/(2π)^(s+1)·Σₙ
/// cos(π(s+1)/2 − 2πna)/n^(s+1)` (DLMF 25.11.9) turns the combination into
///
/// ```text
/// E(s) = 2·cos(πs/2)·Γ(s+1)/(2π)^(s+1)·(ζ(s+1, ¼) − ζ(s+1, ¾))
/// ```
///
/// (`= 2Γ(s+1)(2/π)^(s+1)·cos(πs/2)·β(s+1)`, Dirichlet's `β`; `E(2) = −1`,
/// `E(4) = 5`), whose Hurwitz sums of positive order do not cancel.  The odd
/// integers are exact zeros (`cos(πs/2)` from an exact reduction), the even
/// ones up to [`crate::base::bernoulli::EXACT_EULER_N`] exact.  The loss of
/// the Hurwitz sums is measured and the evaluation repeated with that many
/// more bits, as for `polylog`.  (Before 0.44 `evalf` had no routine for
/// `euler_number`: `euler_number(1/2)` was refused.)
pub(super) fn euler_function(s: &BigFloat, w: &mut Wp<'_>) -> Result<BigFloat, SymplexError> {
    let prec = w.bits();
    let unevaluable = |reason: &str| SymplexError::Unevaluable {
        reason: reason.into(),
    };
    if s.is_nan() || s.is_inf() {
        return Err(unevaluable("euler_number of a non-finite argument"));
    }
    let s_f = super::f64_of(s)?;
    if s_f.abs() > MAX_EULER_ORDER {
        return Err(unevaluable(
            "euler_number of an argument beyond 10^6 is not evaluated numerically",
        ));
    }
    if s.is_zero() {
        return Ok(BigFloat::from_i32(1, prec));
    }
    if s.is_int() && !s.is_negative() {
        let n = s_f as usize;
        if n % 2 == 1 {
            return Ok(BigFloat::new(prec));
        }
        if n <= crate::base::bernoulli::EXACT_EULER_N {
            let e = crate::base::combinatorics::euler_even(n / 2);
            return Ok(super::round_to(
                crate::base::numeric::bigint_to_bigfloat(&e, prec + 8),
                prec,
            ));
        }
    }
    if *s == BigFloat::from_i32(-1, 64) {
        let pi = w.cc.pi(prec + 8, RM).clone();
        return Ok(super::round_to(
            pi.div(&BigFloat::from_i32(2, 64), prec + 8, RM),
            prec,
        ));
    }
    let cap = super::cancellation_cap(prec);
    let size = s_f.abs() * (s_f.abs().ln().max(0.0) + 2.0) + 2.0;
    let mut extra = (32 + super::magnitude_bits(size)).min(cap);
    loop {
        let wp = prec + extra;
        let (v, largest) = if s_f >= 0.5 {
            euler_positive(s, &mut w.at(wp))?
        } else {
            euler_negative(s, &mut w.at(wp))?
        };
        // |v| ≥ 2^(e − 1).
        let lv = lg(&v) - 1.0;
        let lost = if lv.is_finite() && largest.is_finite() {
            ((largest - lv).max(0.0).ceil() as usize).min(cap + 1)
        } else {
            wp
        };
        if lost + 24 <= extra {
            return Ok(super::round_to(v, prec));
        }
        if extra >= cap || lost + 24 > cap {
            return Err(super::special_exhausted(prec));
        }
        extra = (lost + 40).max(2 * extra).min(cap);
    }
}

/// An upper bound on `log₂|x|` (its binary exponent), `−∞` at 0.
fn lg(x: &BigFloat) -> f64 {
    if x.is_zero() {
        return f64::NEG_INFINITY;
    }
    x.exponent().map_or(f64::NEG_INFINITY, f64::from)
}

/// `E(s)` for `s ≥ ½` by Hurwitz's formula (see [`euler_function`]) at `wp`
/// bits, and `log₂` of the largest term of the Hurwitz sums scaled as the
/// value.
fn euler_positive(s: &BigFloat, w: &mut Wp<'_>) -> Result<(BigFloat, f64), SymplexError> {
    let wp = w.bits();
    let one = BigFloat::from_i32(1, 64);
    let x = s.add(&one, super::exact_bits(s, wp), RM);
    let cos = super::sin_half_pi(&x, w);
    if cos.is_zero() {
        return Ok((BigFloat::new(wp), f64::NEG_INFINITY));
    }
    let quarter = BigFloat::from_f64(0.25, 64);
    let three_quarters = BigFloat::from_f64(0.75, 64);
    let (h1, l1) = super::polylog::hurwitz_real(&x, &quarter, w)?;
    let (h2, l2) = super::polylog::hurwitz_real(&x, &three_quarters, w)?;
    // Γ(s+1)/(2π)^(s+1) = exp(ln Γ(x) − x·ln 2π), the exponent to an absolute
    // error of 2^−wp (the guard bits of `euler_function` cover its size).
    let lw = wp + 16;
    let two_pi =
        w.cc.pi(lw, RM)
            .clone()
            .mul(&BigFloat::from_i32(2, 64), lw, RM);
    let ln_gamma = super::stirling_log_gamma(&x, &mut w.at(lw))?;
    let scale = ln_gamma
        .sub(&x.mul(&two_pi.ln(lw, RM, w.cc), lw, RM), lw, RM)
        .exp(wp, RM, w.cc)
        .mul(&cos, wp, RM)
        .mul(&BigFloat::from_i32(2, 64), wp, RM);
    let diff = h1.sub(&h2, wp, RM);
    Ok((scale.mul(&diff, wp, RM), lg(&scale) + l1.max(l2)))
}

/// `E(s)` for `s < ½`, `s ≠ −1`, by SymPy's formula (see
/// [`euler_function`]) at `wp` bits, and `log₂` of the largest term of the
/// Hurwitz sums scaled as the value.
fn euler_negative(s: &BigFloat, w: &mut Wp<'_>) -> Result<(BigFloat, f64), SymplexError> {
    let wp = w.bits();
    let exact = super::exact_bits(s, wp);
    let sigma = s.neg();
    let half = BigFloat::from_f64(0.5, 64);
    let three_quarters = BigFloat::from_f64(0.75, 64);
    let (h1, l1) = super::polylog::hurwitz_real(&sigma, &half, w)?;
    let (h2, l2) = super::polylog::hurwitz_real(&sigma, &three_quarters, w)?;
    let s1 = s.add(&BigFloat::from_i32(1, 64), exact, RM);
    let p = super::bf_pow(&BigFloat::from_i32(2, 64), &s1, &mut w.at(wp + 16));
    let lp = lg(&p);
    let inner = h1.sub(&p.mul(&h2, wp, RM), wp, RM);
    Ok((p.mul(&inner, wp, RM), lp + l1.max(lp + l2)))
}

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

    /// The first even Bernoulli numbers.
    #[test]
    fn first_values() {
        let q = |n: i64, d: i64| Ratio::new(BigInt::from(n), BigInt::from(d));
        assert_eq!(even(1), q(1, 6));
        assert_eq!(even(2), q(-1, 30));
        assert_eq!(even(3), q(1, 42));
        assert_eq!(even(6), q(691, -2730));
    }
}