dashu_float/

log.rs

use dashu_base::{
    utils::{next_down, next_up},
    Abs, AbsOrd,
    Approximation::*,
    EstimatedLog2, PowerOfTwo, Sign, UnsignedAbs,
};
use dashu_int::{IBig, UBig};

use crate::{
    ball::{ceil_shift, Ball},
    error::{assert_finite, assert_limited_precision, FpError, FpResult},
    fbig::FBig,
    math::cache::{reborrow_cache, ConstCache},
    repr::{Context, Repr, Word},
    round::{mode, ErrorBounds, Round},
};
use core::cmp::Ordering;

impl<const B: Word> EstimatedLog2 for Repr<B> {
    // currently a Word has at most 64 bits, so log2() < f32::MAX
    fn log2_bounds(&self) -> (f32, f32) {
        if self.significand.is_zero() {
            return (f32::NEG_INFINITY, f32::NEG_INFINITY);
        }

        // log(s*B^e) = log(s) + e*log(B)
        let (logs_lb, logs_ub) = self.significand.log2_bounds();
        let (logb_lb, logb_ub) = if B.is_power_of_two() {
            let log = B.trailing_zeros() as f32;
            (log, log)
        } else {
            B.log2_bounds()
        };
        let e = self.exponent as f32;
        let (lb, ub) = if self.exponent >= 0 {
            (logs_lb + e * logb_lb, logs_ub + e * logb_ub)
        } else {
            (logs_lb + e * logb_ub, logs_ub + e * logb_lb)
        };
        (next_down(lb), next_up(ub))
    }

    fn log2_est(&self) -> f32 {
        let logs = self.significand.log2_est();
        let logb = if B.is_power_of_two() {
            B.trailing_zeros() as f32
        } else {
            B.log2_est()
        };
        logs + self.exponent as f32 * logb
    }
}

impl<R: Round, const B: Word> EstimatedLog2 for FBig<R, B> {
    #[inline]
    fn log2_bounds(&self) -> (f32, f32) {
        self.repr.log2_bounds()
    }

    #[inline]
    fn log2_est(&self) -> f32 {
        self.repr.log2_est()
    }
}

impl<R: ErrorBounds, const B: Word> FBig<R, B> {
    /// Calculate the natural logarithm function (`log(x)`) on the float number.
    ///
    /// # Examples
    ///
    /// ```
    /// # use core::str::FromStr;
    /// # use dashu_base::ParseError;
    /// # use dashu_float::DBig;
    /// let a = DBig::from_str("1.234")?;
    /// assert_eq!(a.ln(), DBig::from_str("0.2103")?);
    /// # Ok::<(), ParseError>(())
    /// ```
    #[inline]
    pub fn ln(&self) -> Self {
        self.context.unwrap_fp(self.context.ln(&self.repr, None))
    }

    /// Calculate the natural logarithm function (`log(x+1)`) on the float number
    ///
    /// # Examples
    ///
    /// ```
    /// # use core::str::FromStr;
    /// # use dashu_base::ParseError;
    /// # use dashu_float::DBig;
    /// let a = DBig::from_str("0.1234")?;
    /// assert_eq!(a.ln_1p(), DBig::from_str("0.11636")?);
    /// # Ok::<(), ParseError>(())
    /// ```
    #[inline]
    pub fn ln_1p(&self) -> Self {
        self.context.unwrap_fp(self.context.ln_1p(&self.repr, None))
    }

    /// Calculate the base-2 logarithm (`log2(x)`) on the float number.
    ///
    /// Correctly rounded to the context's precision under any rounding mode. For an exact power
    /// of two the result is the exact integer `log2(x)`.
    ///
    /// # Examples
    ///
    /// ```
    /// # use core::str::FromStr;
    /// # use dashu_base::ParseError;
    /// # use dashu_float::DBig;
    /// let a = DBig::from_str("8")?;
    /// assert_eq!(a.log2(), DBig::from_str("3")?);
    /// # Ok::<(), ParseError>(())
    /// ```
    #[inline]
    pub fn log2(&self) -> Self {
        self.context.unwrap_fp(self.context.log2(&self.repr, None))
    }

    /// Calculate the base-10 logarithm (`log10(x)`) on the float number.
    ///
    /// Correctly rounded to the context's precision under any rounding mode. For an exact power
    /// of ten the result is the exact integer `log10(x)`.
    ///
    /// # Examples
    ///
    /// ```
    /// # use core::str::FromStr;
    /// # use dashu_base::ParseError;
    /// # use dashu_float::DBig;
    /// let a = DBig::from_str("1000")?;
    /// assert_eq!(a.log10(), DBig::from_str("3")?);
    /// # Ok::<(), ParseError>(())
    /// ```
    #[inline]
    pub fn log10(&self) -> Self {
        self.context.unwrap_fp(self.context.log10(&self.repr, None))
    }
}

// `ln2`/`ln10`/`iacoth`/`ln_base`/`ln_compute` are the near-correct logarithm primitives: they
// evaluate the series at a working precision and round once, without a Ziv certification step.
// They live on `R: Round` so that base conversion (`with_base_and_precision`, which only needs a
// near-correct constant `ln(B)`) can use them without inheriting the `ErrorBounds` bound. The
// correctly-rounded public `ln`/`ln_1p` (in the `ErrorBounds` impl below) wrap `ln_compute` in a
// Ziv loop.
impl<R: Round> Context<R> {
    /// Calculate log(2)
    ///
    /// The precision of the output will be larger than self.precision
    #[inline]
    fn ln2<const B: Word>(&self, cache: Option<&mut ConstCache>) -> FBig<R, B> {
        if let Some(c) = cache {
            return c.ln2::<B, R>(self.precision);
        }
        // log(2) = 4L(6) + 2L(99)
        // see formula (24) from Gourdon, Xavier, and Pascal Sebah.
        // "The Logarithmic Constant: Log 2." (2004)
        4 * self.iacoth(6.into()) + 2 * self.iacoth(99.into())
    }

    /// Calculate log(10)
    ///
    /// The precision of the output will be larger than self.precision
    #[inline]
    fn ln10<const B: Word>(&self, cache: Option<&mut ConstCache>) -> FBig<R, B> {
        if let Some(c) = cache {
            return c.ln10::<B, R>(self.precision);
        }
        // log(10) = log(2) + log(5) = 3log(2) + 2L(9)
        3 * self.ln2(None) + 2 * self.iacoth(9.into())
    }

    /// Calculate log(B), for internal use only
    ///
    /// The precision of the output will be larger than self.precision
    #[inline]
    pub(crate) fn ln_base<const B: Word>(&self, cache: Option<&mut ConstCache>) -> FBig<R, B> {
        if let Some(c) = cache {
            return c.ln_base::<B, R>(self.precision);
        }
        match B {
            2 => self.ln2(None),
            10 => self.ln10(None),
            i if i.is_power_of_two() => self.ln2(None) * i.trailing_zeros(),
            _ => {
                // Near-correct ln(B) via the atanh series (no Ziv certification — base conversion
                // only needs a near-correct constant). `ln_compute` is on `R: Round`, so this keeps
                // `ln_base` callable from `R: Round` contexts (base conversion).
                let guard = self.base_guard_digits::<B>() + 2;
                self.ln_compute::<B>(
                    &Repr::new(Repr::<B>::BASE.into(), 0),
                    self.precision + guard,
                    false,
                    None,
                )
                .to_value_radius::<R>()
                .0
            }
        }
    }

