diff --git a/AGENTS.md b/AGENTS.md index d0133a64..d9439210 100644 --- a/AGENTS.md +++ b/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 diff --git a/float/CHANGELOG.md b/float/CHANGELOG.md index 4d450cc2..29b41346 100644 --- a/float/CHANGELOG.md +++ b/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 diff --git a/float/src/exp.rs b/float/src/exp.rs index 9d8152f3..54948119 100644 --- a/float/src/exp.rs +++ b/float/src/exp.rs @@ -109,6 +109,11 @@ impl Context { /// 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(&self, base: &Repr, exp: IBig) -> FpResult> { + // 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 Context { mut cache: Option<&mut ConstCache>, ) -> FpResult> { 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 Context { }; } + 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 Context { 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::::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::::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::::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; + + 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; diff --git a/float/src/fbig_cached.rs b/float/src/fbig_cached.rs index 1ea2a645..87bebc33 100644 --- a/float/src/fbig_cached.rs +++ b/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)); } diff --git a/float/src/fbig_cached_ops.rs b/float/src/fbig_cached_ops.rs index c3f44e4c..3ac82869 100644 --- a/float/src/fbig_cached_ops.rs +++ b/float/src/fbig_cached_ops.rs @@ -666,6 +666,7 @@ macro_rules! forward_to_fbig { impl CachedFBig { forward_to_context!(ln); forward_to_context!(ln_1p); + forward_to_context!(log2); forward_to_context!(exp); forward_to_context!(exp_m1); diff --git a/float/src/log.rs b/float/src/log.rs index 02772914..325dbddf 100644 --- a/float/src/log.rs +++ b/float/src/log.rs @@ -97,6 +97,28 @@ impl FBig { 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 Context { @@ -214,6 +236,48 @@ impl Context { 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::::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( + &self, + x: &Repr, + cache: Option<&mut ConstCache>, + ) -> FpResult> { + 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 Context { mut cache: Option<&mut ConstCache>, ) -> Rounded> { 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 Context { 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 Context { }; 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( + &self, + x: &Repr, + mut cache: Option<&mut ConstCache>, + ) -> Rounded> { + 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::::new(self.precision + guard); + let ln_x = work + .ln_internal(x, false, reborrow_cache(&mut cache)) + .value(); + let ln2 = work.ln2::(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::::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::::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::::new(p); + let hi = Context::::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::::new(p); + for x in [Repr::new(IBig::from(10), 0), Repr::new(IBig::from(3), 13)] { + let down = Context::::new(p).log2::<2>(&x, None).unwrap().value(); + let up = Context::::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::::new(10); diff --git a/float/src/root.rs b/float/src/root.rs index 5e806285..465517f6 100644 --- a/float/src/root.rs +++ b/float/src/root.rs @@ -90,11 +90,13 @@ impl Context { 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); } diff --git a/float/src/round.rs b/float/src/round.rs index 6aaa61ff..265a2833 100644 --- a/float/src/round.rs +++ b/float/src/round.rs @@ -2,7 +2,7 @@ use core::cmp::Ordering; use core::ops::{Add, AddAssign}; -use dashu_base::{AbsOrd, Approximation, BitTest, EstimatedLog2, Sign, Signed, UnsignedAbs}; +use dashu_base::{AbsOrd, Approximation, BitTest, EstimatedLog2, Sign, Signed}; use dashu_int::{IBig, UBig, Word}; use crate::FBig; @@ -104,8 +104,18 @@ pub trait Round: Copy { /// assuming |fract| / X^precision < 1. Return the adjustment. #[inline] fn round_fract(integer: &IBig, fract: IBig, precision: usize) -> Rounding { - // this assertion is costly, so only check in debug mode - debug_assert!(fract.clone().unsigned_abs() < UBig::from_word(B).pow(precision)); + // this assertion is costly, so only check in debug mode. + // Verify the precondition |fract| < B^precision *without* materializing B^precision: + // for a sparse sticky tail produced by aligned subtraction, `precision` is the + // exponent gap and can be astronomically large, so building B^precision (the old + // check) exhausts memory in debug builds. Use log2 bounds instead, and only flag a + // *proven* violation so a valid input never trips the assertion. + debug_assert!({ + let (lb, _ub) = fract.log2_bounds(); // bounds on log2|fract| (−inf when fract == 0) + let (_b_lb, b_ub) = B.log2_bounds(); + // certain violation (|fract| >= B^precision) only when lb >= precision * b_ub + lb < b_ub * precision as f32 + }); if fract.is_zero() { return Rounding::NoOp; diff --git a/integer/CHANGELOG.md b/integer/CHANGELOG.md index a478bd03..63053385 100644 --- a/integer/CHANGELOG.md +++ b/integer/CHANGELOG.md @@ -1,5 +1,18 @@ # Changelog +## Unreleased + +### Fix +- `to_f64` of a magnitude at `DoubleWord::MAX` (e.g. `IBig::from(-(2^128 - 1)).to_f64()` on 64-bit + targets) reported the inexact conversion as `Approximation::Exact`. The exactness test used a + saturating `f as DoubleWord` round-trip, which clamps back to `DoubleWord::MAX` when the value + rounds up to `2^BITS`; that saturation case is now detected before the comparison. The same guard + is applied to the 16/32-bit `RefLarge` → f64 fast path. +- GCD of two large, similarly-sized integers no longer aborts with "internal error: not enough memory + allocated" when a later euclidean step takes the Burnikel-Ziegler path. The scratchpad was reserved + once from the *initial* operand lengths, but each step's division dispatches on the *current* + lengths; it is now sized for the worst case reachable during reduction (`mul::memory_requirement_up_to(rhs_len, rhs_len/2)`). + ## 0.5.0 ### Add diff --git a/integer/src/convert.rs b/integer/src/convert.rs index 8ea85175..fb3383a8 100644 --- a/integer/src/convert.rs +++ b/integer/src/convert.rs @@ -1085,8 +1085,16 @@ mod repr { if self.bit_len() <= 64 { let v: u64 = self.try_to_unsigned().unwrap(); let f = v as f64; - let back = f as u64; - match back.cmp(&v) { + // Same saturation guard as `to_f64_small`: `f as u64` clamps to + // `u64::MAX` when `v` rounds up to 2^64, which would otherwise + // report `Equal` (Exact) for `v == u64::MAX`. + const TWO_POW_64: f64 = (1u64 << 63) as f64 * 2.0; + let cmp = if f >= TWO_POW_64 { + Ordering::Greater + } else { + (f as u64).cmp(&v) + }; + match cmp { Ordering::Greater => Inexact(f, Sign::Positive), Ordering::Equal => Exact(f), Ordering::Less => Inexact(f, Sign::Negative), @@ -1130,9 +1138,19 @@ mod repr { fn to_f64_small(dword: DoubleWord) -> Approximation { const_assert!((DoubleWord::MAX as f64) < f64::MAX); let f = dword as f64; - let back = f as DoubleWord; - match back.partial_cmp(&dword).unwrap() { + // `f as DoubleWord` is a saturating float→int cast: when `dword` rounds up to + // 2^BITS (the f64 just above `DoubleWord::MAX`), the cast clamps back to + // `DoubleWord::MAX == dword` and the round-trip comparison reports `Equal`, + // wrongly tagging the conversion `Exact`. `dword as f64` is at most 2^BITS, so + // reaching 2^BITS marks exactly that saturation case — the value rounded up. + const TWO_POW_BITS: f64 = (1u128 << (DoubleWord::BITS - 1)) as f64 * 2.0; + let cmp = if f >= TWO_POW_BITS { + Ordering::Greater + } else { + (f as DoubleWord).partial_cmp(&dword).unwrap() + }; + match cmp { Ordering::Greater => Inexact(f, Sign::Positive), Ordering::Equal => Exact(f), Ordering::Less => Inexact(f, Sign::Negative), diff --git a/integer/src/gcd/lehmer.rs b/integer/src/gcd/lehmer.rs index 00a4bb4e..40db187a 100644 --- a/integer/src/gcd/lehmer.rs +++ b/integer/src/gcd/lehmer.rs @@ -225,9 +225,21 @@ pub(crate) fn lehmer_step(x: &mut [Word], y: &mut [Word], a: Word, b: Word, c: W /// Temporary memory required for gcd. #[inline] pub fn memory_requirement_up_to(lhs_len: usize, rhs_len: usize) -> Layout { - // Required memory: - // - temporary space for the division in the euclidean step - div::memory_requirement_exact(lhs_len, rhs_len) + // The scratchpad is reserved ONCE from the initial operand lengths, but each + // euclidean step's division dispatches on the *current* lengths. As Lehmer shrinks + // y while x stays large, a later step can divide a still-large x by a much smaller y + // and take the Burnikel-Ziegler path — even when the first (similarly-sized) + // division was "simple" and needed no scratch. `div::memory_requirement_exact` only + // bounds a *single* division with the given lengths, so it under-reserves here. + // + // Reserve the worst case reachable in the loop instead: every in-loop division has a + // divisor R ≤ rhs_len and a quotient split min(R/2, L-R) ≤ rhs_len/2, and + // mul::memory_requirement_up_to is non-decreasing in both arguments, so + // mul::memory_requirement_up_to(rhs_len, rhs_len/2) bounds them all (independent of + // lhs_len). It returns zero_layout for small operands automatically, since mul's own + // schoolbook threshold gates it. + let _ = lhs_len; + mul::memory_requirement_up_to(rhs_len, rhs_len / 2) } pub(crate) fn gcd_in_place( @@ -332,12 +344,14 @@ fn lehmer_ext_step( pub fn memory_requirement_ext_up_to(lhs_len: usize, rhs_len: usize) -> Layout { // Required memory: // - two numbers (t0 & t1) with at most the same size as lhs, add 1 buffer word - // - temporary space for a division (for euclidean step), and later a mulitplication (for coeff update) + // - temporary space for a division (for euclidean step, see memory_requirement_up_to + // for why the loop-wide bound rather than the initial lengths), and later a + // multiplication (for coeff update) let t_words = 2 * lhs_len + 2; memory::add_layout( memory::array_layout::(t_words), memory::max_layout( - div::memory_requirement_exact(lhs_len, rhs_len), // + mul::memory_requirement_up_to(rhs_len, rhs_len / 2), // worst-case in-loop division mul::memory_requirement_up_to(lhs_len, lhs_len / 2), // for coeff update ), ) diff --git a/integer/tests/convert.rs b/integer/tests/convert.rs index 568e0dc8..1d1775a9 100644 --- a/integer/tests/convert.rs +++ b/integer/tests/convert.rs @@ -434,6 +434,28 @@ fn test_to_f64() { ); } +/// Regression test: a magnitude at `DoubleWord::MAX` rounds up to the next power +/// of two as f64, and the saturating `f as DoubleWord` round-trip used to clamp +/// back to `DoubleWord::MAX` and report that inexact conversion as `Exact`. +#[cfg(target_pointer_width = "64")] +#[test] +fn test_to_f64_dword_max_boundary() { + // 2^128 - 1 == DoubleWord::MAX on 64-bit Word targets; it fits a single DoubleWord, + // so this exercises the `to_f64_small` fast path. + let two_pow_128: f64 = 2f64.powi(128); + let max = ubig!(340282366920938463463374607431768211455); + assert_eq!(max.to_f64(), Inexact(two_pow_128, Positive)); + let neg_max = -ibig!(340282366920938463463374607431768211455); + assert_eq!(neg_max.to_f64(), Inexact(-two_pow_128, Negative)); + + // One below the boundary still rounds up to 2^128; the guard must report it Inexact. + assert_eq!((max - ubig!(1)).to_f64(), Inexact(two_pow_128, Positive)); + + // A clearly representable value below the saturation zone is still Exact — the + // guard must not be over-eager. + assert_eq!((ubig!(1) << 127).to_f64(), Exact(2f64.powi(127))); +} + #[cfg(not(any(target_pointer_width = "16", force_bits = "16")))] #[test] fn test_from_u64_i64_const() { diff --git a/integer/tests/gcd.rs b/integer/tests/gcd.rs index e52639bd..03ca7f51 100644 --- a/integer/tests/gcd.rs +++ b/integer/tests/gcd.rs @@ -166,3 +166,31 @@ fn test_gcd_0() { fn test_gcd_ext_0() { let _ = ubig!(0).gcd_ext(ubig!(0)); } + +#[test] +fn test_gcd_large_lopsided_reduction() { + // Regression test: the gcd scratchpad is reserved once from the *initial* operand + // lengths, but each euclidean step's division dispatches on the *current* lengths. + // Before the fix, two large, similarly-sized operands (zero initial scratch) that + // later reduced through a lopsided euclidean step (divisor > 48 words and quotient + // > 32 words, taking the Burnikel-Ziegler + Karatsuba path) aborted with + // "internal error: not enough memory allocated" in both debug and release. + use dashu_int::{UBig, Word}; + + // Build a reduction that starts from similar-sized operands (so the initial division + // is "simple" and reserves zero scratch) but reaches a lopsided Burnikel-Ziegler + // step. Let R (~100 words), Q (~50 words), L = Q*R + 1 (~150 words); then + // gcd(L+R, L): step 1 (L+R)/L has quotient 1 -> next pair (L, R) + // step 2 L/R has quotient Q (>32 words) with divisor R (>48 words) + // -> Burnikel-Ziegler + Karatsuba, the under-reserved step. + // R and Q are fixed all-ones values chosen only for their size; the path is selected + // by operand length, not value. + let r = (UBig::from(1u8) << (100 * Word::BITS as usize)) - ubig!(1); // ~100 words + let q = (UBig::from(1u8) << (50 * Word::BITS as usize)) - ubig!(1); // ~50 words + let l = &q * &r + ubig!(1); + let a = &l + &r; // L + R + let b = l; // L + + // gcd(L+R, L) = gcd(L, R) = gcd(Q*R+1, R) = gcd(1, R) = 1. + assert_eq!((&a).gcd(&b), ubig!(1)); +} diff --git a/rational/CHANGELOG.md b/rational/CHANGELOG.md index 62bcef9c..d6075d45 100644 --- a/rational/CHANGELOG.md +++ b/rational/CHANGELOG.md @@ -1,5 +1,12 @@ # Changelog +## Unreleased + +### Fix +- Constructing or reducing an `RBig` with two large, similarly-sized operands no longer aborts with + "internal error: not enough memory allocated". The underlying integer GCD scratchpad is now sized + for the worst-case division reachable during reduction (see the `dashu-int` changelog). + ## 0.5.0 ### Fix diff --git a/rational/tests/arith_prop.rs b/rational/tests/arith_prop.rs index f9192489..acc3bf57 100644 --- a/rational/tests/arith_prop.rs +++ b/rational/tests/arith_prop.rs @@ -59,3 +59,28 @@ proptest! { prop_assert_eq!(reduced.as_relaxed(), scaled.as_relaxed()); } } + +/// Regression test: constructing an `RBig` reduces numerator/denominator by their +/// gcd, which must not abort with "not enough memory allocated" when that gcd +/// reduction hits a lopsided Burnikel-Ziegler division. The gcd scratchpad is sized +/// from the *initial* operand lengths, but each step dispatches on the *current* +/// lengths, so two large, similarly-sized operands (zero initial scratch) that later +/// reduce through a wide quotient under-reserved and panicked. +#[test] +fn rbig_reduce_large_lopsided() { + use dashu_int::Word; + + // L = Q*R + 1, then gcd(L+R, L) reduces through L/R (huge quotient Q, large divisor R). + // R and Q are fixed all-ones values chosen only for their size; the reduction path is + // selected by operand length, not value. + let r = (UBig::from(1u8) << (100 * Word::BITS as usize)) - UBig::from(1u8); // ~100 words + let q = (UBig::from(1u8) << (50 * Word::BITS as usize)) - UBig::from(1u8); // ~50 words + let l = &q * &r + UBig::from(1u64); + let num = IBig::from(&l + &r); + let den = l; + + let v = RBig::from_parts(num.clone(), den.clone()); + // gcd(L+R, L) = gcd(L, R) = gcd(Q*R+1, R) = gcd(1, R) = 1, so it is already reduced. + assert_eq!(v.numerator(), &num); + assert_eq!(v.denominator(), &den); +}