dashu_float/math/

cache.rs

//! Opt-in cache of mathematical constants, enabling progressive refinement.

use core::fmt;

use dashu_base::{BitTest, EstimatedLog2, Sign, UnsignedAbs};
use dashu_int::{IBig, UBig};

use crate::error::assert_limited_precision;
use crate::fbig::FBig;
use crate::repr::{Context, Repr, Word};
use crate::round::{Round, Rounded};
use crate::utils::ceil_usize;

/// Binary-splitting tree state — exact integers, losslessly extensible.
///
/// Represents `binary_split(start, num_terms)` as the universal triple
/// `(P, Q, T)`, where `start` is 0 for π and 1 for `L(n)` (whose `k=0` term
/// `1/n` is pulled out). These are pure integers: independent of base and
/// rounding mode. To extend to `new_terms > num_terms`, compute the right half
/// over the new range and merge with the universal `T' = T_l·Q_r + P_l·T_r`.
#[derive(Clone)]
pub(crate) struct CachedState {
    pub p: UBig,
    pub q: UBig,
    pub t: IBig,
    pub num_terms: usize,
}

/// An opt-in cache of mathematical constants.
///
/// Holds exact binary-splitting tree state so that repeated calls at increasing
/// precision *extend* prior work instead of recomputing from scratch. For
/// example, computing π at 100 digits and then at 1000 digits only pays for the
/// extra ~900 digits of work.
///
/// The cache is **base-free**: a single [`ConstCache`] serves any base. The base
/// and rounding mode are specified on each method call (e.g.
/// `cache.pi::<10, HalfAway>(100)` for 100 decimal digits).
///
/// `ConstCache` is a plain struct of big integers, so it is `Send + Sync`. The
/// methods take `&mut self` (they extend the cached state on a miss), so a caller
/// either owns one directly, or — to share it across many values and operations —
/// wraps it in `Rc<RefCell<ConstCache>>` as the
/// [`CachedFBig`](crate::CachedFBig) type does. To share one cache across
/// threads, wrap a `ConstCache` (or a `CachedFBig`) in `Arc<Mutex<..>>`.
///
/// # Examples
///
/// ```
/// use dashu_float::ConstCache;
/// use dashu_float::round::mode::HalfAway;
///
/// let mut cache = ConstCache::new();
/// // first call computes from scratch
/// let _pi_100 = cache.pi::<10, HalfAway>(100).value();
/// // second call at higher precision extends the cached state
/// let pi_1000 = cache.pi::<10, HalfAway>(1000).value();
/// assert!(pi_1000.to_string().starts_with("3.141592653589793"));
/// ```
pub struct ConstCache {
    pi: Option<CachedState>,
    /// `L(6)`, `L(9)`, `L(99)` — the sub-series used by ln2 / ln10.
    iacoth_6: Option<CachedState>,
    iacoth_9: Option<CachedState>,
    iacoth_99: Option<CachedState>,
    /// Base-free integer `floor(sqrt(10005) · 2^sqrt_10005_bits)`, reused by π.
    /// Unlike the series slots this holds a plain value (not a `(P,Q,T)` triple) and
    /// is extended by a fresh Karatsuba `UBig::sqrt` — Newton refinement would be no
    /// faster, since `UBig::sqrt` is already O(M(n)).
    sqrt_10005: Option<UBig>,
    sqrt_10005_bits: usize,
}

impl ConstCache {
    /// Create an empty cache.
    pub const fn new() -> Self {
        Self {
            pi: None,
            iacoth_6: None,
            iacoth_9: None,
            iacoth_99: None,
            sqrt_10005: None,
            sqrt_10005_bits: 0,
        }
    }

    /// `floor(sqrt(10005) · 2^bits)` as a base-free integer, cached and extended on
    /// demand. Used by [`pi`](Self::pi). Computed via Karatsuba `UBig::sqrt` (O(M(n))).
    /// Returns the value together with the number of bits it actually corresponds to
    /// (which may be larger than requested, when a higher-precision value is reused).
    fn sqrt_10005(&mut self, bits: usize) -> (UBig, usize) {
        if bits > self.sqrt_10005_bits {
            let n = UBig::from(10005u32) << (2 * bits);
            self.sqrt_10005 = Some(dashu_base::SquareRoot::sqrt(&n));
            self.sqrt_10005_bits = bits;
        }
        (self.sqrt_10005.as_ref().unwrap().clone(), self.sqrt_10005_bits)
    }