    /// `ln(B)` as a [`Ball`], carrying the mechanical radius.
    ///
    /// The cached bases (2, 10, powers of 2) evaluate correctly-rounded constants (error ≤ ½ ulp at
    /// the working precision), so the fixed `8` ulps is a sound loose bound. A generic base falls
    /// back to `ln_compute`'s atanh series, whose radius is ~(series terms + B) ulps — far larger
    /// than 8 — so its mechanical radius is kept instead (a hard-coded `8` would under-bind the
    /// `s·ln(B)` reconstruction error in `exp_compute`'s reduction).
    pub(crate) fn ln_base_ball<const B: Word>(
        &self,
        mut cache: Option<&mut ConstCache>,
    ) -> Ball<B> {
        let ctx = Context::<mode::HalfEven>::new(self.precision);
        match B {
            10 => {
                let logb = ctx.ln_base::<B>(reborrow_cache(&mut cache));
                Ball::with_error(logb, IBig::from(8))
            }
            i if i.is_power_of_two() => {
                let logb = ctx.ln_base::<B>(reborrow_cache(&mut cache));
                Ball::with_error(logb, IBig::from(8))
            }
            _ => {
                // Generic base: no cached sub-series applies, so compute ln(B) directly at the
                // requested precision and keep `ln_compute`'s ball error.
                ctx.ln_compute::<B>(
                    &Repr::new(Repr::<B>::BASE.into(), 0),
                    self.precision,
                    false,
                    reborrow_cache(&mut cache),
                )
            }
        }
    }

    /// Calculate L(n) = acoth(n) = atanh(1/n) = 1/2 log((n+1)/(n-1)), given by the
    /// series
    ///
    /// ```text
    ///                1     n + 1              1
    ///   atanh(1/n) = — log(—————) = Σ   ——————————————————
    ///                2     n - 1   i≥0 n^(2i+1) · (2i+1)
    /// ```
    ///
    /// This method is intended to be used in logarithm calculation,
    /// so the precision of the output will be larger than desired precision.
    ///
    /// Evaluated by binary splitting (see [`iacoth_bs`][crate::math::cache::iacoth_bs]):
    /// the exact integer tree state `(P, Q, T)` over `[1, N)` satisfies
    /// `L(n) = (Q + T)/(n·Q)`, with `Q` kept at O(p) digits by the ratio-form
    /// term recurrence.
    fn iacoth<const B: Word>(&self, n: IBig) -> FBig<R, B> {
        let n: u32 = (&n).try_into().expect("iacoth argument must fit in u32");

        // number of series terms until r_k < B^{-p}:  (2k+1)·log_B(n) > p.
        // The count is generously over-provisioned, so a truncating cast stands in
        // for a ceiling.
        let log_b_n = n.log2_est() / B.log2_est();
        let num_terms = (self.precision as f32 / (2.0 * log_b_n)) as usize + 10;

        let (_p, q, t) = crate::math::cache::iacoth_bs(n, 1, num_terms + 1);

        // L(n) = (Q + T) / (n·Q). Extra guard digits absorb the division's rounding
        // (the binary-splitting state is exact, so only this single round loses anything).
        let guard_digits = self.base_guard_digits::<B>();
        let work_context = Self::new(self.precision + guard_digits + 2);

        let num = work_context.convert_int::<B>(q.as_ibig() + &t).value();
        let denom = work_context.convert_int::<B>(IBig::from(n) * &q).value();
        num / denom
    }

    /// Evaluate `ln(x)` (or `ln(x+1)` when `one_plus`) at `work_precision` via the atanh series,
    /// returning a [`Ball`] whose radius is derived mechanically by Ball arithmetic (error
    /// propagates term-by-term through the series; cancellation and the `s·ln(B)` reconstruction
    /// flow through the Ball scale factors).
    ///
    /// This is the near-correct computation core shared by the public Ziv-backed `ln`/`ln_1p`
    /// (which wrap it in a retry loop) and by `ln_base` (which only needs a near-correct constant
    /// `ln(B)`). It lives on `R: Round` so those near-correct callers don't inherit the
    /// `ErrorBounds` bound.
    pub(crate) fn ln_compute<const B: Word>(
        &self,
        x: &Repr<B>,
        mut work_precision: usize,
        one_plus: bool,
        mut cache: Option<&mut ConstCache>,
    ) -> Ball<B> {
        // Round the input to the working precision; the input's own rounding is the only error
        // introduced here.
        let context = Context::<mode::HalfEven>::new(work_precision);
        let x_ball = Ball::from_rounded(context.repr_round_ref(x).map(|r| FBig::new(r, context)));

        // When one_plus is true and |x| < 1/B, the input is fed into the Maclaurin without scaling.
        let no_scaling = one_plus && x_ball.mid.log2_est() < -B.log2_est();

        let (s, mut x_scaled) = if no_scaling {
            (0, x_ball)
        } else {
            let x_ball = if one_plus {
                x_ball.add(&Ball::exact_int(work_precision, IBig::ONE))
            } else {
                x_ball
            };

            let log2 = x_ball.mid.log2_bounds().0;
            let s = log2 as isize - (log2 < 0.) as isize; // floor(log2(x))

            let mut exact = x_ball.n.is_zero();
            let x_scaled = if B == 2 {
                x_ball.shift(s) // exact (power-of-base shift)
            } else if s > 0 {
                // Exact divisor 2^s: the error shrinks by it directly (no general-division rational).
                x_ball.div_exact(&(IBig::ONE << s as usize))
            } else {
                // Scaling by 2^|s| is exact (finite decimal × power of two). Use the tracking
                // variant so an exact operand keeps n = 0 — otherwise the unconditional `+1`
                // would be amplified by `rescale_precision` into a `B^precision`-sized error
                // count that never shrinks across Ziv retries (the powf-of-base-<1 hang).
                x_ball.scale_int_tracking(&(IBig::ONE << (-s) as usize), &mut exact)
            };
            debug_assert!(x_scaled.mid >= FBig::<mode::HalfEven, B>::ONE);
            (s, x_scaled)
        };

        // The reconstruction 2·sum + s·ln(B) *cancels* for x < 1 (s < 0), so the series runs at
        // double precision to keep the pre-cancellation sum accurate. The finer ulp rescales `n`.
        if s < 0 || x_scaled.mid.repr().sign() == Sign::Negative {
            work_precision += self.precision;
            x_scaled.rescale_precision(self.precision);
        }
        let work_context = Context::<mode::HalfEven>::new(work_precision);

        // after the number is scaled to nearly one, use Maclaurin series on log(x) = 2atanh(z):
        // let z = (x-1)/(x+1) < 1, log(x) = 2atanh(z) = 2Σ(z²ⁱ⁺¹/(2i+1)) for i = 1,3,5,...
        let z = if no_scaling {
            let two = Ball::exact_int(work_precision, IBig::from(2));
            let den = x_scaled.add(&two);
            x_scaled.div(&den)
        } else {
            let one = Ball::exact_int(work_precision, IBig::ONE);
            let num = x_scaled.sub(&one);
            let den = x_scaled.add(&one);
            num.div(&den)
        };
        let z2 = z.mul(&z);
        let mut pow = z.clone();
        let mut sum = z;
        let mut k: usize = 3;
        loop {
            pow = pow.mul(&z2);

            let increase = pow.div_int(k);
            if increase.mid.abs_cmp(&sum.mid.ulp_lb()).is_le() {
                break;
            }

            sum = sum.add(&increase);
            k += 2;
        }

        // Omitted series tail: the first omitted term is ≤ sum.ulp_lb(), and the tail of the
        // atanh series shrinks by z² per step with 1/(1−z²) < B for x_scaled ∈ [1, B), so the
        // tail is < B·sum.ulp_lb() < B ulps of sum.
        sum.inflate(&IBig::from(B));

        // compose the logarithm of the original number
        let sum2 = sum.scale_int(&IBig::from(2));
        if no_scaling {
            sum2
        } else {
            // ln(2) as a ball. The constant evaluates the atanh series via binary splitting at
            // work + guard digits and rounds once to `work_precision`, so its error is a handful
            // of work-precision ulps; 8 is a conservative sound bound for every code path
            // (cached and uncached).
            let ln2 = work_context.ln2::<B>(reborrow_cache(&mut cache));
            let ln2 = Ball::with_error(ln2, IBig::from(8));
            sum2.add(&ln2.scale_int(&IBig::from(s)))
        }
    }

