uucore: num_parser: Update pow_with_context

This is the latest version in https://github.com/akubera/bigdecimal-rs/pull/148 ,
but the discussion sort of stalled, this is really complicated math,
but this new function is a little bit better than the original
(at least I hope so).

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

@@ -5,7 +5,9 @@
 
 //! Utilities for parsing numbers in various formats
 
-// spell-checker:ignore powf copysign prec inity infinit infs bigdecimal extendedbigdecimal biguint underflowed
+// spell-checker:ignore powf copysign prec ilog inity infinit infs bigdecimal extendedbigdecimal biguint underflowed muls
+
+use std::num::NonZeroU64;
 
 use bigdecimal::{
     BigDecimal, Context,
@@ -375,30 +377,66 @@ 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 {
+// Compute bd**exp using exponentiation by squaring algorithm, while maintaining the
+// precision specified in ctx (the number of digits would otherwise explode).
+//
+// Algorithm comes from https://en.wikipedia.org/wiki/Exponentiation_by_squaring
+//
+// TODO: Still pending discussion in https://github.com/akubera/bigdecimal-rs/issues/147,
+// we do lose a little bit of precision, and the last digits may not be correct.
+fn pow_with_context(bd: &BigDecimal, exp: i64, ctx: &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())
+    // When performing a multiplication between 2 numbers, we may lose up to 2 digits
+    // of precision.
+    // "Proof": https://github.com/akubera/bigdecimal-rs/issues/147#issuecomment-2793431202
+    const MARGIN_PER_MUL: u64 = 2;
+    // When doing many multiplication, we still introduce additional errors, add 1 more digit
+    // per 10 multiplications.
+    const MUL_PER_MARGIN_EXTRA: u64 = 10;
+
+    fn trim_precision(bd: BigDecimal, ctx: &Context, margin: u64) -> BigDecimal {
+        let prec = ctx.precision().get() + margin;
+        if bd.digits() > prec {
+            bd.with_precision_round(NonZeroU64::new(prec).unwrap(), 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)
+    // Count the number of multiplications we're going to perform, one per "1" binary digit
+    // in exp, and the number of times we can divide exp by 2.
+    let mut n = exp.unsigned_abs();
+    // Note: 63 - n.leading_zeros() == n.ilog2, but that's only available in recent Rust versions.
+    let muls = (n.count_ones() + (63 - n.leading_zeros()) - 1) as u64;
+    // Note: div_ceil would be nice to use here, but only available in recent Rust versions.
+    let margin_extra = (muls + MUL_PER_MARGIN_EXTRA / 2) / MUL_PER_MARGIN_EXTRA;
+    let mut margin = margin_extra + MARGIN_PER_MUL * muls;
+
+    let mut bd_y: BigDecimal = 1.into();
+    let mut bd_x = if exp >= 0 {
+        bd.clone()
     } else {
-        &bd * pow_with_context(bd.square(), (exp - 1) / 2, ctx)
+        bd.inverse_with_context(&ctx.with_precision(
+            NonZeroU64::new(ctx.precision().get() + margin + MARGIN_PER_MUL).unwrap(),
+        ))
     };
-    trim_precision(ret, ctx)
+
+    while n > 1 {
+        if n % 2 == 1 {
+            bd_y = trim_precision(&bd_x * bd_y, ctx, margin);
+            margin -= MARGIN_PER_MUL;
+            n -= 1;
+        }
+        bd_x = trim_precision(bd_x.square(), ctx, margin);
+        margin -= MARGIN_PER_MUL;
+        n /= 2;
+    }
+    debug_assert_eq!(margin, margin_extra);
+
+    trim_precision(bd_x * bd_y, ctx, 0)
 }
 
 // Construct an ExtendedBigDecimal based on parsed data
@@ -444,22 +482,20 @@ fn construct_extended_big_decimal<'a>(
         let bd = BigDecimal::from_bigint(signed_digits, 0)
             / BigDecimal::from_bigint(BigInt::from(16).pow(scale as u32), 0);
 
-        let abs_exponent = exponent.abs();
-        // Again, pow "only" supports u32 values. Just overflow/underflow if the value provided
-        // is > 2**32 or < 2**-32.
-        if abs_exponent > u32::MAX.into() {
-            return Err(make_error(exponent.is_positive(), negative));
-        }
+        // pow_with_context "only" supports i64 values. Just overflow/underflow if the value provided
+        // is > 2**64 or < 2**-64.
+        let exponent = match exponent.to_i64() {
+            Some(exp) => exp,
+            None => {
+                return Err(make_error(exponent.is_positive(), negative));
+            }
+        };
 
         // Confusingly, exponent is in base 2 for hex floating point numbers.
+        let base: BigDecimal = 2.into();
         // Note: We cannot overflow/underflow BigDecimal here, as we will not be able to reach the
         // maximum/minimum scale (i64 range).
-        let base: BigDecimal = if !exponent.is_negative() {
-            2.into()
-        } else {
-            BigDecimal::from(2).inverse()
-        };
-        let pow2 = pow_with_context(base, abs_exponent.to_u32().unwrap(), &Context::default());
+        let pow2 = pow_with_context(&base, exponent, &Context::default());
 
         bd * pow2
     } else {