    /// π at `precision` base-`B` digits, rounded per `R`. Extends any prior π
    /// state cached in `self`.
    ///
    /// # Panics
    ///
    /// Panics if `precision` is 0.
    #[must_use]
    pub fn pi<const B: Word, R: Round>(&mut self, precision: usize) -> Rounded<FBig<R, B>> {
        assert_limited_precision(precision);

        let bits = bits_for_precision::<B>(precision);
        let num_terms = (bits * 100 / 4708) + 2;

        let (_p, q, t) = extend_or_compute(&mut self.pi, 0, num_terms, chudnovsky_bs);

        // Finalize: π = 426880·√10005·Q / T  (identical to Context::pi)
        let guard_bits = num_terms.bit_len() + 32;
        let work_bits = bits + guard_bits;
        let work_precision = precision_for_bits::<B>(work_bits);
        let work = Context::<R>::new(work_precision);

        // Finalize: π = 426880·√10005·Q / T. With √10005 ≈ isqrt_val·2^(-isqrt_bits)
        // from the base-free cached isqrt, this folds into a single integer ratio
        //   π = (426880 · isqrt_val · Q) / (T · 2^isqrt_bits),
        // avoiding any cross-base conversion of √10005 (convert_int is the fast path,
        // the same one used for Q and T).
        let (isqrt_val, isqrt_bits) = self.sqrt_10005(work_bits);
        let num = IBig::from(426_880) * IBig::from(isqrt_val) * IBig::from(q);
        let den = t << isqrt_bits;
        let num_f = work.convert_int::<B>(num).value();
        let den_f = work.convert_int::<B>(den).value();
        let pi = num_f / den_f;
        pi.with_precision(precision)
    }

    /// `L(n) = acoth(n)` at `precision` base-`B` digits, extending its cached
    /// series state. Only `n ∈ {6, 9, 99}` are cached (the sub-series of ln2 / ln10).
    fn iacoth<const N: u32, const B: Word, R: Round>(&mut self, precision: usize) -> FBig<R, B> {
        // terms until r_k < B^{-p}: (2k+1)·log_B(n) > p. The count is generously
        // over-provisioned (extra terms only add precision), so a plain (truncating)
        // cast suffices in place of a ceiling.
        let log_b_n = N.log2_est() / B.log2_est();
        let required_terms = (precision as f32 / (2.0 * log_b_n)) as usize + 10;

        let slot = match N {
            6 => &mut self.iacoth_6,
            9 => &mut self.iacoth_9,
            99 => &mut self.iacoth_99,
            _ => unreachable!("iacoth only caches n ∈ {{6, 9, 99}}"),
        };
        let (_p, q, t) = extend_or_compute(slot, 1, required_terms, |a, b| iacoth_bs(N, a, b));

        // L(n) = (Q + T) / (n·Q)
        let guard = ceil_usize(precision.log2_est() / B.log2_est()) + 2;
        let work = Context::<R>::new(precision + guard);
        let num = work.convert_int::<B>(q.as_ibig() + &t).value();
        let denom = work.convert_int::<B>(IBig::from(N) * &q).value();
        num / denom
    }

    /// ln(2) at `precision` base-`B` digits, reusing the cached `L(6)` and
    /// `L(99)` sub-series.
    ///
    /// # Panics
    ///
    /// Panics if `precision` is 0.
    #[must_use]
    pub fn ln2<const B: Word, R: Round>(&mut self, precision: usize) -> FBig<R, B> {
        // log(2) = 4·L(6) + 2·L(99)  (Gourdon & Sebah, "Log 2")
        let work = precision + combine_guard::<B>(precision);
        let l6 = self.iacoth::<6, B, R>(work);
        let l99 = self.iacoth::<99, B, R>(work);
        (4u8 * l6 + 2u8 * l99).with_precision(precision).value()
    }

    /// ln(10) at `precision` base-`B` digits, reusing the cached `L(6)`, `L(99)`,
    /// and `L(9)` sub-series.
    ///
    /// # Panics
    ///
    /// Panics if `precision` is 0.
    #[must_use]
    pub fn ln10<const B: Word, R: Round>(&mut self, precision: usize) -> FBig<R, B> {
        // log(10) = 3·log(2) + 2·L(9) = 3·(4·L(6) + 2·L(99)) + 2·L(9)
        //          = 12·L(6) + 6·L(99) + 2·L(9)
        // Flattening avoids the intermediate rounding of ln2 inside the product.
        let work = precision + combine_guard::<B>(precision);
        let l6 = self.iacoth::<6, B, R>(work);
        let l99 = self.iacoth::<99, B, R>(work);
        let l9 = self.iacoth::<9, B, R>(work);
        (12u8 * l6 + 6u8 * l99 + 2u8 * l9)
            .with_precision(precision)
            .value()
    }

    /// ln(B) at `precision` base-`B` digits, reusing the cached ln2 / ln10 where
    /// possible.
    ///
    /// The generic-base arm evaluates `ln(BASE)` through `ln_compute`, whose ball arithmetic can
    /// report a range error. It cannot actually do so for a base that fits a [`Word`]
    /// (`ln(BASE) ≤ Word::BITS · ln 2` and every intermediate stays within a few hundred bits of
    /// it, nowhere near an `isize` exponent limit), so the arm asserts instead of plumbing a
    /// `Result` through this near-correct constant accessor.
    ///
    /// # Panics
    ///
    /// Panics if `precision` is 0, or — for a base that is not 2, 10 or a power of two — if
    /// `ln(BASE)` range-errors (unreachable for any base that fits a [`Word`], see above).
    #[must_use]
    pub fn ln_base<const B: Word, R: Round>(&mut self, precision: usize) -> FBig<R, B> {
        match B {
            2 => self.ln2::<B, R>(precision),
            10 => self.ln10::<B, R>(precision),
            b if b.is_power_of_two() => {
                // ln(2^k) = k·ln(2); evaluate ln2 at elevated precision so the
                // k·ln2 product survives the final round.
                let work = precision + combine_guard::<B>(precision);
                let bits = b.trailing_zeros() as usize;
                (bits * self.ln2::<B, R>(work))
                    .with_precision(precision)
                    .value()
            }
            _ => {
                // generic base: no cached L(n) sub-series applies, so compute ln(B) directly via
                // the near-correct `ln_compute` (the cache only stores a near-correct constant — a
                // cached value is re-rounded on read, and Ziv isn't needed here).
                let ctx = Context::<R>::new(precision);
                ctx.ln_compute::<B>(
                    &Repr::new(Repr::<B>::BASE.into(), 0),
                    // no cache for the generic base (its L(n) isn't cached)
                    precision,
                    false,
                    None,
                )
                .expect("ln(BASE) cannot range-error (the base fits the exponent range)")
                .to_value_radius::<R>(&ctx)
                .0
            }
        }
    }