    /// `ln(1 + arg)` of a *ball* input. [`ln_compute`](Self::ln_compute) evaluates the series on
    /// `arg.mid`; the input ball's own error `|θ| ≤ arg.n·ulp(arg)` then contributes `|θ|/(1+arg)`
    /// to the log. Bound via `(1+arg)`'s ball magnitude: for a mostly-correct argument the factor
    /// `1/(1−|θ|/(1+arg)) ≤ 2` is sound, so the adjustment is `⌈2·n_arg·ulp_arg/((1+arg)·ulp_ln)⌉`.
    pub(crate) fn ln_1p_ball<const B: Word>(
        &self,
        arg: &Ball<B>,
        mut cache: Option<&mut ConstCache>,
    ) -> Ball<B> {
        let mut ln_ball =
            self.ln_compute::<B>(arg.mid.repr(), self.precision, true, reborrow_cache(&mut cache));
        let den = arg.add(&Ball::exact_int(self.precision, IBig::ONE));
        let e_d = den.mid.repr().exponent; // (1+arg) = sig_d·B^(e_d)
        let sig_d = den.mid.repr().significand.clone().abs();
        // `lead_*` is the leading position (`lead_exp`), so `ulp_arg = B^(lead_arg − p_arg)` and
        // `ulp_ln = B^(lead_ln − p_ln)`. The input error propagates as
        //   n_arg·ulp_arg / ((1+arg)·ulp_ln) = n_arg·B^(lead_arg − p_arg − e_d − lead_ln + p_ln)/sig_d;
        // ×2 for the 1/(1−|θ|/(1+arg)) factor.
        // The precision difference is essential: `ln_compute`'s s<0 path runs at double precision,
        // so `ln_ball` sits at 2·self.precision while `arg` stays at self.precision — dropping the
        // `−p_arg+p_ln` term under-bounds the adjust by B^(p_ln−p_arg) (atanh(x<0) near the pole
        // then mis-certifies, e.g. off by 2^13 ulps).
        let lead_arg = Ball::lead_exp(&arg.mid);
        let p_arg = arg.mid.precision();
        let lead_ln = Ball::lead_exp(&ln_ball.mid);
        let p_ln = ln_ball.mid.precision();
        let shift = lead_arg - p_arg as isize - e_d - lead_ln + p_ln as isize;
        let num = ceil_shift::<B>(2 * &arg.n, shift);
        let adjust = (num + &sig_d - IBig::ONE) / sig_d;
        ln_ball.inflate(&adjust);
        ln_ball
    }
}

// `ln`/`ln_1p` are correctly rounded via the Ziv loop, whose containment test needs the rounding
// preimage (`R: ErrorBounds`). They delegate the series to `ln_compute`.
impl<R: ErrorBounds> Context<R> {
    /// Calculate the natural logarithm function (`log(x)`) on the float number under this context.
    ///
    /// # Examples
    ///
    /// ```
    /// # use core::str::FromStr;
    /// # use dashu_base::ParseError;
    /// # use dashu_float::DBig;
    /// use dashu_base::Approximation::*;
    /// use dashu_float::{Context, round::{mode::HalfAway, Rounding::*}};
    ///
    /// let context = Context::<HalfAway>::new(2);
    /// let a = DBig::from_str("1.234")?;
    /// assert_eq!(context.ln(&a.repr(), None), Ok(Inexact(DBig::from_str("0.21")?, NoOp)));
    /// # Ok::<(), ParseError>(())
    /// ```
    #[inline]
    pub fn ln<const B: Word>(
        &self,
        x: &Repr<B>,
        cache: Option<&mut ConstCache>,
    ) -> FpResult<FBig<R, B>> {
        if x.is_infinite() {
            return Err(FpError::InfiniteInput);
        }
        if x.significand.is_zero() {
            // ln(±0) = -inf (a value, not an error)
            return Ok(Exact(FBig::new(Repr::neg_infinity(), *self)));
        }
        if x.sign() == Sign::Negative {
            return Err(FpError::OutOfDomain);
        }
        self.ln_internal(x, false, cache)
    }

    /// Calculate the natural logarithm function (`log(x+1)`) on the float number under this context.
    ///
    /// # Examples
    ///
    /// ```
    /// # use core::str::FromStr;
    /// # use dashu_base::ParseError;
    /// # use dashu_float::DBig;
    /// use dashu_base::Approximation::*;
    /// use dashu_float::{Context, round::{mode::HalfAway, Rounding::*}};
    ///
    /// let context = Context::<HalfAway>::new(2);
    /// let a = DBig::from_str("0.1234")?;
    /// assert_eq!(context.ln_1p(&a.repr(), None), Ok(Inexact(DBig::from_str("0.12")?, AddOne)));
    /// # Ok::<(), ParseError>(())
    /// ```
    #[inline]
    pub fn ln_1p<const B: Word>(
        &self,
        x: &Repr<B>,
        cache: Option<&mut ConstCache>,
    ) -> FpResult<FBig<R, B>> {
        if x.is_infinite() {
            return Err(FpError::InfiniteInput);
        }
        // Domain of ln_1p is x > -1. x == -1 gives -inf; x < -1 is out of domain.
        if x.sign() == Sign::Negative && !x.significand.is_zero() {
            match FBig::<R, B>::new(x.clone(), *self).abs_cmp(&FBig::ONE) {
                Ordering::Greater => return Err(FpError::OutOfDomain), // x < -1
                Ordering::Equal => return Ok(Exact(FBig::new(Repr::neg_infinity(), *self))),
                _ => {}
            }
        }
        self.ln_internal(x, true, cache)
    }

