diff --git a/src/arithmetic/cbrt.rs b/src/arithmetic/cbrt.rs index 45e63b0..3ed58cc 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 as u64); - // 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() + 1); + + 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) } @@ -128,6 +150,32 @@ mod test { 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 { diff --git a/src/arithmetic/mod.rs b/src/arithmetic/mod.rs index c321f70..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 { diff --git a/src/lib.rs b/src/lib.rs index 310d12e..30129af 100644 --- a/src/lib.rs +++ b/src/lib.rs @@ -761,11 +761,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 @@ -1221,6 +1217,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) @@ -1310,6 +1319,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)]