    /// Sum of `num_terms` across all populated cache slots.
    #[inline]
    pub fn total_terms(&self) -> usize {
        let sum = |s: &Option<CachedState>| s.as_ref().map_or(0, |s| s.num_terms);
        sum(&self.pi) + sum(&self.iacoth_6) + sum(&self.iacoth_9) + sum(&self.iacoth_99)
    }

    /// Sum of word counts across all cached big integers (P, Q, T, and the cached
    /// `√10005` isqrt).
    ///
    /// This reflects the underlying storage words used by the cached state.
    #[inline]
    pub fn total_words(&self) -> usize {
        let slot_words = |s: &Option<CachedState>| {
            s.as_ref().map_or(0, |s| {
                s.p.as_words().len() + s.q.as_words().len() + s.t.as_sign_words().1.len()
            })
        };
        slot_words(&self.pi)
            + slot_words(&self.iacoth_6)
            + slot_words(&self.iacoth_9)
            + slot_words(&self.iacoth_99)
            + self.sqrt_10005.as_ref().map_or(0, |s| s.as_words().len())
    }

    /// Clear all cached constant state, freeing the underlying memory.
    ///
    /// After calling `clear()`, the next constant computation will start from scratch
    /// rather than extending the prior cached state.
    #[inline]
    pub fn clear(&mut self) {
        self.pi = None;
        self.iacoth_6 = None;
        self.iacoth_9 = None;
        self.iacoth_99 = None;
        self.sqrt_10005 = None;
        self.sqrt_10005_bits = 0;
    }
}

impl Default for ConstCache {
    #[inline]
    fn default() -> Self {
        Self::new()
    }
}

/// Ensure `slot` holds state for at least `target` terms, then return `(P, Q, T)`
/// covering `target` terms (or more, when an existing higher-precision state
/// already covers `target` — finalize then rounds down to the requested precision).
///
/// `range_bs(a, b)` computes the leaf-merged state over `[a, b)` and must handle
/// `a == b` by returning the identity `(1, 1, 0)`.
fn extend_or_compute<F>(
    slot: &mut Option<CachedState>,
    start: usize,
    target: usize,
    range_bs: F,
) -> (UBig, UBig, IBig)
where
    F: Fn(usize, usize) -> (UBig, UBig, IBig),
{
    match slot {
        // Already have >= target terms: reuse (extra terms only add precision).
        Some(s) if s.num_terms >= target => (s.p.clone(), s.q.clone(), s.t.clone()),
        // Have fewer terms: extend the right half [num_terms, target) and merge.
        Some(s) => {
            let (pr, qr, tr) = range_bs(s.num_terms, target);
            let (p, q, t) = merge(&s.p, &s.q, &s.t, &pr, &qr, &tr);
            *slot = Some(CachedState {
                p: p.clone(),
                q: q.clone(),
                t: t.clone(),
                num_terms: target,
            });
            (p, q, t)
        }
        // Cold: compute from `start`.
        None => {
            let (p, q, t) = range_bs(start, target);
            *slot = Some(CachedState {
                p: p.clone(),
                q: q.clone(),
                t: t.clone(),
                num_terms: target,
            });
            (p, q, t)
        }
    }
}

/// Reborrow an `Option<&mut ConstCache>` so it can be threaded into several
/// sequential sub-calls. `as_deref_mut` is the natural reborrow here; clippy's
/// `needless_option_as_deref` flags it (the deref target equals the referent),
/// so the lint is allowed at this single centralized spot.
#[inline]
#[allow(clippy::needless_option_as_deref)]
pub(crate) fn reborrow_cache<'a>(
    cache: &'a mut Option<&mut ConstCache>,
) -> Option<&'a mut ConstCache> {
    cache.as_deref_mut()
}

/// Number of bits needed to represent `precision` base-`B` digits exactly.
///
/// For power-of-two bases this is exact; for arbitrary bases it uses the upper
/// bound from [`EstimatedLog2`], which is far tighter than `ilog2(B) + 1`.
fn bits_for_precision<const B: Word>(precision: usize) -> usize {
    if B.is_power_of_two() {
        precision.saturating_mul(B.ilog2() as usize)
    } else {
        // ub ≥ log2(B) with error ≤ 2/256.  Multiply in f64 so the product
        // is exact for precision up to 2^53.  +1 guards float rounding.
        let ub = B.log2_bounds().1;
        ceil_usize(precision as f32 * ub) + 1
    }
}