    fn ln_internal<const B: Word>(
        &self,
        x: &Repr<B>,
        one_plus: bool,
        mut cache: Option<&mut ConstCache>,
    ) -> FpResult<FBig<R, B>> {
        assert_finite(x);

        // Exact special cases first: they need no rounding, so a precision-0 (unlimited)
        // value such as `FBig::ONE` or the one from `try_from(0.0)` must still resolve
        // ln/ln_1p exactly rather than tripping the limited-precision assertion below.
        if !one_plus && x.is_one() {
            return Ok(Exact(FBig::ZERO)); // ln(1) = +0
        }
        if one_plus && x.significand.is_zero() {
            // ln_1p(±0) = ±0
            let zero = if x.is_neg_zero() {
                FBig::new(Repr::neg_zero(), *self)
            } else {
                FBig::ZERO
            };
            return Ok(Exact(zero));
        }

        assert_limited_precision(self.precision);

        // Correct rounding via the Ziv loop: `ln_compute` evaluates the atanh series at `p + guard`
        // and reports a provable error radius; the driver retries with more guard digits until the
        // approximation's error interval lies entirely inside one rounding bin. The guard is a
        // *performance* knob (first-attempt hit rate), not a correctness backstop — Ziv certifies
        // the result. (The pre-Ziv `+ 2` is retained: with the conservative radius below it is still
        // needed for the first attempt to clear the half-ulp preimage at typical precisions.)
        let base_guard = self.base_guard_digits::<B>() + 2;
        self.ziv(base_guard + one_plus as usize, |guard| {
            Ok(self
                .ln_compute::<B>(x, self.precision + guard, one_plus, reborrow_cache(&mut cache))
                .to_value_radius::<R>())
        })
    }

    /// Calculate the base-2 logarithm (`log2(x)`) on the float number under this context.
    ///
    /// Correctly rounded to the context's precision under any rounding mode; for an exact power
    /// of two the result is the exact integer `log2(x)`.
    ///
    /// # Domain
    ///
    /// `log2(±0) = −∞` and a negative (non-zero) input is out of domain; an infinite input is an
    /// error (a finite context cannot produce the infinite `log2(+∞) = +∞` exactly).
    ///
    /// # Examples
    ///
    /// ```
    /// # use core::str::FromStr;
    /// # use dashu_base::ParseError;
    /// # use dashu_float::DBig;
    /// use dashu_base::Approximation::*;
    /// use dashu_float::{Context, round::{mode::HalfAway, Rounding::*}};
    ///
    /// let context = Context::<HalfAway>::new(4);
    /// let a = DBig::from_str("10")?;
    /// assert_eq!(context.log2(&a.repr(), None), Ok(Inexact(DBig::from_str("3.322")?, AddOne)));
    /// # Ok::<(), ParseError>(())
    /// ```
    #[inline]
    pub fn log2<const B: Word>(
        &self,
        x: &Repr<B>,
        cache: Option<&mut ConstCache>,
    ) -> FpResult<FBig<R, B>> {
        if x.is_infinite() {
            return Err(FpError::InfiniteInput);
        }
        if x.significand.is_zero() {
            // log2(±0) = -inf (a value, not an error)
            return Ok(Exact(FBig::new(Repr::neg_infinity(), *self)));
        }
        if x.sign() == Sign::Negative {
            return Err(FpError::OutOfDomain);
        }
        self.log2_internal(x, cache)
    }

    fn log2_internal<const B: Word>(
        &self,
        x: &Repr<B>,
        mut cache: Option<&mut ConstCache>,
    ) -> FpResult<FBig<R, B>> {
        assert_finite(x);

        // Exact shortcuts first — they also cover unlimited precision, which the Ziv loop below
        // rejects via its limited-precision assertion.
        if x.is_one() {
            return Ok(Exact(FBig::ZERO)); // log2(1) = +0
        }

        // Exact power-of-two shortcut: if x = 2^k for an integer k, log2(x) = k. This is *required*
        // for directed rounding — the Ziv loop below cannot certify an exactly-representable
        // result whose true value sits on a rounding boundary (its shrinking error interval
        // always straddles the boundary), so without this shortcut log2(2^-159) under `Up` would
        // exhaust the retry cap and return k + 1 ulp instead of the exact k.
        //
        // log2(x) = log2(significand) + exponent·log2(B). With significand = 2^m this is an exact
        // integer whenever log2(B) is integral (B a power of two), or — for a non-power-of-two
        // base — when the exponent is zero.
        let mag = (&x.significand).unsigned_abs();
        if mag.is_power_of_two() && (x.exponent == 0 || B.is_power_of_two()) {
            let m = mag.trailing_zeros().unwrap(); // = log2(significand)
            let log2_b = B.trailing_zeros() as isize;
            let k = IBig::from(m) + IBig::from(x.exponent) * IBig::from(log2_b);
            return Ok(self.convert_int::<B>(k));
        }

        assert_limited_precision(self.precision);

        // log2(x) = ln(x)/ln(2), correctly rounded via the Ziv loop. Both logarithms come from the
        // Ball-based `ln_compute`, and dividing them as Balls composes the radius mechanically:
        // the quotient's error is bounded from the two logarithms' relative errors, with no
        // directed-interval bookkeeping or guard-digit constant. The driver certifies the result
        // against the rounding preimage exactly as before.
        let initial_guard = self.base_guard_digits::<B>() + 4;
        self.ziv(initial_guard, |guard| {
            let work_precision = self.precision + guard;
            let lx = self.ln_compute::<B>(x, work_precision, false, reborrow_cache(&mut cache));
            let two = Repr::new(IBig::from(2), 0);
            let l2 = self.ln_compute::<B>(&two, work_precision, false, reborrow_cache(&mut cache));
            Ok(lx.div(&l2).to_value_radius::<R>())
        })
    }

    /// Calculate the base-10 logarithm (`log10(x)`) on the float number under this context.
    ///
    /// Correctly rounded to the context's precision under any rounding mode; for an exact power
    /// of ten the result is the exact integer `log10(x)`.
    ///
    /// # Domain
    ///
    /// `log10(±0) = −∞` and a negative (non-zero) input is out of domain; an infinite input is an
    /// error (a finite context cannot produce the infinite `log10(+∞) = +∞` exactly).
    ///
    /// # Examples
    ///
    /// ```
    /// # use core::str::FromStr;
    /// # use dashu_base::ParseError;
    /// # use dashu_float::DBig;
    /// use dashu_base::Approximation::*;
    /// use dashu_float::{Context, round::{mode::HalfAway, Rounding::*}};
    ///
    /// let context = Context::<HalfAway>::new(4);
    /// let a = DBig::from_str("100")?;
    /// assert_eq!(context.log10(&a.repr(), None), Ok(Exact(DBig::from_str("2")?)));
    /// # Ok::<(), ParseError>(())
    /// ```
    #[inline]
    pub fn log10<const B: Word>(
        &self,
        x: &Repr<B>,
        cache: Option<&mut ConstCache>,
    ) -> FpResult<FBig<R, B>> {
        if x.is_infinite() {
            return Err(FpError::InfiniteInput);
        }
        if x.significand.is_zero() {
            // log10(±0) = -inf (a value, not an error)
            return Ok(Exact(FBig::new(Repr::neg_infinity(), *self)));
        }
        if x.sign() == Sign::Negative {
            return Err(FpError::OutOfDomain);
        }
        self.log10_internal(x, cache)
    }

