Fix powi spurious overflow for bases near 1 on 32-bit (#95)

float/CHANGELOG.md

@@ -8,6 +8,14 @@
   it, and the complex literal macro uses it.
 
 ### Fix
+- `powi`'s overflow guard no longer reports a spurious overflow when the base is very close to 1
+  (a large significand with a large negative exponent). It estimated `log2(base)` with `log2_est`,
+  which for such a base is the difference of two ~1e3-magnitude terms and catastrophically cancels
+  to ~1e-4 of `f32` noise; scaled by a large exponent that noise crossed the overflow threshold
+  — which on 32-bit targets is only `isize::MAX·log2(B) ≈ 7e9` — and returned a spurious ±inf. This
+  made high-precision base conversion (`FBig::with_base`) panic ("arithmetic operations with the
+  infinity are not allowed!") on 32-bit targets (wasm32, i686). The guard now uses the
+  bit-length-based `log2_bounds`, which does not cancel (#95).
 - `to_f64`/`to_f32` now round the source once, directly to the target's precision at its own
   magnitude (fewer than 53/24 bits for subnormals), instead of through a fixed 53/24-bit
   intermediate that re-rounds into the subnormal grid. This removes a 1-ULP double-rounding error

float/src/convert.rs

@@ -1043,4 +1043,46 @@ mod tests {
         let int_val = FBig::<HalfAway, 10>::new(Repr::new(IBig::from(42), 0), Context::new(8));
         assert_eq!(IBig::try_from(int_val), Ok(IBig::from(42)));
     }
+
+    #[test]
+    fn with_base_high_precision_no_overflow() {
+        // Regression for issue #95: converting a high-precision base-2 float to base
+        // 10 panicked on 32-bit targets ("arithmetic operations with the infinity are
+        // not allowed!"). The base conversion evaluates exp(r) as `sum^(B^n)` through
+        // `powi` with a huge exponent (B^n) on a base (sum) very close to 1; `powi`'s
+        // overflow guard estimated log2(base) with the catastrophically-canceling
+        // `log2_est`, and the ~1e-4 of f32 noise scaled by the exponent crossed the
+        // (much smaller on 32-bit) isize threshold, yielding a spurious ±inf that then
+        // panicked when shifted. See `powi` in exp.rs for the fix.
+        use crate::round::mode::Zero;
+        use core::str::FromStr;
+
+        // The reporter's input: -1.1111…0011 in binary (578 significant bits), written
+        // in the hex form dashu accepts for base-2 floats (`0x1.<hex>…`). The value is
+        // identical to the raw binary literal.
+        let num = FBig::<Zero, 2>::from_str(
+            "-0x1.fffdc8d645194a5a95df4be063472d4406dd096339dd7dc2a8527d208b3da7b9e5c36b4f49a7982cb2ad20a4e7e4c016f858fe8cddea011a6d01fe3823189c4ed4f57a7babc331498",
+        )
+        .unwrap();
+
+        // at the original 578-bit precision the conversion succeeds …
+        let a = num
+            .clone()
+            .with_precision(578)
+            .value()
+            .with_base::<10>()
+            .value();
+        assert!(a.repr().is_finite());
+        // … and so does a slightly higher precision (586), which panicked on 32-bit
+        // (wasm32 / i686). The result matches the value computed on 64-bit. Compared
+        // by value (FBig equality ignores context) rather than via string formatting,
+        // so this works under no_std too.
+        let b = num.with_precision(586).value().with_base::<10>().value();
+        assert!(b.repr().is_finite());
+        let expected = FBig::<Zero, 10>::from_str(
+            "-1.9999661944503703041843468850635057967553124154072485151176192294480158424234268438137612977886891381228704640656094986435381057574477216648567249609280392009533217665484389886",
+        )
+        .unwrap();
+        assert_eq!(b, expected);
+    }
 }

float/src/exp.rs

@@ -133,32 +133,50 @@ impl<R: Round> Context<R> {
             return Ok(repr.map(|v| FBig::new(v, *self)));
         }
 
-        // Guard against exponent overflow for astronomically large results: the result
-        // magnitude has log2 ≈ exp·log2(base); if that exceeds the isize exponent range,
-        // return ±inf (|base| > 1) or 0 (|base| < 1) instead of overflowing mid-computation.
-        let base_log2 = base.log2_est() as f64;
+        // Guard against exponent overflow/underflow for astronomically large results.
+        // The result magnitude has log2 ≈ exp·log2(base); if that leaves the isize
+        // exponent range, return ±inf (|base|>1) or signed 0 (|base|<1) instead of
+        // overflowing the exponent mid-computation.
+        //
+        // Use the *bounds* of log2(base), never the point estimate `log2_est`: when
+        // base is very close to 1 (a large significand with a large negative exponent),
+        // log2(base) is the difference of two ~1e3-magnitude terms and suffers
+        // catastrophic cancellation — `log2_est` returns ~1e-4 of f32 noise rather than
+        // ~0. Scaled by a large exponent that noise crosses the overflow threshold,
+        // which on 32-bit is only isize::MAX·log2(B) ≈ 7e9 (vs ≈3e19 on 64-bit), so the
+        // guard fires spuriously and returns ±inf — see issue #95. The bounds are
+        // derived from the exact significand bit length, so they don't cancel. Declare
+        // overflow only on the lower bound (no false positives) and underflow only on
+        // the upper bound (no false underflows); anything in between is computed.
+        let (base_log2_lb, base_log2_ub) = base.log2_bounds();
+        let base_log2_lb = base_log2_lb as f64;
+        let base_log2_ub = base_log2_ub as f64;
         let threshold = (isize::MAX as f64) * (B.log2_est() as f64);
         let exp_f64 = i64::try_from(&exp).ok().map(|e| e as f64);
         let overflows = match exp_f64 {
-            Some(e) => e * base_log2 > threshold,
-            None => base_log2 != 0.0, // exp doesn't fit i64: overflows unless |base| == 1
+            Some(e) => e * base_log2_lb > threshold,
+            None => base_log2_lb > 0.0, // exp doesn't fit i64: overflows iff base > 1
         };
         if overflows {
-            return if base_log2 > 0.0 {
-                Err(FpError::Overflow(if base.sign() == Sign::Negative {
-                    Sign::Negative
-                } else {
-                    Sign::Positive
-                }))
+            return Err(FpError::Overflow(if base.sign() == Sign::Negative {
+                Sign::Negative
             } else {
-                // |base| < 1 and exponent huge → underflow to signed zero
-                let underflow_sign = if base.sign() == Sign::Negative && exp.bit(0) {
-                    Sign::Negative
-                } else {
-                    Sign::Positive
-                };
-                Err(FpError::Underflow(underflow_sign))
+                Sign::Positive
+            }));
+        }
+        let underflows = !base.significand.is_zero()
+            && match exp_f64 {
+                Some(e) => e * base_log2_ub < -threshold,
+                None => base_log2_ub < 0.0, // exp doesn't fit i64: underflows iff base < 1
+            };
+        if underflows {
+            // |base| < 1 and exponent huge → underflow to signed zero
+            let underflow_sign = if base.sign() == Sign::Negative && exp.bit(0) {
+                Sign::Negative
+            } else {
+                Sign::Positive
             };
+            return Err(FpError::Underflow(underflow_sign));
         }
 
         let work_context = if self.is_limited() {