/// Convert a work-precision expressed in bits back to base-`B` digits.
///
/// For base 2 the identity holds; for power-of-two bases it uses ceiling
/// division; for arbitrary bases it inverts the lower bound from
/// [`EstimatedLog2`] to get a tight ceiling.
fn precision_for_bits<const B: Word>(bits: usize) -> usize {
    if B.is_power_of_two() {
        let log2 = B.ilog2() as usize;
        (bits + log2 - 1) / log2
    } else {
        // lb ≤ log2(B), so 1/lb ≥ 1/log2(B).  +1 guards float rounding.
        let lb = B.log2_bounds().0;
        ceil_usize(bits as f32 / lb) + 1
    }
}

/// Guard digits added when combining sub-series, large enough that the linear
/// combination and its final round to `precision` are unaffected by summation
/// rounding (a few digits cover the constant multipliers and term count).
fn combine_guard<const B: Word>(precision: usize) -> usize {
    ceil_usize(precision.log2_est() / B.log2_est()) + 4
}

/// Universal binary-splitting merge:
/// `combine((P_l,Q_l,T_l), (P_r,Q_r,T_r)) = (P_l·P_r, Q_l·Q_r, T_l·Q_r + P_l·T_r)`.
///
/// This operation is associative, so the `(P, Q, T)` for a range is independent of
/// how the recursion splits it — which is exactly what lets a cached partial tree
/// state be extended by merging in a freshly computed right half.
pub(crate) fn merge(
    pl: &UBig,
    ql: &UBig,
    tl: &IBig,
    pr: &UBig,
    qr: &UBig,
    tr: &IBig,
) -> (UBig, UBig, IBig) {
    let p = pl * pr;
    let q = ql * qr;
    // re-interpret the magnitudes as signed without cloning the big integers
    let t = qr.as_ibig() * tl + pl.as_ibig() * tr;
    (p, q, t)
}

/// Binary splitting implementation for the Chudnovsky series.
/// Returns `(P, Q, T)` for the range `[a, b)`. An empty range `[a, a)` yields the
/// identity `(1, 1, 0)`, so callers may safely merge a cached left state with a
/// right half that starts exactly where the left one ended.
pub(crate) fn chudnovsky_bs(a: usize, b: usize) -> (UBig, UBig, IBig) {
    if a >= b {
        return (UBig::ONE, UBig::ONE, IBig::ZERO);
    }
    if b - a == 1 {
        const COEFF1: IBig = IBig::from_parts_const(Sign::Positive, 13591409);
        const COEFF2: IBig = IBig::from_parts_const(Sign::Positive, 545140134);

        // Base case: calculate single term
        if a == 0 {
            return (UBig::ONE, UBig::ONE, COEFF1);
        }

        let k = a as u64;
        let p = UBig::from(6 * k - 5) * (2 * k - 1) * (6 * k - 1);
        let q = UBig::from(k).pow(3) * 10_939_058_860_032_000u64;
        let t_val = COEFF1 + COEFF2 * k;
        let t_abs = &p * t_val.unsigned_abs();
        let t = IBig::from(t_abs) * Sign::from(a % 2 == 1);
        return (p, q, t);
    }

    // Recursive step
    let mid = (a + b) / 2;
    let (p_l, q_l, t_l) = chudnovsky_bs(a, mid);
    let (p_r, q_r, t_r) = chudnovsky_bs(mid, b);

    // T = T_L * Q_R + T_R * P_L  (the universal merge)
    merge(&p_l, &q_l, &t_l, &p_r, &q_r, &t_r)
}

/// Binary splitting for `L(n) = acoth(n) = Σ_{k≥0} 1/(n^{2k+1}(2k+1))` over `[1, b)`.
///
/// Term ratio (k≥1): `r_k/r_{k-1} = p_k/q_k` with `p_k = 2k-1`, `q_k = (2k+1)·n²`.
/// The `k=0` term `r_0 = 1/n` is pulled out; over `[1, b)` the tree state satisfies
/// `T/Q = n · Σ_{k=1}^{b-1} 1/((2k+1)·n^{2k+1})`, hence `L(n) = (Q + T)/(n·Q)`.
///
/// Using the ratio form (rather than independent `1/q_k` terms) keeps
/// `Q = Π(2k+1)·n²` at O(p) digits: each leaf multiplies only small integers
/// (with `n²` folded in), no `n.pow(2k+1)` per leaf.
pub(crate) fn iacoth_bs(n: u32, a: usize, b: usize) -> (UBig, UBig, IBig) {
    debug_assert!(a >= 1, "iacoth_bs leaf index must be >= 1");
    if a >= b {
        return (UBig::ONE, UBig::ONE, IBig::ZERO); // identity
    }
    // Precomputed initial block [1, 1+K): skip its K leaves on every fresh
    // computation. Because the merge is associative, the constant triple is
    // identical to the recursively computed state regardless of split order.
    // It only applies at the series start (a == 1); recursive/extend calls have
    // a >= 1 + K and never reach this branch.
    if a == 1 {
        if let Some((k, p0, q0, t0)) = iacoth_initial_block(n) {
            // the precomputed block covers [1, 1+k); use it when [a, b) reaches
            // past its end (b > k)
            if b > k {
                let (pr, qr, tr) = iacoth_bs(n, 1 + k, b);
                return merge(&p0, &q0, &t0, &pr, &qr, &tr);
            }
        }
    }
    if b - a == 1 {
        // leaf at k = a (a >= 1): (p_a, q_a, p_a), p_a = 2a-1, q_a = (2a+1)·n²
        let pa = UBig::from(2 * a - 1);
        let n2 = UBig::from(n).pow(2);
        let qa = UBig::from(2 * a + 1) * n2;
        let ta = IBig::from(pa.clone());
        return (pa, qa, ta);
    }
    let mid = (a + b) / 2;
    let (pl, ql, tl) = iacoth_bs(n, a, mid);
    let (pr, qr, tr) = iacoth_bs(n, mid, b);
    merge(&pl, &ql, &tl, &pr, &qr, &tr) // universal merge
}

