mirror of
https://github.com/RGBCube/uutils-coreutils
synced 2025-07-28 11:37:44 +00:00
Merge pull request #8028 from drinkcat/seq-latest-pow-with-context
uucore: num_parser: Update pow_with_context
This commit is contained in:
Daniel Hofstetter
GitHub
commit
fd37324d72
1 changed files with 71 additions and 32 deletions
|
|
@ -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,69 @@ 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.
|
||||
// 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();
|
||||
}
|
||||
|
||||
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.
|
||||
// (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;
|
||||
|
||||
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 +485,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() {
|
||||
// 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 {
|
||||
|
|
@ -960,14 +999,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")
|
||||
|
|
@ -975,21 +1014,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
|
||||
))
|
||||
|
|
|
|||
Loading…
Add table
Add a link
Reference in a new issue