    fn log10_internal<const B: Word>(
        &self,
        x: &Repr<B>,
        mut cache: Option<&mut ConstCache>,
    ) -> FpResult<FBig<R, B>> {
        assert_finite(x);

        // Exact shortcuts first — they also cover unlimited precision, which the Ziv loop below
        // rejects via its limited-precision assertion.
        if x.is_one() {
            return Ok(Exact(FBig::ZERO)); // log10(1) = +0
        }

        // Exact power-of-ten shortcut: if x = 10^m for an integer m, log10(x) = m. This is *required*
        // for directed rounding — the Ziv loop below cannot certify an exactly-representable result
        // whose true value sits on a rounding boundary (its shrinking error interval always
        // straddles the one-sided preimage), so without this shortcut log10(10^-159) under `Up`
        // would exhaust the retry cap and return m + 1 ulp instead of the exact m.
        if let Some(m) = exact_pow10_log::<B>(&x.significand, x.exponent) {
            return Ok(self.convert_int::<B>(IBig::from(m)));
        }

        assert_limited_precision(self.precision);

        // log10(x) = ln(x)/ln(10), correctly rounded via the Ziv loop. Both logarithms come from the
        // Ball-based `ln_compute`, and dividing them as Balls composes the radius mechanically:
        // the quotient's error is bounded from the two logarithms' relative errors, with no
        // directed-interval bookkeeping or guard-digit constant. The driver certifies the result
        // against the rounding preimage exactly as before.
        let initial_guard = self.base_guard_digits::<B>() + 4;
        self.ziv(initial_guard, |guard| {
            let work_precision = self.precision + guard;
            let lx = self.ln_compute::<B>(x, work_precision, false, reborrow_cache(&mut cache));
            let ten = Repr::new(IBig::from(10), 0);
            let l10 = self.ln_compute::<B>(&ten, work_precision, false, reborrow_cache(&mut cache));
            Ok(lx.div(&l10).to_value_radius::<R>())
        })
    }
}

/// If `x = sig·B^e` is exactly `10^m` for some integer `m`, return `m`; otherwise `None`.
///
/// `10^m = 2^m·5^m`, so `x` is a power of ten iff the 2-valuation and 5-valuation of `sig·B^e`
/// coincide and `x` has no other prime factor. [`UBig::remove_word`] divides out all 2s and 5s from
/// both the base and the significand, returning each valuation as the removed exponent (and leaving
/// any non-{2,5} prime factor behind as a cofactor ≠ 1). The base `B = 2^p·5^q·s` (with `s` coprime
/// to 10) contributes `p·e` to the 2-valuation and `q·e` to the 5-valuation, and `s` must not
/// appear (unless `e = 0`). The valuations may be negative (`x` a negative power of ten).
fn exact_pow10_log<const B: Word>(sig: &IBig, e: isize) -> Option<isize> {
    // base: divide out all 2s and 5s, getting (p, q, leftover)
    let mut rest = UBig::from_word(B);
    let p = rest.remove_word(2)? as isize; // B ≥ 2, so never None
    let q = rest.remove_word(5).unwrap() as isize;
    // the leftover would give `x` a non-{2,5} prime factor when e ≠ 0
    if !rest.is_one() && e != 0 {
        return None;
    }

    // significand: divide out all 2s and 5s, getting (v2, v5, cofactor)
    let mut sig_abs = sig.unsigned_abs();
    let v2_sig = sig_abs.remove_word(2)? as isize; // non-zero upstream, so never None
    let v5 = sig_abs.remove_word(5).unwrap() as isize;
    // after removing every 2 and 5 the significand must be 1 (no other prime factor)
    if sig_abs != UBig::ONE {
        return None;
    }

    // `p*e` / `q*e` can overflow `isize` for an exotic base (B a large power of 2) at extreme
    // exponents on 32-bit; the true log10 would not fit `isize` then, so bail rather than compare
    // wrapped values that could falsely match.
    let v2 = v2_sig.checked_add(p.checked_mul(e)?)?;
    let v5 = v5.checked_add(q.checked_mul(e)?)?;
    (v2 == v5).then_some(v2)
}

#[cfg(test)]
mod tests {
    use super::*;
    use crate::round::mode;
    use alloc::vec::Vec;
    use dashu_base::BitTest;

    #[test]
    fn test_log10_domain() {
        let ctx = Context::<mode::HalfEven>::new(53);
        // log10(±0) = -inf (a value, not an error)
        let r = ctx.log10::<2>(&Repr::<2>::zero(), None).unwrap().value();
        assert!(r.repr.is_infinite());
        assert_eq!(r.repr.sign(), Sign::Negative);
        // log10(negative) is out of domain
        assert!(matches!(
            ctx.log10::<2>(&Repr::new((-1).into(), 0), None),
            Err(FpError::OutOfDomain)
        ));
        // an infinite input is rejected
        assert!(matches!(ctx.log10::<2>(&Repr::infinity(), None), Err(FpError::InfiniteInput)));
    }

    #[test]
    fn test_log10_exact_power_of_ten() {
        // log10(10^k) = k exactly under every rounding mode. Regression for the directed-rounding
        // defect the power-of-ten shortcut exists for: rounding ln(x) and ln(10) each toward the
        // mode and dividing once does not bound the quotient, so previously log10(10^-159) under
        // `Up` returned -159 + 1 ulp.
        let p = 53;
        for k in [0isize, 1, -1, 5, -159, 1000, -1000] {
            let x = Repr::<10>::new(IBig::from(1), k); // 10^k (base 10: significand 1)
            let r_down = Context::<mode::Down>::new(p)
                .log10::<10>(&x, None)
                .unwrap()
                .value();
            let r_up = Context::<mode::Up>::new(p)
                .log10::<10>(&x, None)
                .unwrap()
                .value();
            let r_zero = Context::<mode::Zero>::new(p)
                .log10::<10>(&x, None)
                .unwrap()
                .value();
            let r_he = Context::<mode::HalfEven>::new(p)
                .log10::<10>(&x, None)
                .unwrap()
                .value();
            assert_eq!(r_down.repr, r_he.repr, "Down != HalfEven for log10(10^{k})");
            assert_eq!(r_up.repr, r_he.repr, "Up != HalfEven for log10(10^{k})");
            assert_eq!(r_zero.repr, r_he.repr, "Zero != HalfEven for log10(10^{k})");
            assert_eq!(r_he.to_int().value(), IBig::from(k), "value for log10(10^{k})");
        }
    }

    #[test]
    fn test_exact_pow10_log_overflow() {
        // `p*e` overflows isize for a base that is a large power of two at an extreme exponent
        // (on 32-bit and 64-bit alike); the exact log wouldn't fit isize, so it must be None, not
        // a wrapped value that could falsely compare equal and return a wrong `Some`.
        assert_eq!(exact_pow10_log::<256>(&IBig::ONE, isize::MAX), None);
        assert_eq!(exact_pow10_log::<256>(&IBig::ONE, isize::MIN), None);
        // sanity: the normal path is unaffected.
        assert_eq!(exact_pow10_log::<10>(&IBig::ONE, 5), Some(5));
    }

    #[test]
    fn test_log10_exact_power_of_ten_binary_base() {
        // In base 2, 10^k is a float only via a 5^k-significand (e.g. 100 = 25·2^2); the
        // valuation-based shortcut must still detect the exact log10.
        let p = 53;
        for (sig, e, want) in [(25i32, 2isize, 2i64), (5, 1, 1), (125, 3, 3), (50, 1, 2)] {
            // 50·2^1 = 100 = 10^2 too (a non-normalized significand)
            let x = Repr::<2>::new(IBig::from(sig), e);
            let r_down = Context::<mode::Down>::new(p)
                .log10::<2>(&x, None)
                .unwrap()
                .value();
            let r_up = Context::<mode::Up>::new(p)
                .log10::<2>(&x, None)
                .unwrap()
                .value();
            let r_zero = Context::<mode::Zero>::new(p)
                .log10::<2>(&x, None)
                .unwrap()
                .value();
            let r_he = Context::<mode::HalfEven>::new(p)
                .log10::<2>(&x, None)
                .unwrap()
                .value();
            assert_eq!(r_down.to_int().value(), IBig::from(want), "log10({sig}·2^{e}) Down");
            assert_eq!(r_up.to_int().value(), IBig::from(want), "log10({sig}·2^{e}) Up");
            assert_eq!(r_zero.to_int().value(), IBig::from(want), "log10({sig}·2^{e}) Zero");
            assert_eq!(r_he.to_int().value(), IBig::from(want), "log10({sig}·2^{e}) HalfEven");
        }
    }