/// Precomputed binary-splitting state for `L(n) = acoth(n)` over the first `K`
/// terms `[1, 1+K)`, stored as `(K, P, Q, T)`. `K` is chosen (per `n`) so that
/// `P`, `Q` and `|T|` each fit in a `u32`. Since `DoubleWord` is `u32`/`u64`/`u128`
/// for `Word = u16`/`u32`/`u64`, a `u32`-sized literal is accepted by
/// [`UBig::from_dword`] / [`IBig::from_parts_const`] on **every** configuration —
/// so this single set of constants is portable without needing to detect the
/// `Word` width (which is internal to `dashu-int`). The constants also use the
/// inline small-integer representation, so instantiating them never allocates.
///
/// Only the sub-series that back ln2 / ln10 (`n ∈ {6, 9, 99}`) are precomputed;
/// π cannot use this trick because its 2-term `T` already overflows a `u32`.
fn iacoth_initial_block(n: u32) -> Option<(usize, UBig, UBig, IBig)> {
    match n {
        // L(6) over [1, 5): 4 leaves.
        6 => Some(IACOTH_6_INITIAL),
        // L(9) over [1, 4): 3 leaves.
        9 => Some(IACOTH_9_INITIAL),
        // L(99) over [1, 3): 2 leaves.
        99 => Some(IACOTH_99_INITIAL),
        _ => None,
    }
}

/// `(K, P, Q, T)` for `L(6)` over `[1, 5)` (4 leaves).
const IACOTH_6_INITIAL: (usize, UBig, UBig, IBig) = (
    4,
    UBig::from_word(105),
    UBig::from_dword(1587237120),
    IBig::from_parts_const(Sign::Positive, 14946549),
);
/// `(K, P, Q, T)` for `L(9)` over `[1, 4)` (3 leaves).
const IACOTH_9_INITIAL: (usize, UBig, UBig, IBig) = (
    3,
    UBig::from_word(15),
    UBig::from_dword(55801305),
    IBig::from_parts_const(Sign::Positive, 231351),
);
/// `(K, P, Q, T)` for `L(99)` over `[1, 3)` (2 leaves).
const IACOTH_99_INITIAL: (usize, UBig, UBig, IBig) = (
    2,
    UBig::from_word(3),
    UBig::from_dword(1440894015),
    IBig::from_parts_const(Sign::Positive, 49008),
);

/// Binary splitting for `e = Σ_{k≥0} 1/k!` over the index range `[a, b)` with `a ≥ 1`.
///
/// The `k = 0` term `1/0! = 1` is added by the caller. For `k ≥ 1` the term ratio
/// is `tₖ/tₖ₋₁ = 1/k`, so each leaf is `(pₖ, qₖ, pₖ) = (1, k, 1)` and the universal
/// [`merge`] applies verbatim (with `P ≡ 1`, it reduces to `T' = Tₗ·Qᵣ + Tᵣ`).
/// Over `[1, N+1)` the triple satisfies `T/Q = Σ_{k=1}^{N} 1/k!`, hence the
/// finalization `e = (Q + T)/Q` in [`compute_e`].
///
/// This is the optimal algorithm for *e*: there is no Chudnovsky-analog, and the
/// factorial series with binary splitting is `O(M(n) log n)` — faster than π —
/// beating both the AGM approach and `exp(1)` (which pays for argument reduction
/// and `√p` powering that this series needs neither of).
///
/// An empty range yields the identity `(1, 1, 0)`.
pub(crate) fn e_bs(a: usize, b: usize) -> (UBig, UBig, IBig) {
    if a >= b {
        return (UBig::ONE, UBig::ONE, IBig::ZERO);
    }
    if b - a == 1 {
        // leaf at k = a (a ≥ 1): (pₖ, qₖ, pₖ) = (1, a, 1)
        return (UBig::ONE, UBig::from(a), IBig::ONE);
    }
    let mid = (a + b) / 2;
    let (pl, ql, tl) = e_bs(a, mid);
    let (pr, qr, tr) = e_bs(mid, b);
    merge(&pl, &ql, &tl, &pr, &qr, &tr) // universal merge
}

