dashu_float/

log.rs

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

use crate::{
    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))
    }
}

// `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,
                )
                .0
            }
        }
    }

    /// 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 `(value, error_radius)`.
    ///
    /// 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. The radius is a provable upper bound on `|value − true|`, derived from
    /// the term count (every series step is correctly rounded; the truncated tail is `< 1 ulp` by
    /// the break test).
    pub(crate) fn ln_compute<const B: Word>(
        &self,
        x: &Repr<B>,
        mut work_precision: usize,
        one_plus: bool,
        mut cache: Option<&mut ConstCache>,
    ) -> (FBig<R, B>, FBig<R, B>) {
        // log(x) = log(x·B⁻ˢ) + s·log(B), with s = floor(log_B(x)) so x·B⁻ˢ ∈ [1, B).
        let context = Context::<R>::new(work_precision);
        let x = FBig::new(context.repr_round_ref(x).value(), 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.log2_est() < -B.log2_est();

        let (s, mut x_scaled) = if no_scaling {
            (0, x)
        } else {
            let x = if one_plus { x + FBig::ONE } else { x };

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

            let x_scaled = if B == 2 {
                x >> s
            } else if s > 0 {
                x / (IBig::ONE << s as usize)
            } else {
                x * (IBig::ONE << (-s) as usize)
            };
            debug_assert!(x_scaled >= FBig::<R, B>::ONE);
            (s, x_scaled)
        };

        if s < 0 || x_scaled.repr.sign() == Sign::Negative {
            // when s or x_scaled is negative, the final addition is actually a subtraction,
            // therefore we need to double the precision to get the correct result
            work_precision += self.precision;
            x_scaled.context.precision = work_precision;
        }
        let work_context = Context::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 d = &x_scaled + (FBig::ONE + FBig::ONE);
            x_scaled / d
        } else {
            (&x_scaled - FBig::ONE) / (x_scaled + FBig::ONE)
        };
        let z2 = z.sqr();
        let mut pow = z.clone();
        let mut sum = z;
        let mut terms: usize = 1; // the leading z term

        let mut k: usize = 3;
        loop {
            pow *= &z2;

            let increase = &pow / work_context.convert_int::<B>(k.into()).value();
            if increase.abs_cmp(&sum.ulp_lb()).is_le() {
                break;
            }

            sum += increase;
            k += 2;
            terms += 1;
        }

        // compose the logarithm of the original number
        let result: FBig<R, B> = if no_scaling {
            2 * sum.clone()
        } else {
            2 * sum.clone() + (s * work_context.ln2::<B>(reborrow_cache(&mut cache)))
        };

        // Provable error radius, expressed at the *work-precision* scale. Each series step rounds
        // once (< 1 ULP of the running sum) and the truncated tail is < 1 ULP by the break test, so
        // |sum − true| < (terms + 2)·ulp(sum); result = 2·sum + s·ln2 amplifies by ~2 and adds a few
        // reconstruction ULPs — `(4·terms + 12)·ulp` carries a generous margin.
        //
        // The ulp MUST be taken at the work precision, not `result.ulp()`: the `s·ln(B)` term calls
        // `ln_base`, which over-delivers (~`work_precision + guard` digits), and that inflated
        // *context* becomes `result`'s context — even when the `s·ln(B)` operand is zero (x ∈ [1, B),
        // so s = 0), `result` still carries `work + guard` precision while its value has only ~`work`
        // digits. `result.ulp()` is then ~`B^guard` finer than the work-precision ulp, and widening
        // by the *value's* extra digits (often ~0 here) does not compensate — only the *context*
        // inflation does. An unsound radius lets Ziv certify the wrong neighbor for a log result
        // within ~1 work-ulp of a rounding boundary (e.g. directed `ln`/`log2` in [1, B)). Widening
        // by `result.precision() − work_precision` is a no-op when `result` is already at work
        // precision.
        //
        // A second soundness hole: when `s < 0` the reconstruction `2·sum + s·ln(B)` *cancels* (for
        // x just above 1, `log2_bounds` can classify s as −1, so `2·sum ≈ |s|·ln(B)` and `result` is
        // their tiny difference). Subtraction preserves absolute error, so the series error stays at
        // `sum`'s magnitude (O(1)) while `result`'s magnitude collapses — `result.ulp()` then vastly
        // under-estimates the error. The radius must therefore cover the pre-cancellation (`sum`)
        // scale as well; taking the max is a no-op when the two magnitudes agree (s ≥ 0, no cancel).
        let inflation = result.precision().saturating_sub(work_precision);
        let work_ulp = result.ulp().with_precision(0).value() << inflation as isize;
        let sum_ulp = sum.ulp().with_precision(0).value();
        let radius = work_ulp.max(sum_ulp) * (4 * terms + 12);
        (result, radius)
    }
}

