uucore: num_parser: Limit precision when computing 2**exp

src/uucore/src/lib/features/parser/num_parser.rs

@@ -8,7 +8,7 @@
 // spell-checker:ignore powf copysign prec inity infinit bigdecimal extendedbigdecimal biguint underflowed
 
 use bigdecimal::{
-    BigDecimal,
+    BigDecimal, Context,
     num_bigint::{BigInt, BigUint, Sign},
 };
 use num_traits::Signed;
@@ -286,6 +286,32 @@ fn make_error<'a>(overflow: bool, negative: bool) -> ExtendedParserError<'a, Ext
     }
 }
 
+/// Compute bd**exp using exponentiation by squaring algorithm, while maintaining the
+/// precision specified in ctx (the number of digits would otherwise explode).
+// TODO: We do lose a little bit of precision, and the last digits may not be correct.
+// TODO: Upstream this to bigdecimal-rs.
+fn pow_with_context(bd: BigDecimal, exp: u32, ctx: &bigdecimal::Context) -> BigDecimal {
+    if exp == 0 {
+        return 1.into();
+    }
+
+    fn trim_precision(bd: BigDecimal, ctx: &bigdecimal::Context) -> BigDecimal {
+        if bd.digits() > ctx.precision().get() {
+            bd.with_precision_round(ctx.precision(), ctx.rounding_mode())
+        } else {
+            bd
+        }
+    }
+
+    let bd = trim_precision(bd, ctx);
+    let ret = if exp % 2 == 0 {
+        pow_with_context(bd.square(), exp / 2, ctx)
+    } else {
+        &bd * pow_with_context(bd.square(), (exp - 1) / 2, ctx)
+    };
+    trim_precision(ret, ctx)
+}
+
 // Construct an ExtendedBigDecimal based on parsed data
 fn construct_extended_big_decimal<'a>(
     digits: BigUint,
@@ -339,13 +365,14 @@ fn construct_extended_big_decimal<'a>(
         // Confusingly, exponent is in base 2 for hex floating point numbers.
         // Note: We cannot overflow/underflow BigDecimal here, as we will not be able to reach the
         // maximum/minimum scale (i64 range).
-        let pow2 = BigDecimal::from_bigint(BigInt::from(2).pow(abs_exponent.to_u32().unwrap()), 0);
-
-        if !exponent.is_negative() {
-            bd * pow2
+        let base: BigDecimal = if !exponent.is_negative() {
+            2.into()
         } else {
-            bd / pow2
-        }
+            BigDecimal::from(2).inverse()
+        };
+        let pow2 = pow_with_context(base, abs_exponent.to_u32().unwrap(), &Context::default());
+
+        bd * pow2
     } else {
         // scale != 0, which means that integral_only is not set, so only base 10 and 16 are allowed.
         unreachable!();
@@ -771,6 +798,26 @@ mod tests {
             ExtendedBigDecimal::extended_parse("0xf.fffffffffffffffffffff")
         );
 
+        // Test very large exponents (they used to take forever as we kept all digits in the past)
+        // Wolfram Alpha can get us (close to?) these values with a bit of log trickery:
+        // 2**3000000000 = 10**log_10(2**3000000000) = 10**(3000000000 * log_10(2))
+        // TODO: We do lose a little bit of precision, and the last digits are not be correct.
+        assert_eq!(
+            Ok(ExtendedBigDecimal::BigDecimal(
+                // Wolfram Alpha says 9.8162042336235053508313854078782835648991393286913072670026492205522618203568834202759669215027003865... × 10^903089986
+                BigDecimal::from_str("9.816204233623505350831385407878283564899139328691307267002649220552261820356883420275966921514831318e+903089986").unwrap()
+            )),
+            ExtendedBigDecimal::extended_parse("0x1p3000000000")
+        );
+        assert_eq!(
+            Ok(ExtendedBigDecimal::BigDecimal(
+                // Wolfram Alpha says 1.3492131462369983551036088935544888715959511045742395978049631768570509541390540646442193112226520316... × 10^-9030900
+                BigDecimal::from_str("1.349213146236998355103608893554488871595951104574239597804963176857050954139054064644219311222656999e-9030900").unwrap()
+            )),
+            // Couldn't get a answer from Wolfram Alpha for smaller negative exponents
+            ExtendedBigDecimal::extended_parse("0x1p-30000000")
+        );
+
         // ExtendedBigDecimal overflow/underflow
         assert!(matches!(
             ExtendedBigDecimal::extended_parse(&format!("0x1p{}", u32::MAX as u64 + 1)),