/// Number of series terms `N` for `e = Σ 1/k!` so that the truncation tail is
/// safely below the target precision.
///
/// The tail after `k = N` is `Σ_{k>N} 1/k! < 2/(N+1)!`, so we need
/// `log₂((N+1)!) > bits`. Stirling's lower bound `m! ≥ (m/e)ᵐ` gives
/// `log₂(m!) ≥ m·(log₂ m − log₂ e)`; we seed `m` from the closed-form inverse of
/// that bound and grow until it holds, so the result is always sufficient.
/// Over-provisioning is harmless: extra terms only add precision, which the final
/// `with_precision` discards.
fn e_term_count(bits: usize) -> usize {
    // `log2_bounds` (not `f64::log2`, which is std-only) keeps this `no_std`-clean;
    // using its lower bound makes `lb` a valid lower bound on the Stirling LB, so the
    // loop never terminates too early.
    let target = (bits + 64) as f64; // margin for the tail's factor of 2 and final rounding
    let log2_e = core::f64::consts::LOG2_E;
    let lb = |m: usize| -> f64 { m as f64 * (m.log2_bounds().0 as f64 - log2_e) };
    let denom = (bits.log2_bounds().1 as f64 - log2_e).max(1.0);
    let mut m = (target / denom) as usize + 1;
    while lb(m) < target {
        m += 1;
    }
    // `m` corresponds to `N+1`; return `N`.
    m.saturating_sub(1)
}

/// One-shot computation of *e* (Euler's number) at `precision` base-`B` digits,
/// rounded per `R`.
///
/// Evaluates `e = Σ_{k=0}^∞ 1/k!` by exact-integer binary splitting (the
/// [`e_bs`] kernel) and finalizes `e = (Q + T)/Q`. Unlike π, this is fully
/// self-contained — it depends on no cached sub-constant and is reused by no
/// other operation — so it deliberately does **not** live in [`ConstCache`] and is
/// not progressively refined across calls. (See `Context::e` for the rationale.)
pub(crate) fn compute_e<const B: Word, R: Round>(precision: usize) -> Rounded<FBig<R, B>> {
    assert_limited_precision(precision);

    let bits = bits_for_precision::<B>(precision);
    let n = e_term_count(bits);
    // Σ_{k=1}^{n} 1/k! = T/Q; the k=0 term 1/0! = 1 is folded in as (Q + T)/Q.
    let (_p, q, t) = e_bs(1, n + 1);

    // Evaluate (Q + T)/Q at an elevated working precision so the integer→float
    // conversion and the division round correctly down to `precision`. The guard
    // mirrors `ConstCache::pi`: a few bits for the term count, plus a fixed margin.
    let guard_bits = n.bit_len() + 32;
    let work_precision = precision_for_bits::<B>(bits + guard_bits);
    let work = Context::<R>::new(work_precision);
    let num = work.convert_int::<B>(q.as_ibig() + &t).value();
    let den = work.convert_int::<B>(IBig::from(q)).value();
    let e = num / den;
    e.with_precision(precision)
}

impl fmt::Debug for ConstCache {
    fn fmt(&self, f: &mut fmt::Formatter<'_>) -> fmt::Result {
        f.debug_struct("ConstCache")
            .field("pi", &DebugSlot(&self.pi))
            .field("iacoth_6", &DebugSlot(&self.iacoth_6))
            .field("iacoth_9", &DebugSlot(&self.iacoth_9))
            .field("iacoth_99", &DebugSlot(&self.iacoth_99))
            .finish()
    }
}

/// Newtype so we can implement `Debug` for `&Option<CachedState>` via the
/// big-integer `Debug` formatters.
struct DebugSlot<'a>(&'a Option<CachedState>);

impl fmt::Debug for DebugSlot<'_> {
    fn fmt(&self, f: &mut fmt::Formatter<'_>) -> fmt::Result {
        match self.0 {
            Some(s) => f
                .debug_struct("CachedState")
                .field("num_terms", &s.num_terms)
                .field("p", &s.p)
                .field("q", &s.q)
                .field("t", &s.t)
                .finish(),
            None => f.write_str("None"),
        }
    }
}

#[cfg(test)]
mod tests {
    use super::*;
    use crate::round::mode;
    use crate::DBig;
    use alloc::format;
    use alloc::string::ToString; // not in the prelude under no_std

    /// Independently (left-fold) merge the first `k` leaves of `L(n)` and check
    /// that the result matches the precomputed constant triple. Guards against
    /// transcription errors in the `IACOTH_*_INITIAL` literals; correctness of the
    /// finalized `L(n)` values is covered by `log::tests`.
    #[test]
    fn test_iacoth_initial_blocks() {
        fn check(n: u32, expected: &(usize, UBig, UBig, IBig)) {
            let k = expected.0;
            // independently (left-fold) merge the first k leaves
            let mut acc = (UBig::ONE, UBig::ONE, IBig::ZERO);
            for kk in 1..=k {
                let pa = UBig::from(2 * kk as u64 - 1);
                let qa = UBig::from(2 * kk as u64 + 1) * UBig::from(n).pow(2);
                let ta = IBig::from(pa.clone());
                let (p, q, t) = merge(&acc.0, &acc.1, &acc.2, &pa, &qa, &ta);
                acc = (p, q, t);
            }
            assert_eq!(acc.0, expected.1, "P mismatch for n={n}");
            assert_eq!(acc.1, expected.2, "Q mismatch for n={n}");
            assert_eq!(acc.2, expected.3, "T mismatch for n={n}");
            // the seed branch must reproduce the same state as the full recursion
            assert_eq!(iacoth_bs(n, 1, 1 + k), (acc.0, acc.1, acc.2));
        }
        check(6, &IACOTH_6_INITIAL);
        check(9, &IACOTH_9_INITIAL);
        check(99, &IACOTH_99_INITIAL);
    }