// `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),
            ))
        })
    }

    /// 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. Rounding ln(x) and ln(2)
        // separately and dividing once is only *near*-correct: under directed rounding, rounding
        // both operands toward the mode does not bound the quotient (enlarging a positive
        // denominator shrinks it). Instead each `ln_compute` reports a provable error radius, and
        // the two radii are carried through the division as an outward-rounded interval [lo, hi]
        // that is guaranteed to contain the true log2(x); the driver certifies once that interval
        // lies inside a single rounding bin.
        let initial_guard = self.base_guard_digits::<B>() + 4;
        self.ziv(initial_guard, |guard| {
            let work_precision = self.precision + guard;
            let (lx, ex) =
                self.ln_compute::<B>(x, work_precision, false, reborrow_cache(&mut cache));
            // ln(2) via the same near-correct primitive so it carries a provable radius too.
            let two = Repr::new(IBig::from(2), 0);
            let (l2, e2) =
                self.ln_compute::<B>(&two, work_precision, false, reborrow_cache(&mut cache));

            // True ln(x) ∈ [lx−ex, lx+ex] and true ln(2) ∈ [l2−e2, l2+e2] ⊆ (0, ∞). With a
            // positive denominator the quotient ln(x)/ln(2) is minimized by the low numerator
            // over the high denominator and maximized by the converse. Directing each endpoint's
            // rounding outward (lo down, hi up) keeps [lo, hi] a true containing interval.
            //
            // The interval bounds are computed at `work_precision + INTERVAL_GUARD` so the per-step
            // rounding (sub/add/div, each ≤ 1 ulp at the bound's magnitude) is well below one
            // *work* ulp. At work precision alone those roundings can shrink [lo, hi] enough that
            // it no longer contains the true quotient, leaving the radius unsound — which lets Ziv
            // certify the wrong neighbor for a result within ~1 work-ulp of a power of two (e.g.
            // log2(f64::MAX) ≈ 1024 under Down at p=53 returned 1024 instead of 1024 − 2^-42).
            // The extra digits make lo ≤ true ≤ hi rigorous, so `span = hi − lo` soundly bounds the
            // ln_compute spread and the `+ ulp_w` term covers only `value`'s own rounding.
            const INTERVAL_GUARD: usize = 16;
            let ip = work_precision + INTERVAL_GUARD;
            let down = Context::<mode::Down>::new(ip);
            let up = Context::<mode::Up>::new(ip);
            let nx_lo = down.sub(&lx.repr, &ex.repr).unwrap().value();
            let nx_hi = up.add(&lx.repr, &ex.repr).unwrap().value();
            let d_lo = down.sub(&l2.repr, &e2.repr).unwrap().value();
            let d_hi = up.add(&l2.repr, &e2.repr).unwrap().value();
            debug_assert!(
                d_lo.repr.sign() == Sign::Positive,
                "ln(2) lower bound must stay positive (guard digits keep e2 ≪ ln 2 ≈ 0.693)"
            );
            let lo = down.div(&nx_lo.repr, &d_hi.repr).unwrap().value();
            let hi = up.div(&nx_hi.repr, &d_lo.repr).unwrap().value();

            // Working-precision estimate; the driver re-rounds it to the target precision, so the
            // mode used here is immaterial to correctness.
            let value = Context::<R>::new(work_precision)
                .div(&lx.repr, &l2.repr)
                .unwrap()
                .value();

            // Radius: a provable bound on |value − true|. The true value lies in [lo, hi], and
            // `value` is within one working ulp of lx/l2 ∈ [lo, hi], so |value − true| ≤
            // (hi − lo) + ulp_w. Computed at unlimited precision so the bound arithmetic is exact
            // (no rounding that could under-report it), yet scaled by the working-precision span
            // and ulp so it shrinks as the guard grows and the loop converges. `lo`/`hi` were
            // rounded under Down/Up; their *values* are mode-independent, so rebuild them in the
            // target mode R via their reprs to keep the arithmetic single-mode.
            let unlim = Context::<R>::new(0);
            let span = FBig::new(hi.repr.clone(), unlim) - FBig::new(lo.repr.clone(), unlim);
            let ulp_w = value.ulp().with_precision(0).value();
            let radius = span + ulp_w;
            Ok((value, radius))
        })
    }
}

#[cfg(test)]
mod tests {
    use super::*;
    use crate::round::mode;

    #[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);
            }
        }
    }
}