dashu_float/

log.rs

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

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

/// Extra digits for the reconstruction constants (see `ln_compute`): a scale factor `|k|`
/// amplifies a constant's 8-ulp radius by `|k|`; `extra > log_B(|k|)` absorbs the `|k|`
/// amplification, leaving the scaled radius at ≤ 8 work-ulps (a constant fraction of the
/// result's ulp, sub-target once the Ziv guard is added; for the smallest bases it is a few
/// work-ulps, absorbed by the loop's initial guard). Computed from the factor's bit length,
/// integer-only: `log_B(|k|) < bit_len(|k|) / ilog2(B)`.
fn const_extra_digits<const B: Word>(mag_bit_len: usize) -> usize {
    mag_bit_len / (Word::BITS - 1 - B.leading_zeros()) as usize + 1
}

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,
                )
                .expect("ln(BASE) cannot range-error (the base fits the exponent range)")
                .to_value_radius::<R>(&Context::<R>::new(self.precision + guard))
                .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>,
    ) -> Result<Ball<B>, FpError> {
        let ctx = Context::<mode::HalfEven>::new(self.precision);
        match B {
            10 => {
                let logb = ctx.ln_base::<B>(reborrow_cache(&mut cache));
                let rad = ulps::<B>(&logb.repr, self.precision, 8);
                Ok(Ball::with_error(logb.into_repr(), rad))
            }
            i if i.is_power_of_two() => {
                let logb = ctx.ln_base::<B>(reborrow_cache(&mut cache));
                let rad = ulps::<B>(&logb.repr, self.precision, 8);
                Ok(Ball::with_error(logb.into_repr(), rad))
            }
            _ => {
                // 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.
    ///
    /// The magnitude is reduced onto the unit binade in two exact stages (a base-power
    /// exponent re-tag, then a small power-of-two division), so no power of an exponent *gap*
    /// is ever materialized — the reduction stays O(1) for any input exponent.
    ///
    /// The reconstruction `2·sum + s2·ln(2) + e_base·ln(B)` cancels for x < 1, so the series
    /// runs at double working precision — decided *before* the series, since a value-space
    /// radius needs no re-tagging when the precision changes (the mechanism the old ulp-count
    /// `rescale_precision` existed for is gone entirely).
    pub(crate) fn ln_compute<const B: Word>(
        &self,
        x: &Repr<B>,
        mut work_precision: usize,
        one_plus: bool,
        mut cache: Option<&mut ConstCache>,
    ) -> Result<Ball<B>, FpError> {
        // 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), work_precision);

        // 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();

        // The magnitude is reduced onto the unit binade in two exact stages, so that no power
        // of the exponent *gap* is ever materialized (a gap of ~10^9 digits would otherwise
        // allocate and divide integers of ~10^9 bits — `ln(1e1000000000)` never returned).
        //
        // 1. Split at a base-power boundary — a pure exponent re-tag (O(1), `Ball::shift`):
        //    with `d = digits_ub(sig) − 1`, the mid becomes `sig·B^(−d) ∈ [B^(−2), B)` and the
        //    removed factor is `B^(e+d)`. Unlike a single power-of-two reduction over the whole
        //    magnitude, this stage is exact integer exponent arithmetic; the f32 `log2_bounds`
        //    estimate loses ~hundreds of *bits* of accuracy once |log2 x| approaches 10^9, so
        //    deriving the whole reduction from it left the scaled value far outside [1, 2).
        //    Applied only for |e| ≥ 2: within |e| ≤ 1 the power-of-two reduction below shifts by
        //    at most `log2(sig) + 2·log2(B)` bits — proportional to the input — and (crucially)
        //    needs no `ln(B)` term in the reconstruction. That matters because `ln(B)` itself
        //    normalizes to `1·B^1`, so a base split of it would recurse onto its own
        //    reconstruction term (an infinite recursion on generic bases).
        // 2. Finish the reduction onto [1, 2) by a power of two: within this bounded range the
        //    f32 floor is off by at most one, and `2^|s2|` is a couple of words at most.
        //
        let (s2, e_base, x_scaled, below_one) = if no_scaling {
            (0, 0, x_ball, false)
        } else {
            let x_ball = if one_plus {
                x_ball.add(&Ball::exact_int(IBig::ONE, work_precision), work_precision)?
            } else {
                x_ball
            };
            // A zero midpoint means the `ln_1p` argument rounded onto the pole at `−1` and the sum
            // cancelled exactly (`1 + x` is within `rad` of zero, so it may as well be negative —
            // no logarithm exists). `log2_bounds` saturates at `−∞` there, which the reduction
            // below cannot use. The whole-line radius asks the Ziv loop for a higher guard, at
            // which the argument is representable and the sum is a nonzero difference again.
            if x_ball.mid.significand().is_zero() {
                return Ok(Ball::with_error(Repr::zero(), Mag::INFINITY));
            }
            // Does the (rounded) input sit below one? The exact cancellation condition for
            // the reconstruction below — one cheap comparison, immune to how the two-stage
            // reduction happens to split the exponent (master's single-stage `s < 0` read
            // the same predicate off its split).
            let below_one = x_ball.mid.cmp(&Repr::<B>::one()) == Ordering::Less;

            let e = x_ball.mid.exponent();
            let (e_base, x1) = if e >= 2 || e <= -2 {
                let d = x_ball.mid.digits_ub() as isize - 1;
                let e_base = e.saturating_add(d);
                (e_base, x_ball.shift(-e_base)) // value sig·B^(−d) ∈ [B^(−2), B)
            } else {
                (0, x_ball)
            };

            let log2 = x1.mid.log2_bounds().0;
            // The f32 bound is `next_down`-ed, so a midpoint that is *exactly* a power of two
            // (common: the base split above lands `B^e` inputs on `x1 = 1.0`) floors one step
            // too low — `s2 = −1` would scale the ball up to `2` for nothing. An exact
            // comparison settles the power-of-two boundary.
            let mut s2 = log2 as isize - (log2 < 0.) as isize; // floor(log2(x1))
            if s2 < 0 && x1.mid.cmp(&Repr::<B>::one()) != Ordering::Less {
                s2 = 0;
            }
            let x_scaled = if s2 > 0 {
                // Divide by the exact power 2^s2 as a ball division: an exact divisor shrinks
                // the radius by its value (`rad_b = 0` in the division rule) — no special case.
                let k = Ball::exact(Repr::new(IBig::ONE << s2 as usize, 0));
                x1.div(&k, work_precision)?
            } else if s2 < 0 {
                // Scaling by 2^|s2| is exact (finite decimal × power of two); an exact operand
                // keeps `rad = 0` through the multiplication's conditional ε.
                x1.scale_int(&(IBig::ONE << (-s2) as usize), work_precision)?
            } else {
                x1
            };
            debug_assert!(x_scaled.mid.cmp(&Repr::<B>::one()) != Ordering::Less);
            (s2, e_base, x_scaled, below_one)
        };

        // The reconstruction `2·sum + s2·ln(2) + e_base·ln(B)` *cancels* when the exact terms
        // of the scale factors dominate the result — precisely when the input is below one
        // (the series sum is then small against the constant terms), so the series runs at
        // double precision to keep the pre-cancellation sum accurate. A split with `s2 < 0`
        // but `x ≥ 1` (e.g. the power-two stage landing `x1 = 0.999…` below one) does not
        // cancel and stays at single precision.
        if below_one || x_scaled.mid.sign() == Sign::Negative {
            work_precision += self.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(IBig::from(2), work_precision);
            let den = x_scaled.add(&two, work_precision)?;
            x_scaled.div(&den, work_precision)?
        } else {
            let one = Ball::exact_int(IBig::ONE, work_precision);
            let num = x_scaled.sub(&one, work_precision)?;
            let den = x_scaled.add(&one, work_precision)?;
            num.div(&den, work_precision)?
        };
        let z2 = z.mul(&z, work_precision)?;
        let mut pow = z.clone();
        let mut sum = z;
        let mut k: usize = 3;
        loop {
            pow = pow.mul(&z2, work_precision)?;

            let increase = pow.div_int(k, work_precision)?;
            if increase.mid_le_ulp_lb(&sum, work_precision) {
                break;
            }

            sum = sum.add(&increase, work_precision)?;
            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. A zero sum means z (hence every term) is
        // exactly zero — the tail is 0.
        if !sum.mid.significand().is_zero() {
            sum.add_error(ulps::<B>(&sum.mid, work_precision, B as usize));
        }

        // Compose the logarithm of the original number: ln(x) = 2·sum + s2·ln(2) + e_base·ln(B).
        // The constants carry a handful of ulps of error (8 is a conservative sound bound for
        // every code path, cached and uncached), and the exact integer scale factors amplify
        // that error by |s2|/|e_base| while the result's ulp grows only by B^log_B(|scale|) —
        // so the constants are evaluated with enough extra digits to keep the amplified radius
        // below one work-ulp divided by the constant's own 8-ulp slack — i.e. `8/B` work-ulps,
        // sub-ulp for the base-10 origin of the construction and a few work-ulps on the
        // smallest bases (sound either way; Ziv certifies, and the initial guard absorbs it).
        // (The same construction as `exp_compute`'s `extra`.) Without them, ln of a large base
        // power sits at a sizable fraction of the target half-ulp and pays a Ziv retry
        // whenever the work value also lands near a rounding boundary.
        let sum2 = sum.scale_int(&IBig::from(2), work_precision)?;
        if no_scaling {
            Ok(sum2)
        } else if B.is_power_of_two() {
            // ln(B) = log2(B)·ln(2) exactly, so both scale factors fold into one ln(2) term.
            let k = IBig::from(s2) + IBig::from(e_base) * IBig::from(B.trailing_zeros());
            if k.is_zero() {
                return Ok(sum2);
            }
            let const_prec = work_precision + const_extra_digits::<B>(k.bit_len());
            let ln2 =
                Context::<mode::HalfEven>::new(const_prec).ln2::<B>(reborrow_cache(&mut cache));
            let rad = ulps::<B>(&ln2.repr, const_prec, 8);
            let ln2_ball = Ball::with_error(ln2.into_repr(), rad);
            sum2.add(&ln2_ball.scale_int(&k, const_prec)?, const_prec)
        } else {
            let scale = Ord::max(s2.unsigned_abs(), e_base.unsigned_abs());
            let const_prec = work_precision + const_extra_digits::<B>(scale.bit_len());
            let mut result = sum2;
            if s2 != 0 {
                let ln2 =
                    Context::<mode::HalfEven>::new(const_prec).ln2::<B>(reborrow_cache(&mut cache));
                let rad = ulps::<B>(&ln2.repr, const_prec, 8);
                let ln2_ball = Ball::with_error(ln2.into_repr(), rad);
                result =
                    result.add(&ln2_ball.scale_int(&IBig::from(s2), const_prec)?, const_prec)?;
            }
            if e_base != 0 {
                let ln_base_ball = Context::<mode::HalfEven>::new(const_prec)
                    .ln_base_ball::<B>(reborrow_cache(&mut cache))?;
                result = result
                    .add(&ln_base_ball.scale_int(&IBig::from(e_base), const_prec)?, const_prec)?;
            }
            Ok(result)
        }
    }

    /// `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.rad` then contributes
    /// `|θ| / LB(1 + arg_lo)` to the log — the derivative bound at the ball's lower endpoint,
    /// with the endpoint built by exact `Repr` arithmetic (the log1p rule as validated in
    /// dashu-ball, tighter than the old `×2` ulp-domain formula). An endpoint crossing zero
    /// means the ball straddles the pole: the whole-line radius forces a Ziv retry at higher
    /// precision.
    pub(crate) fn ln_1p_ball<const B: Word>(
        &self,
        arg: &Ball<B>,
        mut cache: Option<&mut ConstCache>,
    ) -> Result<Ball<B>, FpError> {
        // The ball must stay clear of the pole at `1 + arg = 0`: there is no logarithm to compute
        // (a zero sum has no `log2_bounds`, and a negative one is out of domain), and the
        // derivative fold below is unbounded. The whole-line radius asks the Ziv loop for a higher
        // guard, where the caller's argument stops rounding onto the pole — `atanh`'s `2x/(1−x)`
        // reaching `±∞` at the first attempt is exactly this case. A whole-line *radius* (another
        // operator's retry signal, e.g. a root that rounded onto zero) is the same situation.
        if !arg.mid.is_finite() || arg.rad.is_infinite() {
            return Ok(Ball::with_error(Repr::zero(), Mag::INFINITY));
        }
        let lo = &(Repr::<B>::one() + &arg.mid) - &arg.rad.to_repr::<B>();
        let denom = if lo.sign() == Sign::Negative {
            Mag::ZERO
        } else {
            Mag::from_repr_lower(&lo)
        };
        if denom.is_zero() {
            return Ok(Ball::with_error(Repr::zero(), Mag::INFINITY));
        }
        let mut ln_ball =
            self.ln_compute::<B>(&arg.mid, self.precision, true, reborrow_cache(&mut cache))?;
        ln_ball.add_error(arg.rad.div(&denom));
        Ok(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. `-1` is a pole,
        // so the boundary is compared *exactly* (no precision argument): `-1 + 10⁻¹⁵⁵` at 100
        // digits rounds onto `-1` but has the finite value `ln(10⁻¹⁵⁵)`.
        if x.sign() == Sign::Negative && !x.significand.is_zero() {
            match repr_cmp_same_base::<B, true>(x, &Repr::<B>::one(), None) {
                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>(&Context::<R>::new(self.precision + guard)))
        })
    }

    /// 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, work_precision)?
                .to_value_radius::<R>(&Context::<R>::new(work_precision)))
        })
    }

    /// 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, work_precision)?
                .to_value_radius::<R>(&Context::<R>::new(work_precision)))
        })
    }
}