    #[test]
    fn test_pi_matches_context() {
        // Cache miss must reproduce Context::pi exactly.
        for &precision in &[10usize, 50, 100] {
            let mut cache = ConstCache::new();
            let cached = cache.pi::<10, mode::HalfEven>(precision).value();
            let direct = Context::<mode::HalfEven>::new(precision)
                .pi::<10>(None)
                .value();
            assert_eq!(cached, direct, "pi mismatch at precision {precision}");
        }
    }

    #[test]
    fn test_pi_lower_precision_reuses() {
        // Compute at high precision, then a lower-precision request must round
        // down from the cached state and still be correct.
        let mut cache = ConstCache::new();
        let _pi_high = cache.pi::<10, mode::HalfEven>(200).value();
        // the slot now holds >=200 terms; a 50-digit request reuses it
        let pi_50 = cache.pi::<10, mode::HalfEven>(50).value();
        let direct = Context::<mode::HalfEven>::new(50).pi::<10>(None).value();
        assert_eq!(pi_50, direct);
    }

    #[test]
    fn test_pi_extension_matches_scratch() {
        // Extending 100 -> 1000 must be bit-identical to a from-scratch 1000-digit compute.
        let mut cache = ConstCache::new();
        let _pi_100 = cache.pi::<10, mode::HalfAway>(100).value();
        let pi_1000_extended = cache.pi::<10, mode::HalfAway>(1000).value();

        let direct = Context::<mode::HalfAway>::new(1000).pi::<10>(None).value();
        assert_eq!(pi_1000_extended, direct);
    }

    #[test]
    fn test_iacoth_matches_context() {
        let mut cache = ConstCache::new();
        // ln2 / ln10 via cache must match ln(2)/ln(10) computed independently
        // through Context::ln (a different, atanh-based algorithm) at several precisions.
        for &precision in &[20usize, 45, 80] {
            let cached_ln2 = cache
                .ln2::<10, mode::Zero>(precision)
                .with_precision(precision)
                .value();
            let ln2_ctx = Context::<mode::Zero>::new(precision);
            let direct_ln2 = ln2_ctx.unwrap_fp(ln2_ctx.ln::<10>(&Repr::new(2.into(), 0), None));
            assert_eq!(cached_ln2, direct_ln2, "ln2 mismatch at precision {precision}");

            let cached_ln10 = cache
                .ln10::<10, mode::Zero>(precision)
                .with_precision(precision)
                .value();
            let ln10_ctx = Context::<mode::Zero>::new(precision);
            let direct_ln10 = ln10_ctx.unwrap_fp(ln10_ctx.ln::<10>(&Repr::new(10.into(), 0), None));
            assert_eq!(cached_ln10, direct_ln10, "ln10 mismatch at precision {precision}");
        }
    }

    #[test]
    fn test_iacoth_extension_matches_scratch() {
        // Extend ln2 from low to high precision; result must match from-scratch.
        let mut cache = ConstCache::new();
        let _ln2_low = cache.ln2::<10, mode::HalfAway>(20);
        let ln2_high = cache.ln2::<10, mode::HalfAway>(120);

        let mut fresh = ConstCache::new();
        let direct = fresh.ln2::<10, mode::HalfAway>(120);
        assert_eq!(ln2_high, direct);
    }

    #[test]
    fn test_ln_base() {
        // binary base: ln(base) == ln(2)
        let mut cache = ConstCache::new();
        let ln_base = cache.ln_base::<2, mode::HalfAway>(50);
        let ln2 = cache.ln2::<2, mode::HalfAway>(50);
        assert_eq!(ln_base, ln2);

        // power-of-two base: ln(8) = 3·ln(2)
        let ln8 = cache.ln_base::<8, mode::HalfAway>(50);
        let expected = 3u8 * cache.ln2::<8, mode::HalfAway>(50);
        assert_eq!(ln8.with_precision(50).value(), expected.with_precision(50).value());
    }

    #[test]
    fn test_debug_shows_bigint_head_tail() {
        let mut cache = ConstCache::new();
        let _pi = cache.pi::<10, mode::HalfAway>(100); // populate the cache (value unused)
        let s = format!("{:?}", cache);
        assert!(s.contains("pi"));
        assert!(s.contains("num_terms"));
        // UBig/IBig Debug prints head..tail, so the output stays compact
        assert!(s.contains(".."), "Debug output should use head..tail truncation");
        assert!(s.len() < 512);
    }

