diff --git a/Cargo.lock b/Cargo.lock index 9245ced..c8aaa0b 100644 --- a/Cargo.lock +++ b/Cargo.lock @@ -8,7 +8,7 @@ checksum = "0c4b4d0bd25bd0b74681c0ad21497610ce1b7c91b1022cd21c80c6fbdd9476b0" [[package]] name = "bigdecimal" -version = "0.4.5" +version = "0.4.6" dependencies = [ "autocfg", "libm", diff --git a/Cargo.toml b/Cargo.toml index 52e5989..1f233ef 100644 --- a/Cargo.toml +++ b/Cargo.toml @@ -1,6 +1,6 @@ [package] name = "bigdecimal" -version = "0.4.5" +version = "0.4.6" authors = ["Andrew Kubera"] description = "Arbitrary precision decimal numbers" documentation = "https://docs.rs/bigdecimal" @@ -53,3 +53,6 @@ std = ["num-bigint/std", "num-integer/std", "num-traits/std"] # BENCH: [[bench]] # BENCH: name = "arithmetic" # BENCH: harness = false + +[lints.rust] +unexpected_cfgs = { level = "allow", check-cfg = ['cfg(no_track_caller)'] } diff --git a/benches/arithmetic.rs b/benches/arithmetic.rs index 04ff995..8f7af2e 100644 --- a/benches/arithmetic.rs +++ b/benches/arithmetic.rs @@ -331,14 +331,17 @@ pub fn criterion_benchmark(c: &mut Criterion) { }, criterion::BatchSize::SmallInput)); - c.bench_function( - "sqrt", - |b| b.iter_batched( - || { - random_decimal.next() - }, - |a| { - black_box(a.sqrt()); - }, - criterion::BatchSize::SmallInput)); + let mut sqrt_group = c.benchmark_group("sqrt"); + for prec in [10, 100, 1000, 10000] { + sqrt_group.bench_function( + format!("sqrt_prec_{prec}"), + |b| b.iter_batched( + || { + (bigdecimal::Context::default().with_precision(std::num::NonZeroU64::new(prec).unwrap()), random_decimal.next()) + }, + |(ctx, a)| { + black_box(a.sqrt_with_context(&ctx)); + }, + criterion::BatchSize::SmallInput)); + } } diff --git a/build.rs b/build.rs index 414eba9..a92fc24 100644 --- a/build.rs +++ b/build.rs @@ -15,6 +15,9 @@ fn main() { let ac = autocfg::new(); ac.emit_rustc_version(1, 70); + // slice::fill + ac.emit_rustc_version(1, 50); + // Option::zip ac.emit_rustc_version(1, 46); diff --git a/src/arithmetic/addition.rs b/src/arithmetic/addition.rs index 11f04b9..bf7136b 100644 --- a/src/arithmetic/addition.rs +++ b/src/arithmetic/addition.rs @@ -9,18 +9,12 @@ pub(crate) fn add_bigdecimals( mut b: BigDecimal, ) -> BigDecimal { if b.is_zero() { - if a.scale < b.scale { - a.int_val *= ten_to_the((b.scale - a.scale) as u64); - a.scale = b.scale; - } + a.extend_scale_to(b.scale); return a; } if a.is_zero() { - if b.scale < a.scale { - b.int_val *= ten_to_the((a.scale - b.scale) as u64); - b.scale = a.scale; - } + b.extend_scale_to(a.scale); return b; } @@ -51,23 +45,35 @@ fn add_aligned_bigdecimals( pub(crate) fn add_bigdecimal_refs<'a, 'b, Lhs, Rhs>( lhs: Lhs, rhs: Rhs, + ctx: Option<&Context>, ) -> BigDecimal where Rhs: Into>, Lhs: Into>, { + use stdlib::cmp::Ordering::*; + let lhs = lhs.into(); let rhs = rhs.into(); if rhs.is_zero() { - return lhs.to_owned(); + let scale_diff = rhs.scale.saturating_sub(lhs.scale).max(0).min(15); + return lhs.to_owned_with_scale(lhs.scale + scale_diff); } if lhs.is_zero() { - return rhs.to_owned(); + let scale_diff = lhs.scale.saturating_sub(rhs.scale).max(0).min(15); + return rhs.to_owned_with_scale(rhs.scale + scale_diff); } - if lhs.scale >= rhs.scale { - lhs.to_owned() + rhs - } else { - rhs.to_owned() + lhs + + match lhs.scale.cmp(&rhs.scale) { + Equal => { + add_aligned_bigdecimal_ref_ref(lhs, rhs) + } + Greater => { + add_unaligned_bigdecimal_ref_ref(lhs, rhs, ctx) + } + Less => { + add_unaligned_bigdecimal_ref_ref(rhs, lhs, ctx) + } } } @@ -108,3 +114,40 @@ pub(crate) fn addassign_bigdecimal_ref<'a, T: Into>>( } } } + +/// Add BigDecimal references which have the same scale (integer addition) +fn add_aligned_bigdecimal_ref_ref( + lhs: BigDecimalRef, rhs: BigDecimalRef +) -> BigDecimal { + debug_assert!(!lhs.is_zero() && !rhs.is_zero()); + debug_assert_eq!(lhs.scale, rhs.scale); + + if lhs.digits.bits() >= rhs.digits.bits() { + lhs.to_owned() + rhs + } else { + rhs.to_owned() + lhs + } +} + +fn add_unaligned_bigdecimal_ref_ref( + lhs: BigDecimalRef, rhs: BigDecimalRef, _ctx: Option<&Context>, +) -> BigDecimal { + debug_assert!(!lhs.is_zero() && !rhs.is_zero()); + debug_assert!(lhs.scale >= rhs.scale); + + let scale_diff = (lhs.scale - rhs.scale) as u64; + + let shifted_rhs_digits = rhs.digits * ten_to_the_uint(scale_diff); + let shifted_rhs_int = BigInt::from_biguint(rhs.sign, shifted_rhs_digits); + let shifted_rhs = BigDecimal::new(shifted_rhs_int, lhs.scale); + + shifted_rhs + lhs +} + + +#[cfg(test)] +mod test { + use super::*; + + include!("addition.tests.rs"); +} diff --git a/src/arithmetic/addition.tests.rs b/src/arithmetic/addition.tests.rs new file mode 100644 index 0000000..0fe1d82 --- /dev/null +++ b/src/arithmetic/addition.tests.rs @@ -0,0 +1,52 @@ +mod add_bigdecimals { + use super::*; + use paste::paste; + + macro_rules! impl_case { + ( $name:ident: $a:literal + $b:literal = $c:literal ) => { + + #[test] + fn $name() { + let lhs: BigDecimal = $a.parse().unwrap(); + let rhs: BigDecimal = $b.parse().unwrap(); + + let l_plus_r = add_bigdecimals(lhs.clone(), rhs.clone()); + let r_plus_l = add_bigdecimals(rhs, lhs); + + let expected: BigDecimal = $c.parse().unwrap(); + assert_eq!(expected.int_val, l_plus_r.int_val); + assert_eq!(expected.scale, l_plus_r.scale); + + assert_eq!(expected.int_val, r_plus_l.int_val); + assert_eq!(expected.scale, r_plus_l.scale); + } + + paste! { + #[test] + fn [< $name _refs >]() { + let lhs: BigDecimal = $a.parse().unwrap(); + let rhs: BigDecimal = $b.parse().unwrap(); + + let l_plus_r = add_bigdecimal_refs(&lhs, &rhs, None); + let r_plus_l = add_bigdecimal_refs(&rhs, &lhs, None); + + let expected: BigDecimal = $c.parse().unwrap(); + assert_eq!(expected.int_val, l_plus_r.int_val); + assert_eq!(expected.scale, l_plus_r.scale); + + assert_eq!(expected.int_val, r_plus_l.int_val); + assert_eq!(expected.scale, r_plus_l.scale); + } + } + }; + } + + impl_case!(case_1d2345_123d45: "1.2345" + "123.45" = "124.6845"); + impl_case!(case_123d43e5_1d2345: "123.43e5" + "1.2345" = "12343001.2345"); + + impl_case!(case_0_0: "0" + "0" = "0"); + impl_case!(case_0_0d00: "0" + "0.00" = "0.00"); + impl_case!(case_10_0d00: "10" + "0.00" = "10.00"); + impl_case!(case_22132e2_0d0000: "22132e2" + "0.0000" = "2213200.0000"); + impl_case!(case_n316d79_0en6: "-316.79" + "0e-6" = "-316.790000"); +} diff --git a/src/arithmetic/cbrt.rs b/src/arithmetic/cbrt.rs index 8887fb9..a75b0cf 100644 --- a/src/arithmetic/cbrt.rs +++ b/src/arithmetic/cbrt.rs @@ -2,85 +2,107 @@ use crate::*; use num_bigint::BigUint; +use rounding::NonDigitRoundingData; +use stdlib::num::NonZeroU64; -/// implementation of cuberoot - always positive -pub(crate) fn impl_cbrt_uint_scale(n: &BigUint, scale: i64, ctx: &Context) -> BigDecimal { - // make guess based on number of bits in the number - let guess = make_cbrt_guess(n.bits(), scale); +pub(crate) fn impl_cbrt_int_scale(n: &BigInt, scale: i64, ctx: &Context) -> BigDecimal { + let rounding_data = NonDigitRoundingData { + sign: n.sign(), + mode: ctx.rounding_mode(), + }; - let three = BigInt::from(3); + impl_cbrt_uint_scale((n.magnitude(), scale).into(), ctx.precision(), rounding_data) +} - let n = BigInt::from_biguint(Sign::Plus, n.clone()); +/// implementation of cuberoot - always positive +pub(crate) fn impl_cbrt_uint_scale( + n: WithScale<&BigUint>, + precision: NonZeroU64, + // contains sign and rounding mode + rounding_data: NonDigitRoundingData, +) -> BigDecimal { + if n.is_zero() { + let biguint = BigInt::from_biguint(Sign::Plus, n.value.clone()); + return BigDecimal::new(biguint, n.scale / 3); + } - let max_precision = ctx.precision().get(); + // count number of digits in the decimal + let integer_digit_count = count_decimal_digits_uint(n.value); - let next_iteration = move |r: BigDecimal| { - let sqrd = r.square(); - let tmp = impl_division( - n.clone(), - &sqrd.int_val, - scale - sqrd.scale, - max_precision + 1, - ); - let tmp = tmp + r.double(); - impl_division(tmp.int_val, &three, tmp.scale, max_precision + 3) - }; + // extra digits to use for rounding + let extra_rounding_digit_count = 4; - // result initial - let mut running_result = next_iteration(guess); + // required number of digits for precision and rounding + let required_precision = precision.get() + extra_rounding_digit_count; + let required_precision = 3 * required_precision; - let mut prev_result = BigDecimal::one(); - let mut result = BigDecimal::zero(); + // number of extra zeros to add to end of integer_digits + let mut exp_shift = required_precision.saturating_sub(integer_digit_count); - // TODO: Prove that we don't need to arbitrarily limit iterations - // and that convergence can be calculated - while prev_result != result { - // store current result to test for convergence - prev_result = result; + // effective scale after multiplying by 10^exp_shift + // (we've added that many trailing zeros after) + let shifted_scale = n.scale + exp_shift as i64; - running_result = next_iteration(running_result); + let (mut new_scale, remainder) = shifted_scale.div_rem(&3); - // result has clipped precision, running_result has full precision - result = if running_result.digits() > max_precision { - running_result.with_precision_round(ctx.precision(), ctx.rounding_mode()) - } else { - running_result.clone() - }; + if remainder > 0 { + new_scale += 1; + exp_shift += (3 - remainder) as u64; + } else if remainder < 0 { + exp_shift += remainder.neg() as u64; } - return result; -} + // clone-on-write copy of digits + let mut integer_digits = stdlib::borrow::Cow::Borrowed(n.value); -/// Find initial cbrt guess based on number of bits in integer and the scale -/// -/// ```math -/// 2^bit_count * 10^-scale <= *n* < 2^(bit_count+1) * 10^-scale -/// -/// cbrt(n2^bit_count * 10^-scale) -/// cbrt(2^bit_count * 10^-scale) -/// => Exp2[1/3 * Log2[2^bit_count * 10^-scale]] -/// => Exp2[1/3 * (bit_count - scale * Log2[10]] -/// ``` -/// -fn make_cbrt_guess(bit_count: u64, scale: i64) -> BigDecimal { - // weight of cube root average above minimum within range: 3/4*2^(4/3) - let magic_guess_scale = 1.1398815748423097_f64; - - let bit_count = bit_count as f64; - let scale = scale as f64; - - let initial_guess = (bit_count - scale * LOG2_10) / 3.0; - let res = magic_guess_scale * exp2(initial_guess); - - match BigDecimal::try_from(res) { - Ok(res) => res, - Err(_) => { - // can't guess with float - just guess magnitude - let scale = (scale - bit_count / LOG2_10).round() as i64; - BigDecimal::new(BigInt::from(1), scale / 3) - } + // add required trailing zeros to integer_digits + if exp_shift > 0 { + arithmetic::multiply_by_ten_to_the_uint( + integer_digits.to_mut(), exp_shift + ); + } + + let result_digits = integer_digits.nth_root(3); + let result_digits_count = count_decimal_digits_uint(&result_digits); + debug_assert!(result_digits_count > precision.get()); + + let digits_to_trim = result_digits_count - precision.get(); + debug_assert_ne!(digits_to_trim, 0); + debug_assert!((result_digits_count as i64 - count_decimal_digits_uint(&integer_digits) as i64 / 3).abs() < 2); + + new_scale -= digits_to_trim as i64; + + let divisor = ten_to_the_uint(digits_to_trim); + let (mut result_digits, remainder) = result_digits.div_rem(&divisor); + + let remainder_digits = remainder.to_radix_le(10); + let insig_digit0; + let trailing_digits; + if remainder_digits.len() < digits_to_trim as usize { + // leading zeros + insig_digit0 = 0; + trailing_digits = remainder_digits.as_slice(); + } else { + let (&d, rest) = remainder_digits.split_last().unwrap(); + insig_digit0 = d; + trailing_digits = rest; } + + let insig_data = rounding::InsigData::from_digit_and_lazy_trailing_zeros( + rounding_data, insig_digit0, || { trailing_digits.iter().all(Zero::is_zero) } + ); + + // lowest digit to round + let sig_digit = (&result_digits % 10u8).to_u8().unwrap(); + let rounded_digit = insig_data.round_digit(sig_digit); + + let rounding_term = rounded_digit - sig_digit; + result_digits += rounding_term; + + let result = BigInt::from_biguint(rounding_data.sign, result_digits); + + BigDecimal::new(result, new_scale) } @@ -89,46 +111,72 @@ mod test { use super::*; use stdlib::num::NonZeroU64; - #[test] - fn test_cbrt() { - let vals = vec![ - ("0.00", "0"), - ("1.00", "1"), - ("1.001", "1.000333222283909495175449559955220102010284758197360454054345461242739715702641939155238095670636841"), - ("10", "2.154434690031883721759293566519350495259344942192108582489235506346411106648340800185441503543243276"), - ("13409.179789484375", "23.7575"), - ("-59283293e25", "-84006090355.84281237113712383191213626687332139035750444925827809487776780721673264524620270275301685"), - ("94213372931e-127", "2.112049945275324414051072540210070583697242797173805198575907094646677475250362108901530353886613160E-39"), - ]; - for &(x, y) in vals.iter() { - let a = BigDecimal::from_str(x).unwrap().cbrt(); - let b = BigDecimal::from_str(y).unwrap(); - assert_eq!(a, b); - } + macro_rules! impl_test { + ($name:ident; $input:literal => $expected:literal) => { + #[test] + fn $name() { + let n: BigDecimal = $input.parse().unwrap(); + let value = n.cbrt(); + + let expected = $expected.parse().unwrap(); + assert_eq!(value, expected); + } + }; + ($name:ident; prec=$prec:literal; round=$round:ident; $input:literal => $expected:literal) => { + #[test] + fn $name() { + let ctx = Context::new(NonZeroU64::new($prec).unwrap(), RoundingMode::$round); + let n: BigDecimal = $input.parse().unwrap(); + let value = n.cbrt_with_context(&ctx); + + let expected = $expected.parse().unwrap(); + assert_eq!(value, expected); + } + }; } + mod default { + use super::*; - #[test] - fn test_cbrt_prec15() { - let vals = vec![ - ("0.00", "0"), - ("1.00", "1"), - ("1.001", "1.00033322228390"), - ("10", "2.15443469003188"), - ("13409.179789484375", "23.7575"), - ("-59283293e25", "-84006090355.8428"), - ("94213372931e-127", "2.11204994527532E-39"), - ]; - - let ctx = Context::new(NonZeroU64::new(15).unwrap(), RoundingMode::Down); - - for &(x, y) in vals.iter() { - let a = BigDecimal::from_str(x).unwrap().cbrt_with_context(&ctx); - let b = y.parse().unwrap(); - assert_eq!(a, b); - } + impl_test!(case_0; "0.00" => "0"); + impl_test!(case_1; "1.00" => "1"); + impl_test!(case_1d001; "1.001" => "1.000333222283909495175449559955220102010284758197360454054345461242739715702641939155238095670636841"); + impl_test!(case_10; "10" => "2.154434690031883721759293566519350495259344942192108582489235506346411106648340800185441503543243276"); + impl_test!(case_13409d179789484375; "13409.179789484375" => "23.7575"); + impl_test!(case_n59283293e25; "-59283293e25" => "-84006090355.84281237113712383191213626687332139035750444925827809487776780721673264524620270275301685"); + impl_test!(case_94213372931en127; "94213372931e-127" => "2.112049945275324414051072540210070583697242797173805198575907094646677475250362108901530353886613160E-39"); } + impl_test!(case_prec15_down_10; prec=15; round=Down; "10" => "2.15443469003188"); + impl_test!(case_prec6_up_0d979970546636727; prec=6; round=Up; "0.979970546636727" => "0.993279"); + + impl_test!(case_1037d495615705321421375_full; "1037.495615705321421375" => "10.123455"); + impl_test!(case_1037d495615705321421375_prec7_halfdown; prec=7; round=HalfDown; "1037.495615705321421375" => "10.12345"); + impl_test!(case_1037d495615705321421375_prec7_halfeven; prec=7; round=HalfEven; "1037.495615705321421375" => "10.12346"); + impl_test!(case_1037d495615705321421375_prec7_halfup; prec=7; round=HalfUp; "1037.495615705321421375" => "10.12346"); + + impl_test!(case_0d014313506928855520728400001_full; "0.014313506928855520728400001" => "0.242800001"); + impl_test!(case_0d014313506928855520728400001_prec6_down; prec=6; round=Down; "0.014313506928855520728400001" => "0.242800"); + impl_test!(case_0d014313506928855520728400001_prec6_up; prec=6; round=Up; "0.014313506928855520728400001" => "0.242801"); + + impl_test!(case_4151902e20_prec16_halfup; prec=16; round=HalfUp; "4151902e20" => "746017527.6855992"); + impl_test!(case_4151902e20_prec16_up; prec=16; round=Up; "4151902e20" => "746017527.6855993"); + impl_test!(case_4151902e20_prec17_up; prec=17; round=Up; "4151902e20" => "746017527.68559921"); + impl_test!(case_4151902e20_prec18_up; prec=18; round=Up; "4151902e20" => "746017527.685599209"); + // impl_test!(case_4151902e20_prec18_up; prec=18; round=Up; "4151902e20" => "746017527.685599209"); + + impl_test!(case_1850846e201_prec14_up; prec=16; round=Up; "1850846e201" => "1.227788123885769e69"); + + impl_test!(case_6d3797558642427987505823530913e85_prec16_up; prec=160; round=Up; "6.3797558642427987505823530913E+85" => "3995778017e19"); + + impl_test!(case_88573536600476899341824_prec20_up; prec=20; round=Up; "88573536600476899341824" => "44576024"); + impl_test!(case_88573536600476899341824_prec7_up; prec=7; round=Up; "88573536600476899341824" => "4457603e1"); + + impl_test!(case_833636d150970875_prec5_up; prec=5; round=Up; "833636.150970875" => "94.115000"); + impl_test!(case_833636d150970875_prec5_halfup; prec=5; round=HalfUp; "833636.150970875" => "94.115"); + impl_test!(case_833636d150970875_prec4_halfup; prec=4; round=HalfUp; "833636.150970875" => "94.12"); + impl_test!(case_833636d150970875_prec20_up; prec=20; round=Up; "833636.150970875" => "94.115000"); + #[cfg(property_tests)] mod prop { use super::*; diff --git a/src/arithmetic/mod.rs b/src/arithmetic/mod.rs index c6df8fc..e96efb3 100644 --- a/src/arithmetic/mod.rs +++ b/src/arithmetic/mod.rs @@ -62,6 +62,18 @@ pub(crate) fn ten_to_the_uint(pow: u64) -> BigUint { } } +pub(crate) fn multiply_by_ten_to_the_uint(n: &mut T, pow: P) + where T: MulAssign + MulAssign, + P: ToPrimitive +{ + let pow = pow.to_u64().expect("exponent overflow error"); + if pow < 20 { + *n *= 10u64.pow(pow as u32); + } else { + *n *= ten_to_the_uint(pow); + } + +} /// Return number of decimal digits in integer pub(crate) fn count_decimal_digits(int: &BigInt) -> u64 { @@ -123,3 +135,59 @@ pub(crate) fn diff_usize(a: usize, b: usize) -> (Ordering, usize) { Equal => (Equal, 0), } } + + +/// Add carry to given number, returning trimmed value and storing overflow back in carry +/// +pub(crate) fn add_carry(n: u8, carry: &mut u8) -> u8 { + let s = n + *carry; + if s < 10 { + *carry = 0; + s + } else { + debug_assert!(s < 20); + *carry = 1; + s - 10 + } +} + + +/// If n is greater than 10, split and store overflow in carry +/// +/// No action if n is less than 10. +/// +/// Carry is not allowed to be 1 if n is two digits +/// +pub(crate) fn store_carry(n: u8, carry: &mut u8) -> u8 { + if n < 10 { + n + } else { + debug_assert!(n < 20); + debug_assert_eq!(carry, &0); + *carry = 1; + n - 10 + } +} + + +/// Extend destination vector with values in D, adding carry while carry is not zero +/// +/// If carry overflows, it is NOT pushed into the destination vector. +/// +pub(crate) fn extend_adding_with_carry>( + dest: &mut Vec, + mut digits: D, + carry: &mut u8, +) { + while *carry != 0 { + match digits.next() { + Some(d) => { + dest.push(add_carry(d, carry)) + } + None => { + return; + } + } + } + dest.extend(digits); +} diff --git a/src/arithmetic/sqrt.rs b/src/arithmetic/sqrt.rs index d2a56df..d0c0ff8 100644 --- a/src/arithmetic/sqrt.rs +++ b/src/arithmetic/sqrt.rs @@ -3,239 +3,144 @@ use crate::*; -fn impl_division(mut num: BigUint, den: &BigUint, mut scale: i64, max_precision: u64) -> BigDecimal { - use Sign::Plus; - - // quick zero check - if num.is_zero() { - return BigDecimal::new(BigInt::from_biguint(Plus, num), 0); - } - - // shift digits until numerator is larger than denominator (set scale appropriately) - while num < *den { - scale += 1; - num *= 10u8; - } - - // first division - let (mut quotient, mut remainder) = num.div_rem(den); - - // division complete - if remainder.is_zero() { - return BigDecimal { - int_val: BigInt::from_biguint(Plus, quotient), - scale: scale, - }; - } - - let mut precision = count_decimal_digits_uint("ient); - - // shift remainder by 1 decimal; - // quotient will be 1 digit upon next division - remainder *= 10u8; - - while !remainder.is_zero() && precision < max_precision { - let (q, r) = remainder.div_rem(den); - quotient = quotient * 10u8 + q; - remainder = r * 10u8; - - precision += 1; - scale += 1; - } - - if !remainder.is_zero() { - // round final number with remainder - quotient += get_rounding_term_uint(&remainder.div(den)); - } - - BigDecimal::new(BigInt::from_biguint(Plus, quotient), scale) -} - -fn get_rounding_term_uint(num: &BigUint) -> u8 { - if num.is_zero() { - return 0; - } - - let digits = (num.bits() as f64 / LOG2_10) as u64; - let mut n = ten_to_the_uint(digits); - - // loop-method - loop { - if *num < n { - return 1; - } - n *= 5u8; - if *num < n { - return 0; - } - n *= 2u8; - } - - // string-method - // let s = format!("{}", num); - // let high_digit = u8::from_str(&s[0..1]).unwrap(); - // if high_digit < 5 { 0 } else { 1 } -} - pub(crate) fn impl_sqrt(n: &BigUint, scale: i64, ctx: &Context) -> BigDecimal { - // make guess - let guess = { - let magic_guess_scale = 1.1951678538495576_f64; - let initial_guess = (n.bits() as f64 - scale as f64 * LOG2_10) / 2.0; - let res = magic_guess_scale * exp2(initial_guess); - - if res.is_normal() { - BigDecimal::try_from(res).unwrap() - } else { - // can't guess with float - just guess magnitude - let scale = (n.bits() as f64 / -LOG2_10 + scale as f64).round() as i64; - BigDecimal::new(BigInt::from(1), scale / 2) - } - }; - - // // wikipedia example - use for testing the algorithm - // if self == &BigDecimal::from_str("125348").unwrap() { - // running_result = BigDecimal::from(600) - // } - - // TODO: Use context variable to set precision - let max_precision = ctx.precision().get(); - - let next_iteration = move |r: BigDecimal| { - // division needs to be precise to (at least) one extra digit - let tmp = impl_division( - n.clone(), - &r.int_val.magnitude(), - scale - r.scale, - max_precision + 1, - ); - - // half will increase precision on each iteration - (tmp + r).half() - }; - - // calculate first iteration - let mut running_result = next_iteration(guess); - - let mut prev_result = BigDecimal::one(); - let mut result = BigDecimal::zero(); - - // TODO: Prove that we don't need to arbitrarily limit iterations - // and that convergence can be calculated - while prev_result != result { - // store current result to test for convergence - prev_result = result; - - // calculate next iteration - running_result = next_iteration(running_result); - - // 'result' has clipped precision, 'running_result' has full precision - result = if running_result.digits() > max_precision { - running_result.with_prec(max_precision) - } else { - running_result.clone() - }; - } - - result + // Calculate the number of digits and the difference compared to the scale + let num_digits = count_decimal_digits_uint(&n); + let scale_diff = BigInt::from(num_digits) - scale; + + // Calculate the number of wanted digits and the exponent we need to raise the original value to + // We want twice as many digits as the precision because sqrt halves the number of digits + // We add an extra one for rounding purposes + let prec = ctx.precision().get(); + let extra_rounding_digit_count = 5; + let wanted_digits = 2 * (prec + extra_rounding_digit_count); + let exponent = wanted_digits.saturating_sub(num_digits) + u64::from(scale_diff.is_odd()); + let sqrt_digits = (n * ten_to_the_uint(exponent)).sqrt(); + + // Calculate the scale of the result + let result_scale_digits = 2 * (2 * prec - scale_diff) - 1; + let result_scale_decimal: BigDecimal = BigDecimal::new(result_scale_digits, 0) / 4.0; + let mut result_scale = result_scale_decimal.with_scale_round(0, RoundingMode::HalfEven).int_val; + + // Round the value so it has the correct precision requested + result_scale += count_decimal_digits_uint(&sqrt_digits).saturating_sub(prec); + let unrounded_result = BigDecimal::new(sqrt_digits.into(), result_scale.to_i64().unwrap()); + unrounded_result.with_precision_round(ctx.precision(), ctx.rounding_mode()) } - #[cfg(test)] mod test { use super::*; - #[test] - fn test_sqrt() { - let vals = vec![ - ("1e-232", "1e-116"), - ("1.00", "1"), - ("1.001", "1.000499875062460964823258287700109753027590031219780479551442971840836093890879944856933288426795152"), - ("100", "10"), - ("49", "7"), - (".25", ".5"), - ("0.0152399025", ".12345"), - ("152399025", "12345"), - (".00400", "0.06324555320336758663997787088865437067439110278650433653715009705585188877278476442688496216758600590"), - (".1", "0.3162277660168379331998893544432718533719555139325216826857504852792594438639238221344248108379300295"), - ("2", "1.414213562373095048801688724209698078569671875376948073176679737990732478462107038850387534327641573"), - ("125348", "354.0451948551201563108487193176101314241016013304294520812832530590100407318465590778759640828114535"), - ("18446744073709551616.1099511", "4294967296.000000000012799992691725492477397918722952224079252026972356303360555051219312462698703293"), - ("3.141592653589793115997963468544185161590576171875", "1.772453850905515992751519103139248439290428205003682302442979619028063165921408635567477284443197875"), - (".000000000089793115997963468544185161590576171875", "0.000009475922962855041517561783740144225422359796851494316346796373337470068631250135521161989831460407155"), - ("0.7177700109762963922745342343167413624881759290454997218753321040760896053150388903350654937434826216803814031987652326749140535150336357405672040727695124057298138872112244784753994931999476811850580200000000000000000000000000000000", "0.8472130847527653667042980517799020703921106560594525833177762276594388966885185567535692987624493813"), - ("0.01234567901234567901234567901234567901234567901234567901234567901234567901234567901234567901234567901", "0.1111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111"), - ("0.1108890000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000444", "0.3330000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000667"), - ]; - for &(x, y) in vals.iter() { - let a = BigDecimal::from_str(x).unwrap().sqrt().unwrap(); - let b = BigDecimal::from_str(y).unwrap(); - assert_eq!(a, b); - } + macro_rules! impl_case { + ($name:ident; $input:literal => $expected:literal) => { + #[test] + fn $name() { + let n: BigDecimal = $input.parse().unwrap(); + let value = n.sqrt().unwrap(); + + let expected = $expected.parse().unwrap(); + assert_eq!(value, expected); + // assert_eq!(value.scale, expected.scale); + } + }; + ($name:ident; prec=$prec:literal; round=$round:ident; $input:literal => $expected:literal) => { + #[test] + fn $name() { + let ctx = Context::default() + .with_prec($prec).unwrap() + .with_rounding_mode(RoundingMode::$round); + let n: BigDecimal = $input.parse().unwrap(); + let value = n.sqrt_with_context(&ctx).unwrap(); + + let expected = $expected.parse().unwrap(); + assert_eq!(value, expected); + // assert_eq!(value.scale, expected.scale); + } + }; } + impl_case!(case_0d000; "0.000" => "0"); + impl_case!(case_1en232; "1e-232" => "1e-116"); + impl_case!(case_1d00; "1.00" => "1.00"); + impl_case!(case_1d001; "1.001" => "1.000499875062460964823258287700109753027590031219780479551442971840836093890879944856933288426795152"); + impl_case!(case_100d0; "100" => "10.0000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000"); + impl_case!(case_49; "49" => "7.0000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000"); + impl_case!(case_d25; ".25" => ".5"); + impl_case!(case_0d0152399025; "0.0152399025" => ".12345"); + impl_case!(case_0d00400; "0.00400" => "0.06324555320336758663997787088865437067439110278650433653715009705585188877278476442688496216758600590"); + impl_case!(case_0d1; "0.1" => "0.3162277660168379331998893544432718533719555139325216826857504852792594438639238221344248108379300295"); + impl_case!(case_152399025; "152399025" => "12345"); + impl_case!(case_2; "2" => "1.414213562373095048801688724209698078569671875376948073176679737990732478462107038850387534327641573"); + impl_case!(case_125348; "125348" => "354.0451948551201563108487193176101314241016013304294520812832530590100407318465590778759640828114535"); + impl_case!(case_121d000242000121; "121.000242000121000000" => "11.000011000"); + impl_case!(case_0d01234567901234567901234567901234567901234567901234567901234567901234567901234567901234567901234567901; "0.01234567901234567901234567901234567901234567901234567901234567901234567901234567901234567901234567901" => "0.1111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111"); + impl_case!(case_2e70; "2e70" => "141421356237309504880168872420969807.8569671875376948073176679737990732478462107038850387534327641573"); + impl_case!(case_8d9793115997963468544185161590576171875en11; "8.9793115997963468544185161590576171875e-11" => "0.000009475922962855041517561783740144225422359796851494316346796373337470068631250135521161989831460407155"); + impl_case!(case_18446744073709551616d1099511; "18446744073709551616.1099511" => "4294967296.000000000012799992691725492477397918722952224079252026972356303360555051219312462698703293"); + + impl_case!(case_3d1415926; "3.141592653589793115997963468544185161590576171875" => "1.772453850905515992751519103139248439290428205003682302442979619028063165921408635567477284443197875"); + impl_case!(case_0d71777001; "0.7177700109762963922745342343167413624881759290454997218753321040760896053150388903350654937434826216803814031987652326749140535150336357405672040727695124057298138872112244784753994931999476811850580200000000000000000000000000000000" => "0.8472130847527653667042980517799020703921106560594525833177762276594388966885185567535692987624493813"); + impl_case!(case_0d110889ddd444; "0.1108890000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000444" => "0.3330000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000667"); + + impl_case!(case_3e170; "3e170" => "17320508075688772935274463415058723669428052538103806280558069794519330169088000370811.46186757248576"); + impl_case!(case_9e199; "9e199" => "9486832980505137995996680633298155601158665417975650480572514558377783315917714664032744325137900886"); + impl_case!(case_7e200; "7e200" => "2645751311064590590501615753639260425710259183082450180368334459201068823230283627760392886474543611e1"); + impl_case!(case_777e204; "777e204" => "2.787471972953270789531596912111625325974789615194854615319795902911796043681078997362635440358922503E+103"); + impl_case!(case_777e600; "7e600" => "2.645751311064590590501615753639260425710259183082450180368334459201068823230283627760392886474543611E+300"); + impl_case!(case_2e900; "2e900" => "1.414213562373095048801688724209698078569671875376948073176679737990732478462107038850387534327641573E+450"); + impl_case!(case_7e999; "7e999" => "8.366600265340755479781720257851874893928153692986721998111915430804187725943170098308147119649515362E+499"); + impl_case!(case_74908163946345982392040522594123773796e999; "74908163946345982392040522594123773796e999" => "2.736935584670307552030924971360722787091742391079630976117950955395149091570790266754718322365663909E+518"); + impl_case!(case_20e1024; "20e1024" => "4.472135954999579392818347337462552470881236719223051448541794490821041851275609798828828816757564550E512"); + impl_case!(case_3en1025; "3e-1025" => "5.477225575051661134569697828008021339527446949979832542268944497324932771227227338008584361638706258E-513"); + + impl_case!(case_3242053850483855en13_prec11_round_down; prec=11; round=Down; "324.2053850483855" => "18.005704236"); + impl_case!(case_3242053850483855en13_prec11_round_up; prec=11; round=Up; "324.2053850483855" => "18.005704237"); + impl_case!(case_3242053850483855en13_prec31_round_up; prec=31; round=Up; "324.2053850483855" => "18.00570423639090823994825477228"); + + impl_case!(case_5d085019992340351en10_prec25_round_down; prec=25; round=Down; "5.085019992340351e-10" => "0.00002254998889653906459324292"); + + impl_case!(case_3025d13579652399025_prec3_round_up; prec=3; round=Up; "3025.13579652399025" => "55.1"); + + impl_case!(case_3025d13579652399025_prec9_round_down; prec=9; round=Down; "3025.13579652399025" => "55.0012345"); + impl_case!(case_3025d13579652399025_prec9_round_up; prec=9; round=Up; "3025.13579652399025" => "55.0012345"); + + impl_case!(case_3025d13579652399025_prec8_round_halfdown; prec=8; round=HalfDown; "3025.13579652399025" => "55.001234"); + impl_case!(case_3025d13579652399025_prec8_round_halfeven; prec=8; round=HalfEven; "3025.13579652399025" => "55.001234"); + impl_case!(case_3025d13579652399025_prec8_round_halfup; prec=8; round=HalfUp; "3025.13579652399025" => "55.001235"); + #[test] - fn test_big_sqrt() { - use num_bigint::BigInt; + fn test_sqrt_rounding() { let vals = vec![ - (("2", -70), "141421356237309504880168872420969807.8569671875376948073176679737990732478462107038850387534327641573"), - (("3", -170), "17320508075688772935274463415058723669428052538103806280558069794519330169088000370811.46186757248576"), - (("9", -199), "9486832980505137995996680633298155601158665417975650480572514558377783315917714664032744325137900886"), - (("7", -200), "26457513110645905905016157536392604257102591830824501803683344592010688232302836277603928864745436110"), - (("777", -204), "2.787471972953270789531596912111625325974789615194854615319795902911796043681078997362635440358922503E+103"), - (("7", -600), "2.645751311064590590501615753639260425710259183082450180368334459201068823230283627760392886474543611E+300"), - (("2", -900), "1.414213562373095048801688724209698078569671875376948073176679737990732478462107038850387534327641573E+450"), - (("7", -999), "8.366600265340755479781720257851874893928153692986721998111915430804187725943170098308147119649515362E+499"), - (("74908163946345982392040522594123773796", -999), "2.736935584670307552030924971360722787091742391079630976117950955395149091570790266754718322365663909E+518"), - (("20", -1024), "4.472135954999579392818347337462552470881236719223051448541794490821041851275609798828828816757564550E512"), - (("3", 1025), "5.477225575051661134569697828008021339527446949979832542268944497324932771227227338008584361638706258E-513"), + // sqrt(1.21) = 1.1, [Ceiling, Up] should round up + ("1.21", "2", "1", "1", "1", "1", "1", "2"), + // sqrt(2.25) = 1.5, [Ceiling, HalfEven, HalfUp, Up] should round up + ("2.25", "2", "1", "1", "1", "2", "2", "2"), + // sqrt(6.25) = 2.5, [Ceiling, HalfUp, Up] should round up + ("6.25", "3", "2", "2", "2", "2", "3", "3"), + // sqrt(8.41) = 2.9, [Ceiling, HalfDown, HalfEven, HalfUp, Up] should round up + ("8.41", "3", "2", "2", "3", "3", "3", "3"), ]; - for &((s, scale), e) in vals.iter() { - let expected = BigDecimal::from_str(e).unwrap(); - - let sqrt = BigDecimal::new(BigInt::from_str(s).unwrap(), scale).sqrt().unwrap(); - assert_eq!(sqrt, expected); + for &(val, ceiling, down, floor, half_down, half_even, half_up, up) in vals.iter() { + let val = BigDecimal::from_str(val).unwrap(); + let ceiling = BigDecimal::from_str(ceiling).unwrap(); + let down = BigDecimal::from_str(down).unwrap(); + let floor = BigDecimal::from_str(floor).unwrap(); + let half_down = BigDecimal::from_str(half_down).unwrap(); + let half_even = BigDecimal::from_str(half_even).unwrap(); + let half_up = BigDecimal::from_str(half_up).unwrap(); + let up = BigDecimal::from_str(up).unwrap(); + let ctx = Context::default().with_prec(1).unwrap(); + assert_eq!(val.sqrt_with_context(&ctx.with_rounding_mode(RoundingMode::Ceiling)).unwrap(), ceiling); + assert_eq!(val.sqrt_with_context(&ctx.with_rounding_mode(RoundingMode::Down)).unwrap(), down); + assert_eq!(val.sqrt_with_context(&ctx.with_rounding_mode(RoundingMode::Floor)).unwrap(), floor); + assert_eq!(val.sqrt_with_context(&ctx.with_rounding_mode(RoundingMode::HalfDown)).unwrap(), half_down); + assert_eq!(val.sqrt_with_context(&ctx.with_rounding_mode(RoundingMode::HalfEven)).unwrap(), half_even); + assert_eq!(val.sqrt_with_context(&ctx.with_rounding_mode(RoundingMode::HalfUp)).unwrap(), half_up); + assert_eq!(val.sqrt_with_context(&ctx.with_rounding_mode(RoundingMode::Up)).unwrap(), up); } } - #[test] - fn case_sqrt_3242053850483855em13() { - let d: BigDecimal = "324.2053850483855".parse().unwrap(); - - let digitref = d.to_ref(); - let (_, scale, uint) = digitref.as_parts(); - let ctx = Context::default() - .with_prec(11).unwrap() - .with_rounding_mode(RoundingMode::Down); - - let s = impl_sqrt(uint, scale, &ctx); - let expected: BigDecimal = "18.005704236".parse().unwrap(); - assert_eq!(s, expected); - - let ctx = Context::default() - .with_prec(31).unwrap() - .with_rounding_mode(RoundingMode::Up); - - let s = impl_sqrt(uint, scale, &ctx); - let expected: BigDecimal = "18.00570423639090823994825477228".parse().unwrap(); - assert_eq!(s, expected); - } - - #[test] - fn case_sqrt_5085019992340351em25() { - let d: BigDecimal = "5.085019992340351e-10".parse().unwrap(); - - let digitref = d.to_ref(); - let (_, scale, uint) = digitref.as_parts(); - let ctx = Context::default() - .with_prec(25).unwrap() - .with_rounding_mode(RoundingMode::Down); - - let s = impl_sqrt(uint, scale, &ctx); - let expected: BigDecimal = "0.00002254998889653906459324292".parse().unwrap(); - assert_eq!(s, expected); - } - #[cfg(property_tests)] mod prop { use super::*; diff --git a/src/context.rs b/src/context.rs index 32beda9..388d742 100644 --- a/src/context.rs +++ b/src/context.rs @@ -4,11 +4,11 @@ use crate::*; use stdlib::num::NonZeroU64; +use arithmetic::store_carry; + // const DEFAULT_PRECISION: u64 = ${RUST_BIGDECIMAL_DEFAULT_PRECISION} or 100; include!(concat!(env!("OUT_DIR"), "/default_precision.rs")); -// const DEFAULT_ROUNDING_MODE: RoundingMode = ${RUST_BIGDECIMAL_DEFAULT_ROUNDING_MODE} or HalfUp; -include!(concat!(env!("OUT_DIR"), "/default_rounding_mode.rs")); /// Mathematical Context @@ -83,21 +83,34 @@ impl Context { self.rounding } + /// Round decimal to precision in this context, using rounding-mode + pub fn round_decimal(&self, n: BigDecimal) -> BigDecimal { + n.with_precision_round(self.precision(), self.rounding_mode()) + } + + /// Round decimal to precision in this context, using rounding-mode + pub fn round_decimal_ref<'a, D: Into>>(&self, n: D) -> BigDecimal { + let d = n.into().to_owned(); + d.with_precision_round(self.precision(), self.rounding_mode()) + } + /// Round digits x and y with the rounding mode + #[allow(dead_code)] pub(crate) fn round_pair(&self, sign: Sign, x: u8, y: u8, trailing_zeros: bool) -> u8 { self.rounding.round_pair(sign, (x, y), trailing_zeros) } /// Round digits x and y with the rounding mode #[allow(dead_code)] - pub(crate) fn round_pair_with_carry(&self, sign: Sign, x: u8, y: u8, trailing_zeros: bool, carry: &mut u8) -> u8 { - let r = self.round_pair(sign, x, y, trailing_zeros); - if r == 10 { - *carry = 1; - 0 - } else { - r - } + pub(crate) fn round_pair_with_carry( + &self, + sign: Sign, + x: u8, + y: u8, + trailing_zeros: bool, + carry: &mut u8, + ) -> u8 { + self.rounding.round_pair_with_carry(sign, (x, y), trailing_zeros, carry) } } @@ -105,7 +118,7 @@ impl stdlib::default::Default for Context { fn default() -> Self { Self { precision: NonZeroU64::new(DEFAULT_PRECISION).unwrap(), - rounding: DEFAULT_ROUNDING_MODE, + rounding: RoundingMode::default(), } } } @@ -170,4 +183,48 @@ mod test_context { let sum = ctx.with_prec(27).unwrap().with_rounding_mode(RoundingMode::Up).add_refs(&a, neg_b); assert_eq!(sum, "209682.134972197165534775309".parse().unwrap()); } + + mod round_decimal_ref { + use super::*; + + #[test] + fn case_bigint_1234567_prec3() { + let ctx = Context::default().with_prec(3).unwrap(); + let i = BigInt::from(1234567); + let d = ctx.round_decimal_ref(&i); + assert_eq!(d.int_val, 123.into()); + assert_eq!(d.scale, -4); + } + + #[test] + fn case_bigint_1234500_prec4_halfup() { + let ctx = Context::default() + .with_prec(4).unwrap() + .with_rounding_mode(RoundingMode::HalfUp); + let i = BigInt::from(1234500); + let d = ctx.round_decimal_ref(&i); + assert_eq!(d.int_val, 1235.into()); + assert_eq!(d.scale, -3); + } + + #[test] + fn case_bigint_1234500_prec4_halfeven() { + let ctx = Context::default() + .with_prec(4).unwrap() + .with_rounding_mode(RoundingMode::HalfEven); + let i = BigInt::from(1234500); + let d = ctx.round_decimal_ref(&i); + assert_eq!(d.int_val, 1234.into()); + assert_eq!(d.scale, -3); + } + + #[test] + fn case_bigint_1234567_prec10() { + let ctx = Context::default().with_prec(10).unwrap(); + let i = BigInt::from(1234567); + let d = ctx.round_decimal_ref(&i); + assert_eq!(d.int_val, 1234567000.into()); + assert_eq!(d.scale, 3); + } + } } diff --git a/src/impl_fmt.rs b/src/impl_fmt.rs index 058f442..58c7370 100644 --- a/src/impl_fmt.rs +++ b/src/impl_fmt.rs @@ -2,7 +2,9 @@ //! use crate::*; +use rounding::{NonDigitRoundingData, InsigData}; use stdlib::fmt::Write; +use stdlib::num::NonZeroUsize; // const EXPONENTIAL_FORMAT_LEADING_ZERO_THRESHOLD: usize = ${RUST_BIGDECIMAL_FMT_EXPONENTIAL_LOWER_THRESHOLD} or 5; @@ -125,28 +127,33 @@ fn format_full_scale( use stdlib::cmp::Ordering::*; let mut digits = abs_int.into_bytes(); - let mut exp = (this.scale as i128).neg(); + let mut exp = 0; let non_negative = matches!(this.sign, Sign::Plus | Sign::NoSign); debug_assert_ne!(digits.len(), 0); + let rounder = NonDigitRoundingData::default_with_sign(this.sign); + if this.scale <= 0 { - // formatting an integer value (add trailing zeros to the right) + exp = (this.scale as i128).neg(); + // format an integer value by adding trailing zeros to the right zero_right_pad_integer_ascii_digits(&mut digits, &mut exp, f.precision()); } else { let scale = this.scale as u64; - // no-precision behaves the same as precision matching scale (i.e. no padding or rounding) - let prec = f.precision().and_then(|prec| prec.to_u64()).unwrap_or(scale); + + // std::fmt 'precision' has same meaning as bigdecimal 'scale' + // + // interpret 'no-precision' to mean the same as matching the scale + // of the deicmal (i.e. no padding or rounding) + let target_scale = f.precision().and_then(|prec| prec.to_u64()).unwrap_or(scale); if scale < digits.len() as u64 { // format both integer and fractional digits (always 'trim' to precision) - trim_ascii_digits(&mut digits, scale, prec, &mut exp, this.sign); + format_ascii_digits_with_integer_and_fraction(&mut digits, scale, target_scale, rounder); } else { // format only fractional digits - shift_or_trim_fractional_digits(&mut digits, scale, prec, &mut exp, this.sign); + format_ascii_digits_no_integer(&mut digits, scale, target_scale, rounder); } - // never print exp when in this branch - exp = 0; } // move digits back into String form @@ -168,126 +175,221 @@ fn format_full_scale( fn zero_right_pad_integer_ascii_digits( digits: &mut Vec, exp: &mut i128, - precision: Option, + // number of zeros after the decimal point + target_scale: Option, ) { debug_assert!(*exp >= 0); + debug_assert_ne!(digits.len(), 0); - let trailing_zero_count = match exp.to_usize() { + let integer_zero_count = match exp.to_usize() { Some(n) => n, None => { return; } }; - let total_additional_zeros = trailing_zero_count.saturating_add(precision.unwrap_or(0)); + + // did not explicitly request precision, so we'll only + // implicitly right-pad if less than this threshold. + if matches!(target_scale, None) && integer_zero_count > 20 { + // no padding + return; + } + + let fraction_zero_char_count; + let decimal_place_idx; + + if let Some(frac_zero_count) = target_scale.and_then(NonZeroUsize::new) { + // add one char for '.' if target_scale is not zero + fraction_zero_char_count = frac_zero_count.get() + 1; + // indicate we'll need to add a decimal point + decimal_place_idx = Some(digits.len() + integer_zero_count); + } else { + fraction_zero_char_count = 0; + decimal_place_idx = None; + } + + let total_additional_zeros = integer_zero_count.saturating_add(fraction_zero_char_count); + + // no padding if out of bounds if total_additional_zeros > FMT_MAX_INTEGER_PADDING { return; } - // requested 'prec' digits of precision after decimal point - match precision { - None if trailing_zero_count > 20 => { - } - None | Some(0) => { - digits.resize(digits.len() + trailing_zero_count, b'0'); - *exp = 0; - } - Some(prec) => { - digits.resize(digits.len() + trailing_zero_count, b'0'); - digits.push(b'.'); - digits.resize(digits.len() + prec, b'0'); - *exp = 0; - } + digits.resize(digits.len() + total_additional_zeros, b'0'); + if let Some(decimal_place_idx) = decimal_place_idx { + digits[decimal_place_idx] = b'.'; } + + // set exp to zero so it won't be written in `format_full_scale` + *exp = 0; } -/// Fill zeros into utf-8 digits -fn trim_ascii_digits( - digits: &mut Vec, + +/// Insert decimal point into digits_ascii_be, rounding or padding with zeros when necessary +/// +/// (digits_ascii_be, scale) represents a decimal with both integer and fractional digits. +/// +fn format_ascii_digits_with_integer_and_fraction( + digits_ascii_be: &mut Vec, scale: u64, - prec: u64, - exp: &mut i128, - sign: Sign, + target_scale: u64, + rounder: NonDigitRoundingData, ) { - debug_assert!(scale < digits.len() as u64); - // there are both integer and fractional digits - let integer_digit_count = (digits.len() as u64 - scale) - .to_usize() - .expect("Number of digits exceeds maximum usize"); - - if prec < scale { - let prec = prec.to_usize() - .expect("Precision exceeds maximum usize"); - apply_rounding_to_ascii_digits( - digits, exp, integer_digit_count + prec, sign - ); + debug_assert!(scale < digits_ascii_be.len() as u64, "No integer digits"); + let mut digit_scale = scale; + + // decimal has more fractional digits than requested: round (trimming insignificant digits) + if target_scale < scale { + let digit_count_to_remove = (scale - target_scale) + .to_usize() + .expect("Precision exceeds maximum usize"); + + let rounding_idx = NonZeroUsize::new(digits_ascii_be.len() - digit_count_to_remove) + .expect("decimal must have integer digits"); + + // round and trim the digits at the 'rounding index' + let scale_diff = round_ascii_digits(digits_ascii_be, rounding_idx, rounder); + + match scale_diff.checked_sub(digit_scale as usize) { + None | Some(0) => { + digit_scale -= scale_diff as u64; + } + Some(zeros_to_add) => { + debug_assert_eq!(zeros_to_add, 1); + digits_ascii_be.resize(digits_ascii_be.len() + zeros_to_add, b'0'); + digit_scale = 0; + } + } } - if prec != 0 { - digits.insert(integer_digit_count, b'.'); + // ignore the decimal point if the target scale is zero + if target_scale != 0 { + // there are both integer and fractional digits + let integer_digit_count = (digits_ascii_be.len() as u64 - digit_scale) + .to_usize() + .expect("Number of digits exceeds maximum usize"); + + + digits_ascii_be.insert(integer_digit_count, b'.'); } - if scale < prec { - let trailing_zero_count = (prec - scale) + if digit_scale < target_scale { + let trailing_zero_count = (target_scale - digit_scale) .to_usize() .expect("Too Big"); // precision required beyond scale - digits.resize(digits.len() + trailing_zero_count, b'0'); + digits_ascii_be.resize(digits_ascii_be.len() + trailing_zero_count, b'0'); + digit_scale += trailing_zero_count as u64; } + + debug_assert_eq!(digit_scale, target_scale); } -fn shift_or_trim_fractional_digits( - digits: &mut Vec, +/// Insert decimal point into digits_ascii_be, rounding or padding with zeros when necessary +/// +/// (digits_ascii_be, scale) represents a decimal with only fractional digits. +/// +fn format_ascii_digits_no_integer( + digits_ascii_be: &mut Vec, scale: u64, - prec: u64, - exp: &mut i128, - sign: Sign, + target_scale: u64, + rounder: NonDigitRoundingData, ) { - debug_assert!(scale >= digits.len() as u64); - // there are no integer digits - let leading_zeros = scale - digits.len() as u64; + use stdlib::cmp::Ordering::*; - match prec.checked_sub(leading_zeros) { - None => { - digits.clear(); - digits.push(b'0'); - if prec > 0 { - digits.push(b'.'); - digits.resize(2 + prec as usize, b'0'); + debug_assert!(scale >= digits_ascii_be.len() as u64); + let leading_zeros = scale - digits_ascii_be.len() as u64; + + match arithmetic::diff(target_scale, leading_zeros) { + // handle rounding point before the start of digits + (Less, intermediate_zeros) | (Equal, intermediate_zeros) => { + // get insignificant digit + let (insig_digit, trailing_digits) = if intermediate_zeros > 0 { + (0, digits_ascii_be.as_slice()) + } else { + (digits_ascii_be[0] - b'0', &digits_ascii_be[1..]) + }; + + let insig_data = InsigData::from_digit_and_lazy_trailing_zeros( + rounder, insig_digit, || trailing_digits.iter().all(|&d| d == b'0') + ); + + let rounded_digit = insig_data.round_digit(0); + debug_assert_ne!(rounded_digit, 10); + + digits_ascii_be.clear(); + + if target_scale > 0 { + digits_ascii_be.resize(target_scale as usize + 1, b'0'); } - } - Some(0) => { - // precision is at the decimal digit boundary, round one value - let insig_digit = digits[0] - b'0'; - let trailing_zeros = digits[1..].iter().all(|&d| d == b'0'); - let rounded_value = Context::default().round_pair(sign, 0, insig_digit, trailing_zeros); - - digits.clear(); - if leading_zeros != 0 { - digits.push(b'0'); - digits.push(b'.'); - digits.resize(1 + leading_zeros as usize, b'0'); + + digits_ascii_be.push(rounded_digit + b'0'); + + if target_scale > 0 { + digits_ascii_be[1] = b'.'; } - digits.push(rounded_value + b'0'); } - Some(digit_prec) => { - let digit_prec = digit_prec as usize; - let leading_zeros = leading_zeros - .to_usize() - .expect("Number of leading zeros exceeds max usize"); - let trailing_zeros = digit_prec.saturating_sub(digits.len()); - if digit_prec < digits.len() { - apply_rounding_to_ascii_digits(digits, exp, digit_prec, sign); + (Greater, sig_digit_count) => { + let significant_digit_count = sig_digit_count + .to_usize() + .and_then(NonZeroUsize::new) + .expect("Request overflow in sig_digit_count"); + + let mut digit_scale = scale; + + // if 'digits_ascii_be' has more digits than requested, round + if significant_digit_count.get() < digits_ascii_be.len() { + let removed_digit_count = round_ascii_digits( + digits_ascii_be, significant_digit_count, rounder + ); + + digit_scale -= removed_digit_count as u64; + } + + // number of zeros to keep after the significant digits + let trailing_zeros = target_scale - digit_scale; + + // expected length is target scale (number of digits after decimal point) + "0." + let dest_len = target_scale as usize + 2; + + // number of significant digits is whatever is left in the digit vector + let sig_digit_count = digits_ascii_be.len(); + + // index where to store significant digits + let sig_digit_idx = dest_len - trailing_zeros as usize - sig_digit_count; + + // fill output with zeros + digits_ascii_be.resize(dest_len, b'0'); + + // very likely case where there are digits after the decimal point + if digit_scale != 0 { + // copy significant digits to their index location + digits_ascii_be.copy_within(..sig_digit_count, sig_digit_idx); + + // clear copied values + fill_slice(&mut digits_ascii_be[..sig_digit_count.min(sig_digit_idx)], b'0'); + } else { + debug_assert_eq!(sig_digit_count, 1); } - digits.extend_from_slice(b"0."); - digits.resize(digits.len() + leading_zeros, b'0'); - digits.rotate_right(leading_zeros + 2); - // add any extra trailing zeros - digits.resize(digits.len() + trailing_zeros, b'0'); + // add decimal point + digits_ascii_be[1] = b'.'; } } } +#[cfg(rust_1_50)] +fn fill_slice(v: &mut [T], value: T) { + v.fill(value); +} + +#[cfg(not(rust_1_50))] +fn fill_slice(v: &mut [T], value: T) { + for i in v.iter_mut() { + *i = value.clone(); + } +} + /// Format integer as {int}e+{exp} /// @@ -341,21 +443,16 @@ fn format_exponential_bigendian_ascii_digits( if let Some(prec) = f.precision() { // 'prec' is number of digits after the decimal point let total_prec = prec + 1; - let digit_count = digits.len(); - match total_prec.cmp(&digit_count) { - Ordering::Equal => { - // digit count is one more than precision - do nothing - } - Ordering::Less => { - // round to smaller precision - apply_rounding_to_ascii_digits(&mut digits, &mut exp, total_prec, sign); - } - Ordering::Greater => { - // increase number of zeros to add to end of digits - extra_trailing_zero_count = total_prec - digit_count; - } + if total_prec < digits.len() { + // round to smaller precision + let rounder = NonDigitRoundingData::default_with_sign(sign); + let target_scale = NonZeroUsize::new(total_prec).unwrap(); + let delta_exp = round_ascii_digits(&mut digits, target_scale, rounder); + exp += delta_exp as i128; } + + extra_trailing_zero_count = total_prec - digits.len(); } let needs_decimal_point = digits.len() > 1 || extra_trailing_zero_count > 0; @@ -423,6 +520,77 @@ fn format_exponential_bigendian_ascii_digits( } +/// Round big-endian ascii digits to given significant digit count, +/// updating the scale appropriately +/// +/// Returns the number of digits removed; this should be treated as the +/// change in the decimal's scale, and should be subtracted from the scale +/// when appropriate. +/// +fn round_ascii_digits( + // bigendian ascii digits + digits_ascii_be: &mut Vec, + // number of significant digits to keep + significant_digit_count: NonZeroUsize, + // how to round + rounder: NonDigitRoundingData, +) -> usize { + debug_assert!(significant_digit_count.get() < digits_ascii_be.len()); + + let (sig_digits, insig_digits) = digits_ascii_be.split_at(significant_digit_count.get()); + let (&insig_digit, trailing_digits) = insig_digits.split_first().unwrap_or((&b'0', &[])); + + let insig_data = InsigData::from_digit_and_lazy_trailing_zeros( + rounder, insig_digit - b'0', || trailing_digits.iter().all(|&d| d == b'0') + ); + + let rounding_digit_pos = significant_digit_count.get() - 1; + let sig_digit = sig_digits[rounding_digit_pos] - b'0'; + let rounded_digit = insig_data.round_digit(sig_digit); + + // record how many digits to remove (changes the 'scale') + let mut removed_digit_count = insig_digits.len(); + + // discard insignificant digits (and rounding digit) + digits_ascii_be.truncate(rounding_digit_pos); + + if rounded_digit < 10 { + // simple case: no carrying/overflow, push rounded digit + digits_ascii_be.push(rounded_digit + b'0'); + return removed_digit_count; + } + + // handle carrying the 1 + debug_assert!(rounded_digit == 10); + + // carry the one past trailing 9's: replace them with zeros + let next_non_nine_rev_pos = digits_ascii_be.iter().rev().position(|&d| d != b'9'); + match next_non_nine_rev_pos { + // number of nines to remove + Some(backwards_nine_count) => { + let digits_to_trim = backwards_nine_count + 1; + let idx = digits_ascii_be.len() - digits_to_trim; + // increment least significant non-nine zero + digits_ascii_be[idx] += 1; + // remove trailing nines + digits_ascii_be.truncate(idx + 1); + // count truncation + removed_digit_count += digits_to_trim; + } + // all nines! overflow to 1.000 + None => { + digits_ascii_be.clear(); + digits_ascii_be.push(b'1'); + // count digit clearning + removed_digit_count += significant_digit_count.get(); + } + } + + // return the number of removed digits + return removed_digit_count; +} + + #[inline(never)] pub(crate) fn write_scientific_notation(n: &BigDecimal, w: &mut W) -> fmt::Result { if n.is_zero() { @@ -492,67 +660,6 @@ pub(crate) fn write_engineering_notation(n: &BigDecimal, out: &mut W) } -/// Round big-endian digits in ascii -fn apply_rounding_to_ascii_digits( - ascii_digits: &mut Vec, - exp: &mut i128, - prec: usize, - sign: Sign -) { - if ascii_digits.len() < prec { - return; - } - - // shift exp to align with new length of digits - *exp += (ascii_digits.len() - prec) as i128; - - // true if all ascii_digits after precision are zeros - let trailing_zeros = ascii_digits[prec + 1..].iter().all(|&d| d == b'0'); - - let sig_digit = ascii_digits[prec - 1] - b'0'; - let insig_digit = ascii_digits[prec] - b'0'; - let rounded_digit = Context::default().round_pair(sign, sig_digit, insig_digit, trailing_zeros); - - // remove insignificant digits - ascii_digits.truncate(prec - 1); - - // push rounded value - if rounded_digit < 10 { - ascii_digits.push(rounded_digit + b'0'); - return - } - - debug_assert_eq!(rounded_digit, 10); - - // push zero and carry-the-one - ascii_digits.push(b'0'); - - // loop through digits in reverse order (skip the 0 we just pushed) - let digits = ascii_digits.iter_mut().rev().skip(1); - for digit in digits { - if *digit < b'9' { - // we've carried the one as far as it will go - *digit += 1; - return; - } - - debug_assert_eq!(*digit, b'9'); - - // digit was a 9, set to zero and carry the one - // to the next digit - *digit = b'0'; - } - - // at this point all digits have become zero - // just set significant digit to 1 and increase exponent - // - // eg: 9999e2 ~> 0000e2 ~> 1000e3 - // - ascii_digits[0] = b'1'; - *exp += 1; -} - - #[cfg(test)] #[allow(non_snake_case)] mod test { @@ -721,6 +828,32 @@ mod test { impl_case!(fmt_p05d1: "{:+05.7}" => "+123456.0000000"); } + mod dec_d099995 { + use super::*; + + fn test_input() -> BigDecimal { + ".099995".parse().unwrap() + } + + impl_case!(fmt_default: "{}" => "0.099995"); + impl_case!(fmt_d0: "{:.0}" => "0"); + impl_case!(fmt_d1: "{:.1}" => "0.1"); + impl_case!(fmt_d3: "{:.3}" => "0.100"); + } + + mod dec_d9999999 { + use super::*; + + fn test_input() -> BigDecimal { + ".9999999".parse().unwrap() + } + + impl_case!(fmt_default: "{}" => "0.9999999"); + impl_case!(fmt_d0: "{:.0}" => "1"); + impl_case!(fmt_d1: "{:.1}" => "1.0"); + impl_case!(fmt_d3: "{:.3}" => "1.000"); + } + mod dec_9999999 { use super::*; @@ -774,11 +907,31 @@ mod test { impl_case!(fmt_d0: "{:.0}" => "-90037660"); impl_case!(fmt_d3: "{:.3}" => "-90037659.690"); impl_case!(fmt_d4: "{:.4}" => "-90037659.6905"); - impl_case!(fmt_14d4: "{:14.4}" => "-90037659.6905"); + impl_case!(fmt_14d4: "{:14.4}" => "-90037659.6905"); impl_case!(fmt_15d4: "{:15.4}" => " -90037659.6905"); impl_case!(fmt_l17d5: "{:<17.5}" => "-90037659.69050 "); } + mod dec_0d0002394899999500 { + use super::*; + + fn test_input() -> BigDecimal { + "0.0002394899999500".parse().unwrap() + } + + impl_case!(fmt_default: "{}" => "0.0002394899999500"); + impl_case!(fmt_d0: "{:.0}" => "0"); + impl_case!(fmt_d1: "{:.1}" => "0.0"); + impl_case!(fmt_d3: "{:.3}" => "0.000"); + impl_case!(fmt_d4: "{:.4}" => "0.0002"); + impl_case!(fmt_d5: "{:.5}" => "0.00024"); + impl_case!(fmt_d10: "{:.10}" => "0.0002394900"); + impl_case!(fmt_d13: "{:.13}" => "0.0002394900000"); + impl_case!(fmt_d14: "{:.14}" => "0.00023948999995"); + impl_case!(fmt_d20: "{:.20}" => "0.00023948999995000000"); + + } + mod dec_1764031078en13 { use super::*; @@ -848,6 +1001,8 @@ mod test { impl_case!(fmt_default: "{}" => "0.00003102564500"); impl_case!(fmt_d0: "{:.0}" => "0"); + impl_case!(fmt_d1: "{:.1}" => "0.0"); + impl_case!(fmt_d2: "{:.2}" => "0.00"); impl_case!(fmt_d4: "{:.4}" => "0.0000"); impl_case!(fmt_d5: "{:.5}" => "0.00003"); impl_case!(fmt_d10: "{:.10}" => "0.0000310256"); @@ -925,6 +1080,36 @@ mod test { impl_case!(fmt_default: "{}" => "13400476439814628800e+2502"); impl_case!(fmt_d1: "{:.1}" => "13400476439814628800e+2502"); } + + mod dec_d9999 { + use super::*; + + fn test_input() -> BigDecimal { + "0.9999".parse().unwrap() + } + + impl_case!(fmt_default: "{}" => "0.9999"); + impl_case!(fmt_d0: "{:.0}" => "1"); + impl_case!(fmt_d1: "{:.1}" => "1.0"); + impl_case!(fmt_d2: "{:.2}" => "1.00"); + impl_case!(fmt_d3: "{:.3}" => "1.000"); + impl_case!(fmt_d4: "{:.4}" => "0.9999"); + impl_case!(fmt_d6: "{:.6}" => "0.999900"); + } + + mod dec_9d99 { + use super::*; + + fn test_input() -> BigDecimal { + "9.99".parse().unwrap() + } + + impl_case!(fmt_default: "{}" => "9.99"); + impl_case!(fmt_d0: "{:.0}" => "10"); + impl_case!(fmt_d1: "{:.1}" => "10.0"); + impl_case!(fmt_d2: "{:.2}" => "9.99"); + impl_case!(fmt_d3: "{:.3}" => "9.990"); + } } mod fmt_boundaries { diff --git a/src/impl_ops.rs b/src/impl_ops.rs index b7afec1..784a4e0 100644 --- a/src/impl_ops.rs +++ b/src/impl_ops.rs @@ -250,6 +250,7 @@ macro_rules! impl_div_for_primitive { type Output = BigDecimal; #[cfg(rustc_1_70)] // Option::is_some_and + #[allow(clippy::incompatible_msrv)] fn div(self, denom: $t) -> BigDecimal { if denom.is_one() { self diff --git a/src/impl_ops_add.rs b/src/impl_ops_add.rs index f056c95..ac099e4 100644 --- a/src/impl_ops_add.rs +++ b/src/impl_ops_add.rs @@ -46,7 +46,7 @@ impl Add for &'_ BigDecimal { impl<'a, T: Into>> Add for &'_ BigDecimal { type Output = BigDecimal; fn add(self, rhs: T) -> BigDecimal { - arithmetic::addition::add_bigdecimal_refs(self, rhs) + arithmetic::addition::add_bigdecimal_refs(self, rhs, None) } } @@ -72,7 +72,7 @@ impl Add for BigDecimalRef<'_> { impl<'a, T: Into>> Add for BigDecimalRef<'_> { type Output = BigDecimal; fn add(self, rhs: T) -> BigDecimal { - arithmetic::addition::add_bigdecimal_refs(self, rhs) + arithmetic::addition::add_bigdecimal_refs(self, rhs, None) } } diff --git a/src/impl_ops_mul.rs b/src/impl_ops_mul.rs index 181ce9d..a1cbc43 100644 --- a/src/impl_ops_mul.rs +++ b/src/impl_ops_mul.rs @@ -213,72 +213,70 @@ impl MulAssign for BigDecimal { #[cfg(test)] #[allow(non_snake_case)] mod bigdecimal_tests { - use crate::{stdlib, BigDecimal, ToString, FromStr, TryFrom}; + use super::*; use num_traits::{ToPrimitive, FromPrimitive, Signed, Zero, One}; use num_bigint; use paste::paste; - /// Test multiplication of two bigdecimals - #[test] - fn test_mul() { - - let vals = vec![ - ("2", "1", "2"), - ("12.34", "1.234", "15.22756"), - ("2e1", "1", "20"), - ("3", ".333333", "0.999999"), - ("2389472934723", "209481029831", "500549251119075878721813"), - ("1e-450", "1e500", ".1e51"), - ("-995052931372975485719.533153137", "4.523087321", "-4500711297616988541501.836966993116075977"), - ("995052931372975485719.533153137", "-4.523087321", "-4500711297616988541501.836966993116075977"), - ("-8.37664968", "-1.9086963714056968482094712882596748", "15.988480848752691653730876239769592670324064"), - ("-8.37664968", "0", "0"), - ]; - - for &(x, y, z) in vals.iter() { - - let mut a = BigDecimal::from_str(x).unwrap(); - let b = BigDecimal::from_str(y).unwrap(); - let c = BigDecimal::from_str(z).unwrap(); - - assert_eq!(a.clone() * b.clone(), c); - assert_eq!(a.clone() * &b, c); - assert_eq!(&a * b.clone(), c); - assert_eq!(&a * &b, c); - - a *= b; - assert_eq!(a, c); - } + macro_rules! impl_test { + ($name:ident; $a:literal * $b:literal => $expected:literal) => { + #[test] + fn $name() { + let mut a: BigDecimal = $a.parse().unwrap(); + let b: BigDecimal = $b.parse().unwrap(); + let expected = $expected.parse().unwrap(); + + let prod = a.clone() * b.clone(); + assert_eq!(prod, expected); + assert_eq!(prod.scale, expected.scale); + + assert_eq!(a.clone() * &b, expected); + assert_eq!(&a * b.clone(), expected); + assert_eq!(&a * &b, expected); + + a *= b; + assert_eq!(a, expected); + assert_eq!(a.scale, expected.scale); + } + }; + ($name:ident; $bigt:ty; $a:literal * $b:literal => $expected:literal) => { + #[test] + fn $name() { + let a: BigDecimal = $a.parse().unwrap(); + let b: $bigt = $b.parse().unwrap(); + let c = $expected.parse().unwrap(); + + let prod = a.clone() * b.clone(); + assert_eq!(prod, c); + assert_eq!(prod.scale, c.scale); + + assert_eq!(b.clone() * a.clone(), c); + assert_eq!(a.clone() * &b, c); + assert_eq!(b.clone() * &a, c); + assert_eq!(&a * b.clone(), c); + assert_eq!(&b * a.clone(), c); + assert_eq!(&a * &b, c); + assert_eq!(&b * &a, c); + } + }; } - /// Test multiplication between big decimal and big integer - #[test] - fn test_mul_bigint() { - let vals = vec![ - ("2", "1", "2"), - ("8.561", "10", "85.61"), - ("1.0000", "638655273892892437", "638655273892892437"), - ("10000", "638655273892892437", "6386552738928924370000"), - (".0005", "368408638655273892892437473", "184204319327636946446218.7365"), - ("9e-1", "368408638655273892892437473", "331567774789746503603193725.7"), - ("-1.175470587012343730098", "577575785", "-678923347.038065234601180476930"), - ("-1.175470587012343730098", "-76527768352678", "89956140788267.069799533723307502444"), - ("-1.175470587012343730098", "0", "0"), - ]; - - for &(x, y, z) in vals.iter() { - let a = BigDecimal::from_str(x).unwrap(); - let b = num_bigint::BigInt::from_str(y).unwrap(); - let c = BigDecimal::from_str(z).unwrap(); - - assert_eq!(a.clone() * b.clone(), c); - assert_eq!(b.clone() * a.clone(), c); - assert_eq!(a.clone() * &b, c); - assert_eq!(b.clone() * &a, c); - assert_eq!(&a * b.clone(), c); - assert_eq!(&b * a.clone(), c); - assert_eq!(&a * &b, c); - assert_eq!(&b * &a, c); - } - } + impl_test!(case_2_1; "2" * "1" => "2"); + impl_test!(case_12d34_1d234; "12.34" * "1.234" => "15.22756"); + impl_test!(case_2e1_1; "2e1" * "1" => "2e1"); + impl_test!(case_3_d333333; "3" * ".333333" => "0.999999"); + impl_test!(case_2389472934723_209481029831; "2389472934723" * "209481029831" => "500549251119075878721813"); + impl_test!(case_1ed450_1e500; "1e-450" * "1e500" => "0.1e51"); + impl_test!(case_n995052931ddd_4d523087321; "-995052931372975485719.533153137" * "4.523087321" => "-4500711297616988541501.836966993116075977"); + impl_test!(case_995052931ddd_n4d523087321; "995052931372975485719.533153137" * "-4.523087321" => "-4500711297616988541501.836966993116075977"); + impl_test!(case_n8d37664968_n4d523087321; "-8.37664968" * "-1.9086963714056968482094712882596748" => "15.988480848752691653730876239769592670324064"); + impl_test!(case_n8d37664968_0; "-8.37664968" * "0" => "0.00000000"); + + impl_test!(case_8d561_10; BigInt; "8.561" * "10" => "85.610"); + + // Test multiplication between big decimal and big integer + impl_test!(case_10000_638655273892892437; BigInt; "10000" * "638655273892892437" => "6386552738928924370000"); + impl_test!(case_1en10_n9056180052657301; BigInt; "1e-10" * "-9056180052657301" => "-905618.0052657301"); + impl_test!(case_n9en1_n368408638655273892892437473; BigInt; "-9e-1" * "-368408638655273892892437473" => "331567774789746503603193725.7"); + impl_test!(case_n1d175470587012343730098_577575785; BigInt; "-1.175470587012343730098" * "577575785" => "-678923347.038065234601180476930"); } diff --git a/src/lib.rs b/src/lib.rs index e2a3568..f92c6ed 100644 --- a/src/lib.rs +++ b/src/lib.rs @@ -91,7 +91,7 @@ use self::stdlib::Vec; use num_bigint::{BigInt, BigUint, ParseBigIntError, Sign}; use num_integer::Integer as IntegerTrait; -pub use num_traits::{FromPrimitive, Num, One, Signed, ToPrimitive, Zero}; +pub use num_traits::{FromPrimitive, Num, One, Pow, Signed, ToPrimitive, Zero}; use stdlib::f64::consts::LOG2_10; @@ -219,14 +219,28 @@ fn exp2(x: f64) -> f64 { impl BigDecimal { /// Creates and initializes a `BigDecimal`. /// + /// The more explicit method `from_bigint` should be preferred, as new + /// may change in the future. + /// #[inline] pub fn new(digits: BigInt, scale: i64) -> BigDecimal { + BigDecimal::from_bigint(digits, scale) + } + + /// Construct BigDecimal from BigInt and a scale + pub fn from_bigint(digits: BigInt, scale: i64) -> BigDecimal { BigDecimal { int_val: digits, scale: scale, } } + /// Construct positive BigDecimal from BigUint and a scale + pub fn from_biguint(digits: BigUint, scale: i64) -> BigDecimal { + let n = BigInt::from_biguint(Sign::Plus, digits); + BigDecimal::from_bigint(n, scale) + } + /// Make a BigDecimalRef of this value pub fn to_ref(&self) -> BigDecimalRef<'_> { // search for "From<&'a BigDecimal> for BigDecimalRef<'a>" @@ -428,6 +442,13 @@ impl BigDecimal { } } + /// set scale only if new_scale is greater than current + pub(crate) fn extend_scale_to(&mut self, new_scale: i64) { + if new_scale > self.scale { + self.set_scale(new_scale) + } + } + /// Take and return bigdecimal with the given sign /// /// The Sign value `NoSign` is ignored: only use Plus & Minus @@ -490,6 +511,7 @@ impl BigDecimal { /// Return this BigDecimal with the given precision, rounding if needed #[cfg(rustc_1_46)] // Option::zip + #[allow(clippy::incompatible_msrv)] pub fn with_precision_round(&self, prec: stdlib::num::NonZeroU64, round: RoundingMode) -> BigDecimal { let digit_count = self.digits(); let new_prec = prec.get().to_i64(); @@ -740,11 +762,7 @@ impl BigDecimal { return self.clone(); } - let uint = self.int_val.magnitude(); - let result = arithmetic::cbrt::impl_cbrt_uint_scale(uint, self.scale, ctx); - - // always copy sign - result.take_with_sign(self.sign()) + arithmetic::cbrt::impl_cbrt_int_scale(&self.int_val, self.scale, ctx) } /// Compute the reciprical of the number: x-1 @@ -766,96 +784,15 @@ impl BigDecimal { result.take_with_sign(self.sign()) } - /// Return number rounded to round_digits precision after the decimal point + /// Return given number rounded to 'round_digits' precision after the + /// decimal point, using default rounding mode + /// + /// Default rounding mode is `HalfEven`, but can be configured at compile-time + /// by the environment variable: `RUST_BIGDECIMAL_DEFAULT_ROUNDING_MODE` + /// (or by patching _build.rs_ ) + /// pub fn round(&self, round_digits: i64) -> BigDecimal { - // we have fewer digits than we need, no rounding - if round_digits >= self.scale { - return self.with_scale(round_digits); - } - - let (sign, double_digits) = self.int_val.to_radix_le(100); - - let last_is_double_digit = *double_digits.last().unwrap() >= 10; - let digit_count = (double_digits.len() - 1) * 2 + 1 + last_is_double_digit as usize; - - // relevant digit positions: each "pos" is position of 10^{pos} - let least_significant_pos = -self.scale; - let most_sig_digit_pos = digit_count as i64 + least_significant_pos - 1; - let rounding_pos = -round_digits; - - // digits are too small, round to zero - if rounding_pos > most_sig_digit_pos + 1 { - return BigDecimal::zero(); - } - - // highest digit is next to rounding point - if rounding_pos == most_sig_digit_pos + 1 { - let (&last_double_digit, remaining) = double_digits.split_last().unwrap(); - - let mut trailing_zeros = remaining.iter().all(|&d| d == 0); - - let last_digit = if last_is_double_digit { - let (high, low) = num_integer::div_rem(last_double_digit, 10); - trailing_zeros &= low == 0; - high - } else { - last_double_digit - }; - - if last_digit > 5 || (last_digit == 5 && !trailing_zeros) { - return BigDecimal::new(BigInt::one(), round_digits); - } - - return BigDecimal::zero(); - } - - let double_digits_to_remove = self.scale - round_digits; - debug_assert!(double_digits_to_remove > 0); - - let (rounding_idx, rounding_offset) = num_integer::div_rem(double_digits_to_remove as usize, 2); - debug_assert!(rounding_idx <= double_digits.len()); - - let (low_digits, high_digits) = double_digits.as_slice().split_at(rounding_idx); - debug_assert!(high_digits.len() > 0); - - let mut unrounded_uint = num_bigint::BigUint::from_radix_le(high_digits, 100).unwrap(); - - let rounded_uint; - if rounding_offset == 0 { - let high_digit = high_digits[0] % 10; - let (&top, rest) = low_digits.split_last().unwrap_or((&0u8, &[])); - let (low_digit, lower_digit) = num_integer::div_rem(top, 10); - let trailing_zeros = lower_digit == 0 && rest.iter().all(|&d| d == 0); - - let rounding = if low_digit < 5 { - 0 - } else if low_digit > 5 || !trailing_zeros { - 1 - } else { - high_digit % 2 - }; - - rounded_uint = unrounded_uint + rounding; - } else { - let (high_digit, low_digit) = num_integer::div_rem(high_digits[0], 10); - - let trailing_zeros = low_digits.iter().all(|&d| d == 0); - - let rounding = if low_digit < 5 { - 0 - } else if low_digit > 5 || !trailing_zeros { - 1 - } else { - high_digit % 2 - }; - - // shift unrounded_uint down, - unrounded_uint /= num_bigint::BigUint::from_u8(10).unwrap(); - rounded_uint = unrounded_uint + rounding; - } - - let rounded_int = num_bigint::BigInt::from_biguint(sign, rounded_uint); - BigDecimal::new(rounded_int, round_digits) + self.with_scale_round(round_digits, Context::default().rounding_mode()) } /// Return true if this number has zero fractional part (is equal @@ -1241,10 +1178,22 @@ impl BigDecimalRef<'_> { pub fn to_owned_with_scale(&self, scale: i64) -> BigDecimal { use stdlib::cmp::Ordering::*; - let digits = match scale.cmp(&self.scale) { - Equal => self.digits.clone(), - Greater => self.digits * ten_to_the_uint((scale - self.scale) as u64), - Less => self.digits / ten_to_the_uint((self.scale - scale) as u64) + let digits = match arithmetic::diff(self.scale, scale) { + (Equal, _) => self.digits.clone(), + (Less, scale_diff) => { + if scale_diff < 20 { + self.digits * ten_to_the_u64(scale_diff as u8) + } else { + self.digits * ten_to_the_uint(scale_diff) + } + } + (Greater, scale_diff) => { + if scale_diff < 20 { + self.digits / ten_to_the_u64(scale_diff as u8) + } else { + self.digits / ten_to_the_uint(scale_diff) + } + } }; BigDecimal { @@ -1269,6 +1218,19 @@ impl BigDecimalRef<'_> { count_decimal_digits_uint(self.digits) } + /// Return the number of trailing zeros in the referenced integer + #[allow(dead_code)] + fn count_trailing_zeroes(&self) -> usize { + if self.digits.is_zero() || self.digits.is_odd() { + return 0; + } + + let digit_pairs = self.digits.to_radix_le(100); + let loc = digit_pairs.iter().position(|&d| d != 0).unwrap_or(0); + + 2 * loc + usize::from(digit_pairs[loc] % 10 == 0) + } + /// Split into components pub(crate) fn as_parts(&self) -> (Sign, i64, &BigUint) { (self.sign, self.scale, self.digits) @@ -1283,6 +1245,14 @@ impl BigDecimalRef<'_> { } } + /// Create BigDecimal from this reference, rounding to precision and + /// with rounding-mode of the given context + /// + /// + pub fn round_with_context(&self, ctx: &Context) -> BigDecimal { + ctx.round_decimal_ref(*self) + } + /// Take square root of this number pub fn sqrt_with_context(&self, ctx: &Context) -> Option { use Sign::*; @@ -1296,6 +1266,26 @@ impl BigDecimalRef<'_> { } } + /// Take square root of absolute-value of the number + pub fn sqrt_abs_with_context(&self, ctx: &Context) -> BigDecimal { + use Sign::*; + + let (_, scale, uint) = self.as_parts(); + arithmetic::sqrt::impl_sqrt(uint, scale, ctx) + } + + /// Take square root, copying sign of the initial decimal + pub fn sqrt_copysign_with_context(&self, ctx: &Context) -> BigDecimal { + use Sign::*; + + let (sign, scale, uint) = self.as_parts(); + let mut result = arithmetic::sqrt::impl_sqrt(uint, scale, ctx); + if sign == Minus { + result.int_val = result.int_val.neg(); + } + result + } + /// Return if the referenced decimal is zero pub fn is_zero(&self) -> bool { self.digits.is_zero() @@ -1330,6 +1320,53 @@ impl<'a> From<&'a BigInt> for BigDecimalRef<'a> { } } + +/// pair i64 'scale' with some other value +#[derive(Clone, Copy)] +struct WithScale { + pub value: T, + pub scale: i64, +} + +impl fmt::Debug for WithScale { + fn fmt(&self, f: &mut fmt::Formatter) -> fmt::Result { + write!(f, "(scale={} {:?})", self.scale, self.value) + } +} + +impl From<(T, i64)> for WithScale { + fn from(pair: (T, i64)) -> Self { + Self { value: pair.0, scale: pair.1 } + } +} + +impl<'a> From> for BigDecimalRef<'a> { + fn from(obj: WithScale<&'a BigInt>) -> Self { + Self { + scale: obj.scale, + sign: obj.value.sign(), + digits: obj.value.magnitude(), + } + } +} + +impl<'a> From> for BigDecimalRef<'a> { + fn from(obj: WithScale<&'a BigUint>) -> Self { + Self { + scale: obj.scale, + sign: Sign::Plus, + digits: obj.value, + } + } +} + +impl WithScale<&T> { + fn is_zero(&self) -> bool { + self.value.is_zero() + } +} + + #[rustfmt::skip] #[cfg(test)] #[allow(non_snake_case)] @@ -1339,6 +1376,33 @@ mod bigdecimal_tests { use num_bigint; use paste::paste; + + mod from_biguint { + use super::*; + use num_bigint::BigUint; + use num_bigint::Sign; + + macro_rules! impl_case { + ($name:ident; $i:literal; $scale:literal) => { + impl_case!($name; $i.into(); $scale; Plus); + }; + ($name:ident; $i:expr; $scale:literal; $sign:ident) => { + #[test] + fn $name() { + let i: BigUint = $i; + let d = BigDecimal::from_biguint(i.clone(), $scale); + assert_eq!(d.int_val.magnitude(), &i); + assert_eq!(d.scale, $scale); + assert_eq!(d.sign(), Sign::$sign); + } + }; + } + + impl_case!(case_0en3; BigUint::zero(); 3; NoSign); + impl_case!(case_30e2; 30u8; -2); + impl_case!(case_7446124798en5; 7446124798u128; 5); + } + #[test] fn test_fractional_digit_count() { // Zero value @@ -1827,6 +1891,11 @@ mod bigdecimal_tests { ("0.1165085714285714285714285714285714285714", 2, "0.12"), ("0.1165085714285714285714285714285714285714", 5, "0.11651"), ("0.1165085714285714285714285714285714285714", 8, "0.11650857"), + ("-1.5", 0, "-2"), + ("-1.2", 0, "-1"), + ("-0.68", 0, "-1"), + ("-0.5", 0, "0"), + ("-0.49", 0, "0"), ]; for &(x, digits, y) in test_cases.iter() { let a = BigDecimal::from_str(x).unwrap(); diff --git a/src/rounding.rs b/src/rounding.rs index 0a39655..172c5bc 100644 --- a/src/rounding.rs +++ b/src/rounding.rs @@ -1,8 +1,13 @@ //! Rounding structures and subroutines +#![allow(dead_code)] -use crate::Sign; +use crate::*; +use crate::arithmetic::{add_carry, store_carry, extend_adding_with_carry}; use stdlib; +// const DEFAULT_ROUNDING_MODE: RoundingMode = ${RUST_BIGDECIMAL_DEFAULT_ROUNDING_MODE} or HalfUp; +include!(concat!(env!("OUT_DIR"), "/default_rounding_mode.rs")); + /// Determines how to calculate the last digit of the number /// /// Default rounding mode is HalfUp @@ -149,6 +154,18 @@ impl RoundingMode { } } + /// Round digits, and if rounded up to 10, store 1 in carry and return zero + pub(crate) fn round_pair_with_carry( + &self, + sign: Sign, + pair: (u8, u8), + trailing_zeros: bool, + carry: &mut u8, + ) -> u8 { + let r = self.round_pair(sign, pair, trailing_zeros); + store_carry(r, carry) + } + /// Round value at particular digit, returning replacement digit /// /// Parameters @@ -195,6 +212,209 @@ impl RoundingMode { // shift rounded value back to position full * splitter } + + /// Hint used to skip calculating trailing_zeros if they don't matter + fn needs_trailing_zeros(&self, insig_digit: u8) -> bool { + use RoundingMode::*; + + // only need trailing zeros if the rounding digit is 0 or 5 + if matches!(self, HalfUp | HalfDown | HalfEven) { + insig_digit == 5 + } else { + insig_digit == 0 + } + } + +} + +/// Return compile-time constant default rounding mode +/// +/// Defined by RUST_BIGDECIMAL_DEFAULT_ROUNDING_MODE at compile time +/// +impl Default for RoundingMode { + fn default() -> Self { + DEFAULT_ROUNDING_MODE + } +} + + +/// All non-digit information required to round digits +/// +/// Just the mode and the sign. +/// +#[derive(Debug,Clone,Copy)] +pub(crate) struct NonDigitRoundingData { + /// Rounding mode + pub mode: RoundingMode, + /// Sign of digits + pub sign: Sign, +} + +impl NonDigitRoundingData { + /// Round pair of digits, storing overflow (10) in the carry + pub fn round_pair(&self, pair: (u8, u8), trailing_zeros: bool) -> u8 { + self.mode.round_pair(self.sign, pair, trailing_zeros) + } + + /// round-pair with carry-digits + pub fn round_pair_with_carry(&self, pair: (u8, u8), trailing_zeros: bool, carry: &mut u8) -> u8 { + self.mode.round_pair_with_carry(self.sign, pair, trailing_zeros, carry) + } + + /// Use sign and default rounding mode + pub fn default_with_sign(sign: Sign) -> Self { + NonDigitRoundingData { sign, mode: RoundingMode::default() } + } +} + + +/// Relevant information about insignificant digits, used for rounding +/// +/// If rounding at indicated point: +/// +/// ```txt +/// aaaaizzzzzzzz +/// ^ +/// ``` +/// +/// 'a' values are significant, 'i' is the insignificant digit, +/// and trailing_zeros is true if all 'z' are 0. +/// +#[derive(Debug,Clone,Copy)] +pub(crate) struct InsigData { + /// highest insignificant digit + pub digit: u8, + + /// true if all digits more insignificant than 'digit' is zero + /// + /// This is only useful if relevant for the rounding mode, it + /// may be 'wrong' in these cases. + pub trailing_zeros: bool, + + /// rounding-mode and sign + pub rounding_data: NonDigitRoundingData +} + +impl InsigData { + /// Build from insig data and lazily calcuated trailing-zeros callable + pub fn from_digit_and_lazy_trailing_zeros( + rounder: NonDigitRoundingData, + insig_digit: u8, + calc_trailing_zeros: impl FnOnce() -> bool + ) -> Self { + Self { + digit: insig_digit, + trailing_zeros: rounder.mode.needs_trailing_zeros(insig_digit) && calc_trailing_zeros(), + rounding_data: rounder, + } + } + + /// Build from slice of insignificant little-endian digits + pub fn from_digit_slice(rounder: NonDigitRoundingData, digits: &[u8]) -> Self { + match digits.split_last() { + Some((&d0, trailing)) => { + Self::from_digit_and_lazy_trailing_zeros( + rounder, d0, || trailing.iter().all(Zero::is_zero) + ) + } + None => { + Self { + digit: 0, + trailing_zeros: true, + rounding_data: rounder, + } + } + } + } + + /// from sum of overlapping digits, (a is longer than b) + pub fn from_overlapping_digits_backward_sum( + rounder: NonDigitRoundingData, + mut a_digits: stdlib::iter::Rev>, + mut b_digits: stdlib::iter::Rev>, + carry: &mut u8, + ) -> Self { + debug_assert!(a_digits.len() >= b_digits.len()); + debug_assert_eq!(carry, &0); + + // most-significant insignificant digit + let insig_digit; + match (a_digits.next(), b_digits.next()) { + (Some(a), Some(b)) => { + // store 'full', initial sum, we will handle carry below + insig_digit = a + b; + } + (Some(d), None) | (None, Some(d)) => { + insig_digit = *d; + } + (None, None) => { + // both digit slices were empty; all zeros + return Self { + digit: 0, + trailing_zeros: true, + rounding_data: rounder, + }; + } + }; + + // find first non-nine value + let mut sum = 9; + while sum == 9 { + let next_a = a_digits.next().unwrap_or(&0); + let next_b = b_digits.next().unwrap_or(&0); + sum = next_a + next_b; + } + + // if previous sum was greater than ten, + // the one would carry through all the 9s + let sum = store_carry(sum, carry); + + // propagate carry to the highest insignificant digit + let insig_digit = add_carry(insig_digit, carry); + + // if the last 'sum' value isn't zero, or if any remaining + // digit is not zero, then it's not trailing zeros + let trailing_zeros = sum == 0 + && rounder.mode.needs_trailing_zeros(insig_digit) + && a_digits.all(Zero::is_zero) + && b_digits.all(Zero::is_zero); + + Self { + digit: insig_digit, + trailing_zeros: trailing_zeros, + rounding_data: rounder, + } + } + + pub fn round_digit(&self, digit: u8) -> u8 { + self.rounding_data.round_pair((digit, self.digit), self.trailing_zeros) + } + + pub fn round_digit_with_carry(&self, digit: u8, carry: &mut u8) -> u8 { + self.rounding_data.round_pair_with_carry((digit, self.digit), self.trailing_zeros, carry) + } + + pub fn round_slice_into(&self, dest: &mut Vec, digits: &[u8]) { + let (&d0, rest) = digits.split_first().unwrap_or((&0, &[])); + let digits = rest.iter().copied(); + let mut carry = 0; + let r0 = self.round_digit_with_carry(d0, &mut carry); + dest.push(r0); + extend_adding_with_carry(dest, digits, &mut carry); + if !carry.is_zero() { + dest.push(carry); + } + } + + #[allow(dead_code)] + pub fn round_slice_into_with_carry(&self, dest: &mut Vec, digits: &[u8], carry: &mut u8) { + let (&d0, rest) = digits.split_first().unwrap_or((&0, &[])); + let digits = rest.iter().copied(); + let r0 = self.round_digit_with_carry(d0, carry); + dest.push(r0); + + extend_adding_with_carry(dest, digits, carry); + } } diff --git a/src/with_std.rs b/src/with_std.rs index ace3082..0d0ef58 100644 --- a/src/with_std.rs +++ b/src/with_std.rs @@ -13,6 +13,7 @@ mod stdlib { num, ops, iter, + slice, str, string, i8, diff --git a/src/without_std.rs b/src/without_std.rs index deac488..3060568 100644 --- a/src/without_std.rs +++ b/src/without_std.rs @@ -24,6 +24,7 @@ mod stdlib { num, ops, iter, + slice, str, i8, f32,