From 4a0c2a95e804a8213f3d73cf664baf03a43fe157 Mon Sep 17 00:00:00 2001 From: Nicolas Boichat Date: Thu, 29 May 2025 12:07:01 +0200 Subject: [PATCH] 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/lib/features/parser/num_parser.rs | 88 +++++++++++++------ 1 file changed, 62 insertions(+), 26 deletions(-) diff --git a/src/uucore/src/lib/features/parser/num_parser.rs b/src/uucore/src/lib/features/parser/num_parser.rs index 3ee07e413..efd3ba0ed 100644 --- a/src/uucore/src/lib/features/parser/num_parser.rs +++ b/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 {