From 4a0c2a95e804a8213f3d73cf664baf03a43fe157 Mon Sep 17 00:00:00 2001 From: Nicolas Boichat Date: Thu, 29 May 2025 12:07:01 +0200 Subject: [PATCH 1/3] 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 { From 5e3284139c6f280b78a2eb72e1a3bb22fb5c71ec Mon Sep 17 00:00:00 2001 From: Nicolas Boichat Date: Thu, 29 May 2025 14:37:03 +0200 Subject: [PATCH 2/3] num_parser: Fix tests after pow_with_context update We get more precision, and more range, now. --- src/uucore/src/lib/features/parser/num_parser.rs | 12 ++++++------ 1 file changed, 6 insertions(+), 6 deletions(-) diff --git a/src/uucore/src/lib/features/parser/num_parser.rs b/src/uucore/src/lib/features/parser/num_parser.rs index efd3ba0ed..c657b3ade 100644 --- a/src/uucore/src/lib/features/parser/num_parser.rs +++ b/src/uucore/src/lib/features/parser/num_parser.rs @@ -996,14 +996,14 @@ mod tests { assert_eq!( Ok(ExtendedBigDecimal::BigDecimal( // Wolfram Alpha says 9.8162042336235053508313854078782835648991393286913072670026492205522618203568834202759669215027003865... × 10^903089986 - BigDecimal::from_str("9.816204233623505350831385407878283564899139328691307267002649220552261820356883420275966921514831318e+903089986").unwrap() + BigDecimal::from_str("9.816204233623505350831385407878283564899139328691307267002649220552261820356883420275966921502700387e+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() + BigDecimal::from_str("1.349213146236998355103608893554488871595951104574239597804963176857050954139054064644219311222652032e-9030900").unwrap() )), // Couldn't get a answer from Wolfram Alpha for smaller negative exponents ExtendedBigDecimal::extended_parse("0x1p-30000000") @@ -1011,21 +1011,21 @@ mod tests { // ExtendedBigDecimal overflow/underflow assert!(matches!( - ExtendedBigDecimal::extended_parse(&format!("0x1p{}", u32::MAX as u64 + 1)), + ExtendedBigDecimal::extended_parse(&format!("0x1p{}", u64::MAX as u128 + 1)), Err(ExtendedParserError::Overflow(ExtendedBigDecimal::Infinity)) )); assert!(matches!( - ExtendedBigDecimal::extended_parse(&format!("-0x100P{}", u32::MAX as u64 + 1)), + ExtendedBigDecimal::extended_parse(&format!("-0x100P{}", u64::MAX as u128 + 1)), Err(ExtendedParserError::Overflow( ExtendedBigDecimal::MinusInfinity )) )); assert!(matches!( - ExtendedBigDecimal::extended_parse(&format!("0x1p-{}", u32::MAX as u64 + 1)), + ExtendedBigDecimal::extended_parse(&format!("0x1p-{}", u64::MAX as u128 + 1)), Err(ExtendedParserError::Underflow(ebd)) if ebd == ExtendedBigDecimal::zero() )); assert!(matches!( - ExtendedBigDecimal::extended_parse(&format!("-0x0.100p-{}", u32::MAX as u64 + 1)), + ExtendedBigDecimal::extended_parse(&format!("-0x0.100p-{}", u64::MAX as u128 + 1)), Err(ExtendedParserError::Underflow( ExtendedBigDecimal::MinusZero )) From c5b445f6f20fccb6aa3c56e044a2dc6f0492dead Mon Sep 17 00:00:00 2001 From: Nicolas Boichat Date: Sun, 1 Jun 2025 19:27:06 +0200 Subject: [PATCH 3/3] uucore: num_parser: Clarify origin of pow_with_context And why we use an older minimum Rust version in that piece of code. --- src/uucore/src/lib/features/parser/num_parser.rs | 3 +++ 1 file changed, 3 insertions(+) diff --git a/src/uucore/src/lib/features/parser/num_parser.rs b/src/uucore/src/lib/features/parser/num_parser.rs index c657b3ade..84aa82bdd 100644 --- a/src/uucore/src/lib/features/parser/num_parser.rs +++ b/src/uucore/src/lib/features/parser/num_parser.rs @@ -384,6 +384,8 @@ fn make_error<'a>(overflow: bool, negative: bool) -> ExtendedParserError<'a, Ext // // 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. +// Note: This has been copied from the latest revision in https://github.com/akubera/bigdecimal-rs/pull/148, +// so it's using minimum Rust version of `bigdecimal-rs`. fn pow_with_context(bd: &BigDecimal, exp: i64, ctx: &Context) -> BigDecimal { if exp == 0 { return 1.into(); @@ -412,6 +414,7 @@ fn pow_with_context(bd: &BigDecimal, exp: i64, ctx: &Context) -> BigDecimal { // 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. + // (see note above about minimum Rust version in use) let margin_extra = (muls + MUL_PER_MARGIN_EXTRA / 2) / MUL_PER_MARGIN_EXTRA; let mut margin = margin_extra + MARGIN_PER_MUL * muls;