    #[test]
    fn test_log10_fbig_convenience() {
        // FBig::log10 convenience layer: exact powers of ten in both bases.
        let x = FBig::<mode::HalfEven, 10>::from_repr(Repr::new(1000.into(), 0), Context::new(50));
        assert_eq!(x.log10(), FBig::<mode::HalfEven, 10>::from(3u8));
        // base-2 float 100 = 25·2^2
        let x = FBig::<mode::HalfEven, 2>::from_repr(Repr::new(25.into(), 2), Context::new(50));
        assert_eq!(x.log10(), FBig::<mode::HalfEven, 2>::from(2u8));
    }

    /// Fixed inputs for the `log10` oracle differential.
    fn log10_diff_inputs() -> Vec<Repr<2>> {
        let mut v = Vec::new();
        for x in [0.5f64, 1.5, 2.0, 3.0, 10.0, 1000.0, 1e-6, 123.456, 2.5e-10] {
            v.push(FBig::<mode::HalfEven, 2>::try_from(x).unwrap().into_repr());
        }
        for k in [-100isize, -50, -10, -1, 0, 1, 10, 50, 100] {
            v.push(Repr::new(IBig::ONE, k)); // 2^k
        }
        v.push(
            FBig::<mode::HalfEven, 2>::try_from(f64::MAX)
                .unwrap()
                .into_repr(),
        );
        v
    }

    fn check_log10_differential<R: ErrorBounds>(p: usize, x: &Repr<2>, oracle: &Repr<2>) {
        let ctx = Context::<R>::new(p);
        let want = ctx.repr_round_ref(oracle).value();
        let got = ctx.log10_internal::<2>(x, None).unwrap().value();
        assert_eq!(got.repr, want, "p={p} {} x={x:?}", core::any::type_name::<R>(),);
    }

    #[test]
    fn log10_ball_matches_oracle() {
        let inputs = log10_diff_inputs();
        for p in [20usize, 50, 100] {
            for x in &inputs {
                let oracle = Context::<mode::HalfEven>::new(p + 60)
                    .log10::<2>(x, None)
                    .unwrap()
                    .value();
                check_log10_differential::<mode::HalfEven>(p, x, &oracle.repr);
                check_log10_differential::<mode::Down>(p, x, &oracle.repr);
                check_log10_differential::<mode::Up>(p, x, &oracle.repr);
                check_log10_differential::<mode::Zero>(p, x, &oracle.repr);
                check_log10_differential::<mode::Away>(p, x, &oracle.repr);
            }
        }
        // the arbitrary-precision regime: a reduced sweep (directed modes still exercised).
        for x in inputs.iter().step_by(9) {
            let oracle = Context::<mode::HalfEven>::new(560)
                .log10::<2>(x, None)
                .unwrap()
                .value();
            check_log10_differential::<mode::HalfEven>(500, x, &oracle.repr);
            check_log10_differential::<mode::Down>(500, x, &oracle.repr);
            check_log10_differential::<mode::Up>(500, x, &oracle.repr);
        }
    }

    #[test]
    fn test_ln_zero_is_neg_infinity() {
        let ctx = Context::<mode::HalfEven>::new(53);
        let r = ctx.ln::<2>(&Repr::<2>::zero(), None).unwrap().value();
        assert!(r.repr().is_infinite());
        assert_eq!(r.repr().sign(), Sign::Negative);
    }

    #[test]
    fn test_iacoth() {
        let context = Context::<mode::Zero>::new(10);
        let binary_6 = context.iacoth::<2>(6.into()).with_precision(10).value();
        assert_eq!(binary_6.repr.significand, IBig::from(689));
        let decimal_6 = context.iacoth::<10>(6.into()).with_precision(10).value();
        assert_eq!(decimal_6.repr.significand, IBig::from(1682361183));

        let context = Context::<mode::Zero>::new(40);
        let decimal_6 = context.iacoth::<10>(6.into()).with_precision(40).value();
        assert_eq!(
            decimal_6.repr.significand,
            IBig::from_str_radix("1682361183106064652522967051084960450557", 10).unwrap()
        );

        let context = Context::<mode::Zero>::new(201);
        let binary_6 = context.iacoth::<2>(6.into()).with_precision(201).value();
        assert_eq!(
            binary_6.repr.significand,
            IBig::from_str_radix(
                "2162760151454160450909229890833066944953539957685348083415205",
                10
            )
            .unwrap()
        );
    }

    #[test]
    fn test_ln2_ln10() {
        let context = Context::<mode::Zero>::new(45);
        let decimal_ln2 = context.ln2::<10>(None).with_precision(45).value();
        assert_eq!(
            decimal_ln2.repr.significand,
            IBig::from_str_radix("693147180559945309417232121458176568075500134", 10).unwrap()
        );
        let decimal_ln10 = context.ln10::<10>(None).with_precision(45).value();
        assert_eq!(
            decimal_ln10.repr.significand,
            IBig::from_str_radix("230258509299404568401799145468436420760110148", 10).unwrap()
        );

        let context = Context::<mode::Zero>::new(180);
        let binary_ln2 = context.ln2::<2>(None).with_precision(180).value();
        assert_eq!(
            binary_ln2.repr.significand,
            IBig::from_str_radix("1062244963371879310175186301324412638028404515790072203", 10)
                .unwrap()
        );
        let binary_ln10 = context.ln10::<2>(None).with_precision(180).value();
        assert_eq!(
            binary_ln10.repr.significand,
            IBig::from_str_radix("882175346869410758689845931257775553286341791676474847", 10)
                .unwrap()
        );
    }

    #[test]
    fn test_log2_domain() {
        let ctx = Context::<mode::HalfEven>::new(53);
        // log2(±0) = -inf (a value, not an error)
        let r = ctx.log2::<2>(&Repr::<2>::zero(), None).unwrap().value();
        assert!(r.repr.is_infinite());
        assert_eq!(r.repr.sign(), Sign::Negative);
        // log2(negative) is out of domain
        assert!(matches!(
            ctx.log2::<2>(&Repr::new((-1).into(), 0), None),
            Err(FpError::OutOfDomain)
        ));
        // an infinite input is rejected
        assert!(matches!(ctx.log2::<2>(&Repr::infinity(), None), Err(FpError::InfiniteInput)));
    }

    #[test]
    fn test_log2_exact_power_of_two() {
        // log2(2^k) = k exactly under every rounding mode. Regression for the directed-rounding
        // defect: rounding ln(x) and ln(2) each toward the mode and dividing once does not bound
        // the quotient, so previously log2(2^-159) under `Up` returned -159 + 1 ulp.
        let p = 53;
        for k in [0isize, 1, -1, 5, 159, -159, 1000, -1000] {
            let x = Repr::<2>::new(IBig::from(1), k); // 2^k
            let r_down = Context::<mode::Down>::new(p)
                .log2::<2>(&x, None)
                .unwrap()
                .value();
            let r_up = Context::<mode::Up>::new(p)
                .log2::<2>(&x, None)
                .unwrap()
                .value();
            let r_zero = Context::<mode::Zero>::new(p)
                .log2::<2>(&x, None)
                .unwrap()
                .value();
            let r_he = Context::<mode::HalfEven>::new(p)
                .log2::<2>(&x, None)
                .unwrap()
                .value();
            // Every directed mode produces the identical value — no mode-dependent ulp.
            assert_eq!(r_down.repr, r_he.repr, "Down != HalfEven for log2(2^{k})");
            assert_eq!(r_up.repr, r_he.repr, "Up != HalfEven for log2(2^{k})");
            assert_eq!(r_zero.repr, r_he.repr, "Zero != HalfEven for log2(2^{k})");
            // And that value is exactly k.
            assert_eq!(r_he.to_int().value(), IBig::from(k), "value for log2(2^{k})");
        }
    }

