Address opendp-num fuzz findings (DASHU-007/008/015/022/026)

AGENTS.md

@@ -50,6 +50,7 @@ Note: always `--exclude dashu-python` when running workspace-wide commands, sinc
 - When borrowing an algorithm idea from GMP (or any other library), do **not** reference its function names in our docstrings or comments. Describe the algorithm in our own terms and use our own function names (e.g. write `add_mul_dword_same_len_in_place`, never `addmul_2` / `mpn_addmul_2`). External function names must not appear anywhere in the repo.
 - Tests for a specific algorithm/kernel belong in the same source file as the implementation, as a `#[cfg(test)] mod tests` block at the bottom — not in a separate integration test file under `tests/`. Reserve `tests/` for cross-cutting or public-API tests.
 - When debugging or writing test assertions, use `{:?}` (or `{:#?}` for the verbose form with digit/bit counts) to inspect arbitrary precision values. The [`Debug`] format prints a compact head‥tail representation (most significant digits `..` least significant digits) instead of dumping the entire number, making it readable even for thousand-digit integers.
+- In-crate tests (`#[cfg(test)] mod tests` and `tests/`) must use **fixed, deterministic inputs** — never pseudo-random generators, fuzzing harnesses, or property-test styles that draw random data. Tests must fail or pass identically on every run. The only exception is code that tests the `rand` integration itself. For randomized / exploratory input, use the fuzz targets under `fuzz/` (the proper home for that work), not unit tests.
 
 ## Feature flags
 

float/CHANGELOG.md

@@ -6,6 +6,9 @@
 - `Repr::new_const`: a `const`-evaluable, normalized `Repr` constructor from a `DoubleWord`
   significand (the `const` counterpart of `Repr::new`). `FBig::from_parts_const` now delegates to
   it, and the complex literal macro uses it.
+- `FBig::log2` / `Context::log2` / `CachedFBig::log2`: base-2 logarithm, correctly rounded via
+  `ln(x)/ln(2)` evaluated at an elevated working precision. Previously only the f32-precision
+  `log2_bounds` magnitude estimate was available, so directed `log2` was wrong by many ULPs.
 
 ### Fix
 - `to_f64`/`to_f32` now round the source once, directly to the target's precision at its own
@@ -16,6 +19,12 @@
   high-precision inputs: the base-changing conversion now pre-shrinks the source significand before
   dividing, upholding `repr_div`'s dividend-width contract (as `Context::div` already does) instead
   of feeding an oversized dividend into the division.
+- `exp`, `exp_m1`, `sqrt`, `ln`, and `ln_1p` no longer panic on exact zero/one inputs that carry
+  unlimited precision (precision 0), such as `FBig::try_from(0.0)` and the `FBig::ONE`/`ZERO`
+  constants. The exact-result shortcuts now run before the limited-precision assertion.
+- The `round_fract` debug assertion no longer materializes `B^precision`, which for a sparse sticky
+  tail (where `precision` is the exponent gap — e.g. `exp_m1` of a large-magnitude input) could
+  exhaust memory in debug builds. It now checks the precondition with `log2` bounds.
 
 ## 0.5.0
 

float/src/exp.rs