/// 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::ball::ulp_mag;
    use crate::round::mode;
    use alloc::vec::Vec;
    use dashu_base::{Abs, AbsOrd};

    #[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 radius with a
    /// spurious rounding term. An exactly-representable input scaled by a power of two stays
    /// exact through the reduction, so the radius must stay small enough that the Ziv loop
    /// certifies at the target precision (the old ulp-count code inflated `n` by `B^precision`
    /// here and hung `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).unwrap();
            // The radius must certify at the target precision: rad ≤ 1 ulp of the midpoint at
            // precision 50 (i.e. the error interval fits a rounding preimage on this attempt).
            assert!(
                ball.rad <= ulp_mag::<10>(&ball.mid, 50),
                "radius {:?} at guard={guard} does not certify (the s<0 reduction inflated it)",
                ball.rad
            );
        }
    }

    #[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 the radius must still
        // cover the true value after the derivative fold.
        type F = FBig<mode::HalfEven, 10>;
        let ctx = Context::<mode::HalfEven>::new(10);
        // arg mid = −0.9999 at precision 10 (ulp = 1e-10), rad = 5 ulps ⇒ true arg = −0.9999000005.
        let mid = F::from_parts(IBig::from(-9999000000i64), -10)
            .with_precision(10)
            .value();
        let rad = ulps::<10>(&mid.repr, 10, 5);
        let arg = Ball::<10>::with_error(mid.into_repr(), rad);
        let ln_ball = ctx.ln_1p_ball::<10>(&arg, None).unwrap();
        // 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 mid_f = F::new(ln_ball.mid.clone(), Context::<mode::HalfEven>::new(0));
        let diff = (mid_f - true_ln).abs();
        let bound = F::new(ln_ball.rad.to_repr::<10>(), Context::<mode::HalfEven>::new(0));
        assert!(
            diff.abs_cmp(&bound).is_le(),
            "ln_1p_ball: |mid − true| = {diff} > rad = {bound}"
        );
    }

    /// Regression (issue #103): `ln` of a huge-exponent argument must reduce as
    /// `ln(m) + k·ln(B)` — the reduction may never materialize a power of the exponent *gap*
    /// (`ln(1e1000000000)` previously ballooned to gigabytes of integer and never returned).
    #[test]
    fn ln_huge_exponent_arguments() {
        use alloc::string::ToString;
        let ctx = Context::<mode::HalfEven>::new(34);
        // ln(10^1000000000) = 10^9·ln(10), correctly rounded to 34 digits
        let x = Repr::<10>::new(IBig::ONE, 1_000_000_000);
        let r = ctx.ln::<10>(&x, None).unwrap().value();
        assert_eq!(r.to_string(), "2302585092.994045684017991454684364");
        // ln(10^-1000000000) = −ln(10^1000000000) — the s<0 path (double work precision)
        let x = Repr::<10>::new(IBig::ONE, -1_000_000_000);
        let r = ctx.ln::<10>(&x, None).unwrap().value();
        assert_eq!(r.to_string(), "-2302585092.994045684017991454684364");
        // a mid-size exponent on the same path (the previously-working regime)
        let x = Repr::<10>::new(IBig::ONE, 1_000_000);
        let r = ctx.ln::<10>(&x, None).unwrap().value();
        assert_eq!(r.to_string(), "2302585.092994045684017991454684364");
    }

    /// Huge-exponent arguments under directed modes must match a high-precision oracle
    /// re-rounded under the same mode — exercising both reduction stages (the base-exponent
    /// re-tag and the bounded power-of-two finish) on bases 2 and 10, with exact and inexact
    /// significands.
    #[test]
    fn ln_directed_matches_oracle_huge_exponents() {
        fn check<const B: Word, R: ErrorBounds>(p: usize, x: &Repr<B>) {
            let oracle = Context::<mode::HalfEven>::new(p + 60)
                .ln::<B>(x, None)
                .unwrap()
                .value();
            let want = Context::<R>::new(p).repr_round_ref(&oracle.repr).value();
            let got = Context::<R>::new(p).ln::<B>(x, None).unwrap().value();
            assert_eq!(got.repr, want, "p={p} {} x={x:?} base={B}", core::any::type_name::<R>());
        }
        let inputs10 = [
            Repr::<10>::new(IBig::ONE, 1_000_000_000),
            Repr::<10>::new(IBig::from(3), 999_999_998), // 3e999999998
            Repr::<10>::new(IBig::ONE, -1_000_000_000),
            Repr::<10>::new(IBig::from(7), -1_000_000_000), // 7e-1000000000
            Repr::<10>::new(IBig::from(12345), -9_999_995), // 1.2345e-9999991
        ];
        for p in [20usize, 34, 50, 100] {
            for x in &inputs10 {
                check::<10, mode::HalfEven>(p, x);
                check::<10, mode::Down>(p, x);
                check::<10, mode::Up>(p, x);
                check::<10, mode::Zero>(p, x);
            }
        }
        let inputs2 = [
            Repr::<2>::new(IBig::ONE, 1 << 30),
            Repr::<2>::new(IBig::from(3), -(1 << 30)),
            Repr::<2>::new(IBig::ONE, -(1 << 20)),
        ];
        for p in [20usize, 50, 100, 500] {
            for x in &inputs2 {
                check::<2, mode::HalfEven>(p, x);
                check::<2, mode::Down>(p, x);
                check::<2, mode::Up>(p, x);
                check::<2, mode::Zero>(p, x);
            }
        }
    }

    // The reconstruction constants' 8-ulp radii are amplified by the base-power scale factor
    // (`100·rad(ln 10)` here), which used to sit at a third of the target half-ulp at 6 digits —
    // one Ziv retry whenever the work value also landed near the rounding boundary. The constants
    // now carry enough extra digits to keep the amplified radius sub-ulp (see `ln_compute`).
    #[cfg(feature = "std")]
    #[test]
    fn ln_large_base_power_certifies_first_attempt() {
        let ctx = Context::<mode::HalfEven>::new(6);
        crate::ziv_retries_reset();
        let _ = ctx.ln::<10>(&Repr::<10>::new(1.into(), 100), None).unwrap();
        assert_eq!(
            crate::ziv_retries(),
            0,
            "DBig ln(1e100) @6 should certify on the first attempt"
        );
    }

    /// `ln_1p` compares its `−1` boundary *exactly*: `−1 + 10⁻ⁿ` is the finite `ln(10⁻ⁿ)`, not `−∞`.
    ///
    /// Regression: `ln_compute` was handed the `1 + x` cancellation with a *zero* midpoint, where
    /// `log2_bounds` saturates to `−∞` — the reduction panicked on the cast (`attempt to subtract
    /// with overflow`, wrapping to a nonsense exponent in release). The zero now reports the
    /// whole-line radius, so the loop retries at a guard where the argument is representable.
    #[test]
    fn test_ln_1p_endpoint_compare_exactly() {
        for &p in &[20usize, 50, 100, 500] {
            let n = 3 * p / 2 + 5;
            // −1 + 10^-n = −(10^n − 1)·10^-n
            let x = Repr::<10>::new(-(IBig::from(10).pow(n) - IBig::from(1)), -(n as isize));
            let ctx = Context::<mode::HalfEven>::new(p);
            let v = ctx.ln_1p::<10>(&x, None).unwrap().value();
            assert!(v.repr().is_finite(), "ln_1p(-1 + 10^-{n}) at p={p} is finite, not -inf");
            let hi = Context::<mode::HalfEven>::new(p + 60)
                .ln_1p::<10>(&x, None)
                .unwrap()
                .value();
            assert_eq!(v, hi.with_precision(p).value(), "ln_1p(-1 + 10^-{n}) at p={p}");
            // exactly −1 is still the pole
            let at_pole = ctx
                .ln_1p::<10>(&Repr::<10>::neg_one(), None)
                .unwrap()
                .value();
            assert!(at_pole.repr().is_infinite() && at_pole.repr().sign() == Sign::Negative);
        }
    }
}