    #[test]
    fn test_log2_exact_power_of_two_decimal_base() {
        // In a non-power-of-two base the shortcut still fires when the exponent is zero: a
        // significand that is itself a power of two makes x = 2^m exactly.
        let p = 53;
        for (sig, want) in [(8i32, 3isize), (1024, 10), (2, 1), (32, 5)] {
            let x = Repr::<10>::new(IBig::from(sig), 0);
            let r_down = Context::<mode::Down>::new(p)
                .log2::<10>(&x, None)
                .unwrap()
                .value();
            let r_up = Context::<mode::Up>::new(p)
                .log2::<10>(&x, None)
                .unwrap()
                .value();
            let r_he = Context::<mode::HalfEven>::new(p)
                .log2::<10>(&x, None)
                .unwrap()
                .value();
            assert_eq!(r_down.repr, r_he.repr, "Down != HalfEven for log2({sig}) base 10");
            assert_eq!(r_up.repr, r_he.repr, "Up != HalfEven for log2({sig}) base 10");
            assert_eq!(r_he.to_int().value(), IBig::from(want), "value for log2({sig}) base 10");
        }
    }

    /// For a non-power-of-two significand `sig` (so `log2` is irrational and never lands on a
    /// rounding boundary), each directed result must equal a high-precision oracle rounded to the
    /// target precision under the same mode — the definition of correct rounding.
    fn check_log2_directed_matches_oracle<const B: Word>(sig: u32, p: usize) {
        let oracle_ctx = Context::<mode::HalfEven>::new(p + 40);
        let x = Repr::<B>::new(IBig::from(sig), 0);
        let oracle = oracle_ctx.log2::<B>(&x, None).unwrap().value();

        let want_down = Context::<mode::Down>::new(p)
            .repr_round_ref(&oracle.repr)
            .value();
        let want_up = Context::<mode::Up>::new(p)
            .repr_round_ref(&oracle.repr)
            .value();
        let want_he = Context::<mode::HalfEven>::new(p)
            .repr_round_ref(&oracle.repr)
            .value();

        let got_down = Context::<mode::Down>::new(p)
            .log2::<B>(&x, None)
            .unwrap()
            .value();
        let got_up = Context::<mode::Up>::new(p)
            .log2::<B>(&x, None)
            .unwrap()
            .value();
        let got_he = Context::<mode::HalfEven>::new(p)
            .log2::<B>(&x, None)
            .unwrap()
            .value();

        assert_eq!(got_down.repr, want_down, "log2({sig}) base {B} under Down");
        assert_eq!(got_up.repr, want_up, "log2({sig}) base {B} under Up");
        assert_eq!(got_he.repr, want_he, "log2({sig}) base {B} under HalfEven");
    }

    #[test]
    fn test_log2_directed_matches_oracle() {
        let p = 24;
        for sig in [3u32, 7, 10, 12345, 65537] {
            check_log2_directed_matches_oracle::<2>(sig, p);
        }
        // Exercise a non-power-of-two base through the Ziv interval path too.
        for sig in [3u32, 7, 10, 12345] {
            check_log2_directed_matches_oracle::<10>(sig, p);
        }
    }

    // log2 of a value whose result sits within ~1 work-ulp of a power of two must still round to
    // the correct neighbor under directed modes. log2(f64::MAX) ≈ 1024 − 2^-53/ln2 sits just below
    // 1024; under Down at p=53 the answer is 1024 − 2^-42 (the largest p=53 value ≤ it), but an
    // unsound radius previously let Ziv certify 1024 on the first attempt.
    #[test]
    fn test_log2_just_below_power_of_two_directed() {
        let x = FBig::<mode::HalfEven, 2>::try_from(f64::MAX).unwrap();
        // High-precision oracle, then re-rounded to the target precision under each mode.
        let oracle = Context::<mode::HalfEven>::new(200)
            .log2::<2>(x.repr(), None)
            .unwrap()
            .value();
        for p in [24usize, 40, 53, 64] {
            let want_down = Context::<mode::Down>::new(p)
                .repr_round_ref(&oracle.repr)
                .value();
            let want_up = Context::<mode::Up>::new(p)
                .repr_round_ref(&oracle.repr)
                .value();
            let got_down = Context::<mode::Down>::new(p)
                .log2::<2>(x.repr(), None)
                .unwrap()
                .value();
            let got_up = Context::<mode::Up>::new(p)
                .log2::<2>(x.repr(), None)
                .unwrap()
                .value();
            assert_eq!(got_down.repr(), &want_down, "p={p} Down");
            assert_eq!(got_up.repr(), &want_up, "p={p} Up");
            // Directed invariant: Up ≥ Down.
            assert!(got_up.repr() >= got_down.repr(), "p={p} Up < Down");
        }
    }

    /// Directed `ln` of `x ∈ [1, 2)` must match a high-precision oracle re-rounded under the same
    /// mode. This binade (s = 0) is where the radius under-estimated the error: `result` inherits
    /// `ln_base`'s over-delivered context, and for `x` just above 1 the scaling even classifies
    /// `s = −1`, so `2·sum + s·ln2` cancels and the error stays at `sum`'s magnitude while
    /// `result`'s collapses — both make `result.ulp()` the wrong scale for the radius.
    fn check_ln_directed_in_unit_binade(k: usize, p: usize) {
        // x = (2^k + 1) * 2^-k = 1 + 2^-k, exactly representable at precision p when k < p.
        let x = Repr::<2>::new(IBig::from(1i64 << k) + IBig::ONE, -(k as isize));
        let oracle = Context::<mode::HalfEven>::new(p + 60)
            .ln::<2>(&x, None)
            .unwrap()
            .value();
        let want_down = Context::<mode::Down>::new(p)
            .repr_round_ref(&oracle.repr)
            .value();
        let want_up = Context::<mode::Up>::new(p)
            .repr_round_ref(&oracle.repr)
            .value();
        let got_down = Context::<mode::Down>::new(p)
            .ln::<2>(&x, None)
            .unwrap()
            .value();
        let got_up = Context::<mode::Up>::new(p)
            .ln::<2>(&x, None)
            .unwrap()
            .value();
        assert_eq!(got_down.repr(), &want_down, "ln(1+2^-{k}) p={p} Down");
        assert_eq!(got_up.repr(), &want_up, "ln(1+2^-{k}) p={p} Up");
        assert!(got_up.repr() >= got_down.repr(), "ln(1+2^-{k}) p={p} Up < Down");
    }

    #[test]
    fn test_ln_directed_near_one() {
        // Sweep the near-1 binade at low precision, including the k close to p cases that
        // classify as s = −1 and cancel.
        for p in [24usize, 40, 53] {
            for k in 1..p.saturating_sub(1) {
                check_ln_directed_in_unit_binade(k, p);
            }
        }
    }