@@ -109,6 +109,11 @@ impl<R: Round> Context<R> {
     /// Panics if the precision is unlimited and the exponent is negative. In this case, the exact
     /// result is likely to have infinite digits.
     pub fn powi<const B: Word>(&self, base: &Repr<B>, exp: IBig) -> FpResult<FBig<R, B>> {
+        // TODO: range handling has three known limitations at the exponent extremes:
+        // (1) the overflow guard below estimates the result magnitude with an f64 and misclassifies
+        // representable boundaries near 2^63; (2) genuine Overflow/Underflow is unwrapped mode-blindly
+        // to ±inf / signed zero; (3) the negative-exponent reciprocal path can panic. None affects
+        // ordinary inputs; fixing requires mode-aware range saturation.
         if base.is_infinite() {
             return Err(FpError::InfiniteInput);
         }
@@ -335,11 +340,13 @@ impl<R: Round> Context<R> {
         mut cache: Option<&mut ConstCache>,
     ) -> FpResult<FBig<R, B>> {
         assert_finite(x);
-        assert_limited_precision(self.precision);
         let input_sign = x.sign();
 
         if x.significand.is_zero() {
-            // exp(±0) = 1; exp_m1(±0) = ±0 (IEEE 754 §9.2.1 preserves the sign of zero)
+            // exp(±0) = 1; exp_m1(±0) = ±0 (IEEE 754 §9.2.1 preserves the sign of zero).
+            // These exact results need no rounding, so handle them before the
+            // limited-precision assertion: a precision-0 (unlimited) FBig such as the
+            // one produced by `try_from(0.0)` must still compute exp/exp_m1 exactly.
             return match minus_one {
                 false => Ok(Exact(FBig::ONE)),
                 true => {
@@ -353,6 +360,8 @@ impl<R: Round> Context<R> {
             };
         }
 
+        assert_limited_precision(self.precision);
+
         // A simple algorithm:
         // - let r = (x - s logB) / Bⁿ, where s = floor(x / logB), such that r < B⁻ⁿ.
         // - if the target precision is p digits, then there're only about p/m terms in Tyler series
@@ -399,6 +408,11 @@ impl<R: Round> Context<R> {
                 Err(_) => {
                     // |floor(x / ln B)| overflows isize — x is astronomically large, so the
                     // result is an infinity (x → +∞) or underflows to the limit (x → −∞).
+                    //
+                    // TODO: this branch discards the rounding mode. For a huge
+                    // *negative* x the true exp(x) is positive but below the exponent range, so
+                    // directed Up/Away should return the minimum positive value (and exp_m1 the
+                    // value just above -1) rather than +0 / Exact(-1).
                     return if input_sign == Sign::Positive {
                         Err(FpError::Overflow(Sign::Positive))
                     } else if minus_one {
@@ -473,6 +487,56 @@ mod tests {
         assert_eq!(m1, -FBig::<mode::HalfEven>::ONE);
     }
 
+    // A sharp OOM regression needs an exponent gap large enough that 2^gap exceeds any
+    // memory (gap ≳ 1e11), yet with floor(x/ln2) still fitting isize so the overflow
+    // branch is not taken. That window only exists where isize is 64-bit: on 32-bit,
+    // isize tops out at ~2.1e9 — below any OOM-inducing gap — so the overflow branch
+    // always intervenes first. The fix itself (log2_bounds in round_fract) is
+    // arch-independent; only this dedicated sharp test is 64-bit-only.
+    #[cfg(target_pointer_width = "64")]
+    #[test]
+    fn test_exp_m1_large_negative_no_oom() {
+        // Regression test: exp_m1 of a large-magnitude negative input subtracts
+        // 1 from an astronomically small exp(x), so the aligned subtraction hands
+        // round_fract a sparse sticky tail whose `precision` equals the exponent gap
+        // (~6.6e18 here). The debug assert in round_fract used to materialize B^precision
+        // (2^6.6e18 bits → certain OOM); it now uses log2 bounds and must not allocate.
+        //
+        // x = -2^62: the gap (~6.6e18) dwarfs memory, yet floor(x/ln2) ≈ -6.6e18 still
+        // fits 64-bit isize (≈9.2e18), so this does NOT enter the overflow branch.
+        let ctx = Context::<mode::Up>::new(2);
+        let x = Repr::new(-(IBig::from(1) << 62), 0);
+        let r = ctx.exp_m1::<2>(&x, None).unwrap().value();
+        // exp_m1(-2^62) = -1 + ε (ε > 0 vanishingly small); under Up it rounds above -1.
+        assert!(
+            r > -FBig::<mode::Up, 2>::ONE,
+            "exp_m1(-2^62) under Up should round above -1, got {r:?}"
+        );
+    }
+
+    #[test]
+    fn test_exact_results_on_unlimited_precision() {
+        // Regression test: values carrying precision 0 (unlimited) — produced by
+        // `try_from(0.0)` and the `FBig::ONE`/`ZERO` constants — must still compute
+        // their exact-result special cases instead of panicking in
+        // assert_limited_precision before reaching the shortcut.
+        type F = FBig<mode::HalfEven, 2>;
+
+        let zero = F::try_from(0.0_f64).unwrap();
+        assert_eq!(zero.exp(), F::ONE);
+        assert_eq!(zero.exp_m1(), F::ZERO);
+        assert_eq!(zero.sqrt(), F::ZERO);
+        assert_eq!(zero.ln_1p(), F::ZERO);
+
+        // -0.0 preserves its sign through exp_m1 and sqrt.
+        let neg_zero = F::try_from(-0.0_f64).unwrap();
+        assert!(neg_zero.exp_m1().repr().is_neg_zero());
+        assert!(neg_zero.sqrt().repr().is_neg_zero());
+
+        // FBig::ONE carries unlimited precision; ln(1) = 0 is exact.
+        assert_eq!(F::ONE.ln(), F::ZERO);
+    }
+
     #[test]
     fn test_powf_zero_base() {
         use crate::DBig;

float/src/fbig_cached.rs

@@ -539,6 +539,7 @@ mod tests {
         assert_eq!(x.clone().cos().into_fbig(), y.clone().cos());
         assert_eq!(x.clone().exp_m1().into_fbig(), y.clone().exp_m1());
         assert_eq!(x.clone().ln_1p().into_fbig(), y.clone().ln_1p());
+        assert_eq!(x.clone().log2().into_fbig(), y.clone().log2());
         assert_eq!(x.powf(&x.clone()).into_fbig(), y.clone().powf(&y));
     }
 

float/src/fbig_cached_ops.rs

@@ -666,6 +666,7 @@ macro_rules! forward_to_fbig {
 impl<R: Round, const B: Word> CachedFBig<R, B> {
     forward_to_context!(ln);
     forward_to_context!(ln_1p);
+    forward_to_context!(log2);
     forward_to_context!(exp);
     forward_to_context!(exp_m1);
 

float/src/log.rs

@@ -97,6 +97,28 @@ impl<R: Round, const B: Word> FBig<R, B> {
     pub fn ln_1p(&self) -> Self {
         self.context.unwrap_fp(self.context.ln_1p(&self.repr, None))
     }
+
+    /// Calculate the base-2 logarithm function (`log2(x)`) on the float number.
+    ///
+    /// This is evaluated as `ln(x) / ln(2)` at an elevated working precision and rounded
+    /// once, so it is correctly rounded to this number's precision (unlike
+    /// [`log2_bounds`][dashu_base::EstimatedLog2::log2_bounds], which only gives an
+    /// f32-precision magnitude estimate).
+    ///
+    /// # Examples
+    ///
+    /// ```
+    /// # use core::str::FromStr;
+    /// # use dashu_base::ParseError;
+    /// # use dashu_float::DBig;
+    /// let a = DBig::from_str("8.0")?;
+    /// 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))
+    }
 }
 
 impl<R: Round> Context<R> {
@@ -214,6 +236,48 @@ impl<R: Round> Context<R> {
         Ok(self.ln_internal(x, false, cache))
     }
 
+    /// Calculate the base-2 logarithm function (`log2(x)`) on the float number under this
+    /// context.
+    ///
+    /// This is evaluated as `ln(x) / ln(2)` at an elevated working precision and rounded
+    /// once, so it is correctly rounded to this context's precision (unlike
+    /// [`log2_bounds`][dashu_base::EstimatedLog2::log2_bounds], which only gives an
+    /// f32-precision magnitude estimate).
+    ///
+    /// # Examples
+    ///
+    /// ```
+    /// # use core::str::FromStr;
+    /// # use dashu_base::ParseError;
+    /// # use dashu_float::DBig;
+    /// use dashu_base::Approximation::*;
+    /// use dashu_float::{Context, round::mode::HalfEven};
+    ///
+    /// let context = Context::<HalfEven>::new(8);
+    /// let a = DBig::from_str("10")?;
+    /// // log2(10) ≈ 3.3219281
+    /// assert_eq!(context.log2(&a.repr(), None).unwrap().value(), DBig::from_str("3.3219281")?);
+    /// # 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);
+        }
+        Ok(self.log2_internal(x, cache))
+    }
+
     /// Calculate the natural logarithm function (`log(x+1)`) on the float number under this context.
     ///
     /// # Examples
@@ -257,8 +321,10 @@ impl<R: Round> Context<R> {
         mut cache: Option<&mut ConstCache>,
     ) -> Rounded<FBig<R, B>> {
         assert_finite(x);
-        assert_limited_precision(self.precision);
 
+        // 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 Exact(FBig::ZERO); // ln(1) = +0
         }
@@ -272,6 +338,8 @@ impl<R: Round> Context<R> {
             return Exact(zero);
         }
 
+        assert_limited_precision(self.precision);
+
         // A simple algorithm:
         // - let log(x) = log(x/2^s) + slog2 where s = floor(log2(x))
         // - such that x*2^s is close to but larger than 1 (and x*2^s < 2)
@@ -344,6 +412,32 @@ impl<R: Round> Context<R> {
         };
         result.with_precision(self.precision)
     }
+
+    /// `log2(x) = ln(x) / ln(2)`, correctly rounded to this context's precision.
+    ///
+    /// Both `ln(x)` and `ln(2)` are evaluated at a working precision above the target and
+    /// the quotient is rounded once. The result magnitude `|log2(x)| ≈ |floor(log2 x)|`
+    /// grows with the operand, which tracks the division's error amplification, so the
+    /// per-ulp error is ~`B^(target - work)`: a few guard digits certify the final round
+    /// across the whole magnitude range (no manual scaling needed).
+    fn log2_internal<const B: Word>(
+        &self,
+        x: &Repr<B>,
+        mut cache: Option<&mut ConstCache>,
+    ) -> Rounded<FBig<R, B>> {
+        assert_finite(x);
+        if x.is_one() {
+            return Exact(FBig::ZERO); // log2(1) = +0 (exact; also covers unlimited precision)
+        }
+
+        let guard = ceil_usize(self.precision.log2_est() / B.log2_est()) + 3;
+        let work = Context::<R>::new(self.precision + guard);
+        let ln_x = work
+            .ln_internal(x, false, reborrow_cache(&mut cache))
+            .value();
+        let ln2 = work.ln2::<B>(reborrow_cache(&mut cache));
+        (ln_x / ln2).with_precision(self.precision)
+    }
 }
 
 #[cfg(test)]
@@ -359,6 +453,85 @@ mod tests {
         assert_eq!(r.repr().sign(), Sign::Negative);
     }
 
+    #[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_eq!(ctx.log2::<2>(&Repr::new(IBig::from(-1), 0), None), Err(FpError::OutOfDomain));
+    }
+
+    #[test]
+    fn test_log2_powers_of_two() {
+        // log2(2^k) = k, including log2(1) = 0. The only magnitude estimate available
+        // (log2_bounds) is f32-precision, so deriving directed log2 from it gives a small
+        // *negative* value for log2(1) and is loose by many ULPs elsewhere; the
+        // correctly-rounded log2 returns the exact integer value for every power of two
+        // (log2(1) is tagged Exact via the shortcut; other powers are value-exact).
+        let ctx = Context::<mode::HalfEven>::new(53);
+        for k in [0usize, 1, 2, 3, 10, 50, 100, 200, 1000] {
+            let x = Repr::new(IBig::from(1) << k, 0); // 2^k
+            let r = ctx.log2::<2>(&x, None).unwrap();
+            let expected = ctx.convert_int::<2>(IBig::from(k)).value();
+            assert_eq!(r.value(), expected, "log2(2^{k}) = {k}");
+        }
+        // log2(1) = 0 is certified Exact via the special-case shortcut.
+        assert!(matches!(ctx.log2::<2>(&Repr::new(IBig::from(1), 0), None).unwrap(), Exact(_)));
+    }
+
+    #[test]
+    fn test_log2_correctly_rounded() {
+        // log2(x) must equal the correctly-rounded ln(x)/ln(2) oracle (computed at much
+        // higher precision then rounded down) across magnitude ranges — tiny, normal, and
+        // huge inputs, in both base 2 and base 10.
+        let p = 30;
+        let ctx = Context::<mode::HalfEven>::new(p);
+        let hi = Context::<mode::HalfEven>::new(200);
+
+        let cases_base2 = [
+            Repr::new(IBig::from(3), 0),
+            Repr::new(IBig::from(10), 0),
+            Repr::new(IBig::from(1), -24), // near the fuzzer's 2^-24 region
+            Repr::new(IBig::from(1) << 50, 0), // huge
+            Repr::new(IBig::from(1), -50), // tiny
+        ];
+        for x in cases_base2 {
+            let got = ctx.log2::<2>(&x, None).unwrap().value();
+            let want = (hi.ln::<2>(&x, None).unwrap().value() / hi.ln2::<2>(None))
+                .with_precision(p)
+                .value();
+            assert_eq!(got, want, "base-2 log2 mismatch for x={x:?}");
+        }
+
+        // base 10: log2(x) = ln(x)/ln(2) with a decimal significand
+        let x = Repr::new(IBig::from(123456), 0);
+        let got = ctx.log2::<10>(&x, None).unwrap().value();
+        let want = (hi.ln::<10>(&x, None).unwrap().value() / hi.ln2::<10>(None))
+            .with_precision(p)
+            .value();
+        assert_eq!(got, want, "base-10 log2 mismatch");
+    }
+
+    #[test]
+    fn test_log2_directed_bounds() {
+        // Down rounds toward -inf, Up toward +inf, so the correctly-rounded (HalfEven)
+        // value must lie within [Down, Up].
+        use crate::round::mode::{Down, Up};
+        let p = 20;
+        let mid = Context::<mode::HalfEven>::new(p);
+        for x in [Repr::new(IBig::from(10), 0), Repr::new(IBig::from(3), 13)] {
+            let down = Context::<Down>::new(p).log2::<2>(&x, None).unwrap().value();
+            let up = Context::<Up>::new(p).log2::<2>(&x, None).unwrap().value();
+            let half = mid.log2::<2>(&x, None).unwrap().value();
+            assert!(down <= half, "down <= half failed for x={x:?}: {down:?} vs {half:?}");
+            assert!(half <= up, "half <= up failed for x={x:?}: {half:?} vs {up:?}");
+            assert!(down <= up, "down <= up failed for x={x:?}");
+        }
+    }
+
     #[test]
     fn test_iacoth() {
         let context = Context::<mode::Zero>::new(10);

float/src/root.rs

@@ -90,11 +90,13 @@ impl<R: Round> Context<R> {
         if x.is_infinite() {
             return Err(FpError::InfiniteInput);
         }
-        assert_limited_precision(self.precision);
         if x.significand.is_zero() {
-            // sqrt(+0) = +0, sqrt(-0) = -0 (preserve the sign of zero)
+            // sqrt(+0) = +0, sqrt(-0) = -0 (preserve the sign of zero). Exact, so handle
+            // it before the limited-precision assertion: a precision-0 (unlimited) value
+            // such as the one from `try_from(0.0)` must still compute sqrt(0) exactly.
             return Ok(Approximation::Exact(FBig::new(x.clone(), *self)));
         }
+        assert_limited_precision(self.precision);
         if x.sign() == Sign::Negative {
             return Err(FpError::OutOfDomain);
         }