    #[test]
    fn test_sqrt_10005_cached_and_counted() {
        // Computing π caches the base-free √10005 isqrt; total_words counts it, and
        // clear() frees it.
        let mut cache = ConstCache::new();
        assert_eq!(cache.total_terms(), 0);
        assert_eq!(cache.total_words(), 0);

        let _pi = cache.pi::<10, mode::HalfAway>(200); // fills the cache (value unused)
                                                       // the isqrt is now cached (total_terms stays series-only; words include isqrt)
        assert!(cache.total_words() > 0);

        cache.clear();
        assert_eq!(cache.total_terms(), 0);
        assert_eq!(cache.total_words(), 0);

        // after clear, π recomputes from scratch and still matches the direct value
        let direct = Context::<mode::HalfAway>::new(50).pi::<10>(None).value();
        let after_clear = cache.pi::<10, mode::HalfAway>(50).value();
        assert_eq!(after_clear, direct);
    }

    #[test]
    fn test_sqrt_10005_reuse_higher_precision() {
        // A high-precision π call caches a high-bit isqrt; a later lower-precision
        // call must reuse it (no recompute) and still be correct.
        let mut cache = ConstCache::new();
        let _high = cache.pi::<2, mode::HalfEven>(1000);
        let words_after_high = cache.total_words();

        let low = cache.pi::<2, mode::HalfEven>(100).value();
        // word count unchanged ⇒ isqrt (and series) were reused, not recomputed
        assert_eq!(cache.total_words(), words_after_high);

        let direct = Context::<mode::HalfEven>::new(100).pi::<2>(None).value();
        assert_eq!(low, direct);
    }

    /// `e_bs(1, N+1)` produces `(P, Q, T)` with `T/Q = Σ_{k=1}^{N} 1/k!` and `P ≡ 1`.
    #[test]
    fn test_e_bs_partial_sums() {
        // N=1: 1/1! = 1  →  (1, 1, 1)
        let (p, q, t) = e_bs(1, 2);
        assert_eq!(p, UBig::ONE);
        assert_eq!(q, UBig::from(1u32));
        assert_eq!(t, IBig::from(1));

        // N=3: 1 + 1/2 + 1/6 = 10/6  →  (1, 6, 10)
        let (p, q, t) = e_bs(1, 4);
        assert_eq!(p, UBig::ONE);
        assert_eq!(q, UBig::from(6u32));
        assert_eq!(t, IBig::from(10));

        // N=4: +1/24 → 41/24  →  (1, 24, 41)
        let (p, q, t) = e_bs(1, 5);
        assert_eq!(p, UBig::ONE);
        assert_eq!(q, UBig::from(24u32));
        assert_eq!(t, IBig::from(41));

        // empty range is the identity triple
        let (p, q, t) = e_bs(7, 7);
        assert_eq!(p, UBig::ONE);
        assert_eq!(q, UBig::ONE);
        assert_eq!(t, IBig::ZERO);
    }

    /// `compute_e` (factorial binary splitting) and `exp(1)` (argument reduction,
    /// Taylor series, and powering) are two independent algorithms for the same
    /// value, so both must yield the identical correctly-rounded result.
    #[test]
    fn test_e_matches_exp_one() {
        fn check<const B: Word>(p: usize) {
            let ctx = Context::<mode::HalfAway>::new(p);
            let e_const = ctx.e::<B>().value();
            let exp_one = ctx
                .exp::<B>(&Repr::<B>::new(IBig::from(1), 0), None)
                .unwrap()
                .value();
            assert_eq!(e_const, exp_one, "e != exp(1) at precision {p}, base {B}");
        }
        for &p in &[1usize, 2, 5, 13, 50, 137, 500] {
            check::<10>(p);
            check::<2>(p);
        }
    }

    /// Ground-truth prefix of e (OEIS A001113), independent of the `exp` path.
    /// The 45-decimal prefix is far from the precision-100 rounding boundary.
    #[test]
    fn test_e_known_decimal_prefix() {
        // e = 2.71828182845904523536028747135266249775724709369995957…
        let s = DBig::e(100).to_string();
        assert!(s.starts_with("2.718281828459045235360287471352662497757247093699"), "got {s}");
        // sanity: e is in (2, 3)
        let e = DBig::e(10);
        assert!(e > DBig::from(2u32));
        assert!(e < DBig::from(3u32));
    }

    /// `e_term_count` must return enough terms that Stirling's lower bound on
    /// `(N+1)!` actually exceeds the target bit count.
    #[test]
    fn test_e_term_count_sufficient() {
        for &bits in &[1usize, 10, 100, 1000, 10_000, 100_000] {
            let n = e_term_count(bits);
            let m = (n + 1) as f64;
            let stirling_lb = m * (m.log2() - core::f64::consts::LOG2_E);
            assert!(stirling_lb > bits as f64, "N={n} insufficient for bits={bits}");
        }
    }

    /// π must round correctly at every precision — the Chudnovsky series accuracy
    /// headroom must keep the truncation error well below 1 ulp.
    #[test]
    fn test_pi_no_off_by_one_ulp() {
        let pi_ref = FBig::<mode::Zero, 2>::pi(2048);
        for prec in 10..=1024 {
            let mid = Context::<mode::Zero>::new(prec).pi::<2>(None).value();
            let correct = pi_ref.clone().with_precision(prec).value();
            assert_eq!(mid.into_repr(), correct.into_repr(), "pi off by 1 ulp at prec={prec}");
        }
    }
}