    /// Fixed inputs for the `log2` oracle differential: moderate magnitudes and the
    /// near-boundary regimes the legacy directed-interval implementation was specifically sized for.
    fn log2_diff_inputs() -> Vec<Repr<2>> {
        let mut v = Vec::new();
        for x in [0.5f64, 1.5, 2.0, 3.0, 10.0, 1000.0, 1e-6, 123.456, 2.5e-10] {
            v.push(FBig::<mode::HalfEven, 2>::try_from(x).unwrap().into_repr());
        }
        // Exact powers of two.
        for k in [-100isize, -50, -10, -1, 0, 1, 10, 50, 100] {
            v.push(Repr::new(IBig::ONE, k));
        }
        // Just below the largest f64 (log2 ≈ 1024, the directed-regime case in the old comment).
        v.push(
            FBig::<mode::HalfEven, 2>::try_from(f64::MAX)
                .unwrap()
                .into_repr(),
        );
        // The [1, 2) unit binade and its mirror below 1: 1 ± 2^-k and 2 − 2^-k exercise the
        // s = −1 cancellation (the second-classified-s-−1 case the doubling compensates).
        for k in 1usize..=60 {
            v.push(Repr::new(IBig::from(1u64 << k) + IBig::ONE, -(k as isize))); // 1 + 2^-k
            v.push(Repr::new(IBig::from((1u64 << k) - 1), -(k as isize))); // 1 − 2^-k
            v.push(Repr::new(IBig::from((1u64 << (k + 1)) - 1), -(k as isize)));
            // 2 − 2^-k
        }
        v
    }

    /// The Ball-based `log2` must round exactly like a high-precision oracle (the definition of
    /// correct rounding) across precisions, modes, and the near-boundary inputs.
    ///
    /// The legacy directed-interval implementation is *not* used as the oracle: it has its own
    /// residual 1-ulp bug under directed rounding for `log2(1 − 2^-k)` at p=50 (verified against
    /// an independent high-precision computation) — exactly the class of defect this pilot
    /// replaces.
    fn check_log2_differential<R: ErrorBounds>(p: usize, x: &Repr<2>, oracle: &Repr<2>) {
        let ctx = Context::<R>::new(p);
        let want = ctx.repr_round_ref(oracle).value();
        let got = ctx.log2_internal::<2>(x, None).unwrap().value();
        assert_eq!(got.repr, want, "p={p} {} x={x:?}", core::any::type_name::<R>(),);
    }

    /// Regression: `ln_compute`'s s<0 path (base < 1) must NOT inflate its error count with the
    /// working precision. An exactly-representable input scaled by a power of two is exact, so the
    /// radius must shrink monotonically as the work precision grows — otherwise the composed
    /// `pow_exp_log` chain's radius stays constant and the Ziv loop hangs (powf of a base < 1).
    #[test]
    fn ln_small_base_radius_shrinks_with_guard() {
        let ctx = Context::<mode::HalfEven>::new(50);
        // 0.2668 (base 10): s = floor(log2(0.2668)) = -2, the s < 0 path.
        let x = Repr::<10>::new(IBig::from(2668), -4);
        for guard in [4usize, 12, 40, 120] {
            let ball = ctx.ln_compute::<10>(&x, 50 + guard, false, None);
            // The regression: n must be O(series terms) (~10^5, bit_len < 30), NOT inflated to
            // ~B^50 ≈ 10^50 (bit_len ~166) by the s<0 reduction's spurious +1.
            assert!(
                ball.n.bit_len() < 30,
                "n = {} ({} bits) too large at guard={guard}: the s<0 reduction inflated it",
                ball.n,
                ball.n.bit_len()
            );
            // The radius in target (precision 50) ulps must fit a preimage so Ziv certifies on the
            // first attempt: n·B^(E−p_ball)·B^(50−E) ≤ 1.
            let radius_target = crate::ball::ceil_shift::<10>(
                ball.n.clone(),
                Ball::lead_exp(&ball.mid) - ball.mid.precision() as isize + 50,
            );
            assert!(
                radius_target <= IBig::ONE,
                "radius {radius_target} ulps at guard={guard} does not certify (n={})",
                ball.n
            );
        }
    }

    #[test]
    fn log2_ball_matches_oracle() {
        let inputs = log2_diff_inputs();
        // Moderate precisions: full input sweep, all five modes.
        for p in [20usize, 50, 100] {
            for x in &inputs {
                // The oracle is mode-independent: a high-precision HalfEven value re-rounded
                // under each target mode.
                let oracle = Context::<mode::HalfEven>::new(p + 60)
                    .log2::<2>(x, None)
                    .unwrap()
                    .value();
                check_log2_differential::<mode::HalfEven>(p, x, &oracle.repr);
                check_log2_differential::<mode::Down>(p, x, &oracle.repr);
                check_log2_differential::<mode::Up>(p, x, &oracle.repr);
                check_log2_differential::<mode::Zero>(p, x, &oracle.repr);
                check_log2_differential::<mode::Away>(p, x, &oracle.repr);
            }
        }
        // The arbitrary-precision regime: a reduced sweep (directed modes still exercised).
        for x in inputs.iter().step_by(9) {
            let oracle = Context::<mode::HalfEven>::new(560)
                .log2::<2>(x, None)
                .unwrap()
                .value();
            check_log2_differential::<mode::HalfEven>(500, x, &oracle.repr);
            check_log2_differential::<mode::Down>(500, x, &oracle.repr);
            check_log2_differential::<mode::Up>(500, x, &oracle.repr);
        }
    }

    #[test]
    fn ln_1p_ball_bounds_negative_arg() {
        // Regression: `ln_1p_ball`'s input-error adjust dropped the precision-difference term
        // (−p_arg+p_ln). For an arg with 1+arg ∈ (0, 1) (e.g. atanh(x<0) near the pole),
        // `ln_compute` doubles the work precision (the s<0 path), so `ln_ball` sits at 2p while
        // `arg` stays at p — the missing +p under-bounded the adjust by B^p and the radius no
        // longer covered the true value.
        use crate::fbig::FBig;
        use crate::repr::Context;
        type F = FBig<mode::HalfEven, 10>;
        let ctx = Context::<mode::HalfEven>::new(10);
        // arg mid = −0.9999 at precision 10 (ulp = 1e-10), n = 5 ⇒ true arg = −0.9999000005.
        let mid = F::from_parts(IBig::from(-9999000000i64), -10)
            .with_precision(10)
            .value();
        let arg = Ball::<10>::with_error(mid, IBig::from(5));
        let ln_ball = ctx.ln_1p_ball::<10>(&arg, None);
        // true ln(1+arg) = ln(1 − 0.9999000005) = ln(9.99995e-5), oracle at precision 60.
        let one_plus_true = F::from_parts(IBig::from(999995i64), -10)
            .with_precision(0)
            .value();
        let true_ln = one_plus_true
            .with_precision(60)
            .value()
            .ln()
            .with_precision(0)
            .value();
        let diff = (ln_ball.mid.clone().with_precision(0).value() - true_ln).abs();
        let bound = F::from(ln_ball.n.clone()) * ln_ball.mid.ulp().with_precision(0).value();
        assert!(
            diff <= bound,
            "ln_1p_ball: |mid − true| = {diff} > n·ulp = {bound} (n = {}, missing −p_arg+p_ln?)",
            ln_ball.n
        );
    }
}