bigdecimal/arithmetic/
cbrt.rs

Help
1//! Implementation of cube-root algorithm
2
3use crate::*;
4use num_bigint::BigUint;
5use rounding::NonDigitRoundingData;
6use stdlib::num::NonZeroU64;
7
8
9pub(crate) fn impl_cbrt_int_scale(n: &BigInt, scale: i64, ctx: &Context) -> BigDecimal {
10    let rounding_data = NonDigitRoundingData {
11        sign: n.sign(),
12        mode: ctx.rounding_mode(),
13    };
14
15    impl_cbrt_uint_scale((n.magnitude(), scale).into(), ctx.precision(), rounding_data)
16}
17
18/// implementation of cuberoot - always positive
19pub(crate) fn impl_cbrt_uint_scale(
20    n: WithScale<&BigUint>,
21    precision: NonZeroU64,
22    // contains sign and rounding mode
23    rounding_data: NonDigitRoundingData,
24) -> BigDecimal {
25    if n.is_zero() {
26        let biguint = BigInt::from_biguint(Sign::Plus, n.value.clone());
27        return BigDecimal::new(biguint, n.scale / 3);
28    }
29
30    // count number of digits in the decimal
31    let integer_digit_count = count_decimal_digits_uint(n.value);
32
33    // extra digits to use for rounding
34    let extra_rounding_digit_count = 4;
35
36    // required number of digits for precision and rounding
37    let required_precision = precision.get() + extra_rounding_digit_count;
38    let required_precision = 3 * required_precision;
39
40    // number of extra zeros to add to end of integer_digits
41    let mut exp_shift = required_precision.saturating_sub(integer_digit_count);
42
43    // effective scale after multiplying by 10^exp_shift
44    // (we've added that many trailing zeros after)
45    let shifted_scale = n.scale + exp_shift as i64;
46
47    let (mut new_scale, remainder) = shifted_scale.div_rem(&3);
48
49    match remainder.cmp(&0) {
50        Ordering::Greater => {
51            new_scale += 1;
52            exp_shift += (3 - remainder) as u64;
53        }
54        Ordering::Less => {
55            exp_shift += remainder.neg() as u64;
56        }
57        Ordering::Equal => {
58        }
59    }
60
61    // clone-on-write copy of digits
62    let mut integer_digits = stdlib::borrow::Cow::Borrowed(n.value);
63
64    // add required trailing zeros to integer_digits
65    if exp_shift > 0 {
66        arithmetic::multiply_by_ten_to_the_uint(
67            integer_digits.to_mut(), exp_shift
68        );
69    }
70
71    let result_digits = integer_digits.nth_root(3);
72    let result_digits_count = count_decimal_digits_uint(&result_digits);
73    debug_assert!(result_digits_count > precision.get());
74
75    let digits_to_trim = result_digits_count - precision.get();
76    debug_assert_ne!(digits_to_trim, 0);
77    debug_assert!((result_digits_count as i64 - count_decimal_digits_uint(&integer_digits) as i64 / 3).abs() < 2);
78
79    new_scale -= digits_to_trim as i64;
80
81    let divisor = ten_to_the_uint(digits_to_trim);
82    let (mut result_digits, remainder) = result_digits.div_rem(&divisor);
83
84    let remainder_digits = remainder.to_radix_le(10);
85    let insig_digit0;
86    let trailing_digits;
87    if remainder_digits.len() < digits_to_trim as usize {
88        // leading zeros
89        insig_digit0 = 0;
90        trailing_digits = remainder_digits.as_slice();
91    } else {
92        let (&d, rest) = remainder_digits.split_last().unwrap();
93        insig_digit0 = d;
94        trailing_digits = rest;
95    }
96
97    let insig_data = rounding::InsigData::from_digit_and_lazy_trailing_zeros(
98        rounding_data, insig_digit0, || { trailing_digits.iter().all(Zero::is_zero) }
99    );
100
101    // lowest digit to round
102    let sig_digit = (&result_digits % 10u8).to_u8().unwrap();
103    let rounded_digit = insig_data.round_digit(sig_digit);
104
105    let rounding_term = rounded_digit - sig_digit;
106    result_digits += rounding_term;
107
108    let result = BigInt::from_biguint(rounding_data.sign, result_digits);
109
110    BigDecimal::new(result, new_scale)
111}
112
113
114#[cfg(test)]
115mod test {
116    use super::*;
117    use stdlib::num::NonZeroU64;
118
119    macro_rules! impl_test {
120        ($name:ident; $input:literal => $expected:literal) => {
121            #[test]
122            fn $name() {
123                let n: BigDecimal = $input.parse().unwrap();
124                let value = n.cbrt();
125
126                let expected = $expected.parse().unwrap();
127                assert_eq!(value, expected);
128            }
129        };
130        ($name:ident; prec=$prec:literal; round=$round:ident; $input:literal => $expected:literal) => {
131            #[test]
132            fn $name() {
133                let ctx = Context::new(NonZeroU64::new($prec).unwrap(), RoundingMode::$round);
134                let n: BigDecimal = $input.parse().unwrap();
135                let value = n.cbrt_with_context(&ctx);
136
137                let expected = $expected.parse().unwrap();
138                assert_eq!(value, expected);
139            }
140        };
141    }
142
143    mod default {
144        use super::*;
145
146        impl_test!(case_0; "0.00" => "0");
147        impl_test!(case_1; "1.00" => "1");
148        impl_test!(case_1d001; "1.001" => "1.000333222283909495175449559955220102010284758197360454054345461242739715702641939155238095670636841");
149        impl_test!(case_10; "10" => "2.154434690031883721759293566519350495259344942192108582489235506346411106648340800185441503543243276");
150        impl_test!(case_13409d179789484375; "13409.179789484375" => "23.7575");
151        impl_test!(case_n59283293e25; "-59283293e25" => "-84006090355.84281237113712383191213626687332139035750444925827809487776780721673264524620270275301685");
152        impl_test!(case_94213372931en127; "94213372931e-127" => "2.112049945275324414051072540210070583697242797173805198575907094646677475250362108901530353886613160E-39");
153    }
154
155    impl_test!(case_prec15_down_10; prec=15; round=Down; "10" => "2.15443469003188");
156    impl_test!(case_prec6_up_0d979970546636727; prec=6; round=Up; "0.979970546636727" => "0.993279");
157
158    impl_test!(case_1037d495615705321421375_full; "1037.495615705321421375" => "10.123455");
159    impl_test!(case_1037d495615705321421375_prec7_halfdown; prec=7; round=HalfDown; "1037.495615705321421375" => "10.12345");
160    impl_test!(case_1037d495615705321421375_prec7_halfeven; prec=7; round=HalfEven; "1037.495615705321421375" => "10.12346");
161    impl_test!(case_1037d495615705321421375_prec7_halfup; prec=7; round=HalfUp; "1037.495615705321421375" => "10.12346");
162
163    impl_test!(case_0d014313506928855520728400001_full; "0.014313506928855520728400001" => "0.242800001");
164    impl_test!(case_0d014313506928855520728400001_prec6_down; prec=6; round=Down; "0.014313506928855520728400001" => "0.242800");
165    impl_test!(case_0d014313506928855520728400001_prec6_up; prec=6; round=Up; "0.014313506928855520728400001" => "0.242801");
166
167    impl_test!(case_4151902e20_prec16_halfup; prec=16; round=HalfUp; "4151902e20" => "746017527.6855992");
168    impl_test!(case_4151902e20_prec16_up; prec=16; round=Up; "4151902e20" => "746017527.6855993");
169    impl_test!(case_4151902e20_prec17_up; prec=17; round=Up; "4151902e20" => "746017527.68559921");
170    impl_test!(case_4151902e20_prec18_up; prec=18; round=Up; "4151902e20" => "746017527.685599209");
171    // impl_test!(case_4151902e20_prec18_up; prec=18; round=Up; "4151902e20" => "746017527.685599209");
172
173    impl_test!(case_1850846e201_prec14_up; prec=16; round=Up; "1850846e201" => "1.227788123885769e69");
174
175    impl_test!(case_6d3797558642427987505823530913e85_prec16_up; prec=160; round=Up; "6.3797558642427987505823530913E+85" => "3995778017e19");
176
177    impl_test!(case_88573536600476899341824_prec20_up; prec=20; round=Up; "88573536600476899341824" => "44576024");
178    impl_test!(case_88573536600476899341824_prec7_up; prec=7; round=Up; "88573536600476899341824" => "4457603e1");
179
180    impl_test!(case_833636d150970875_prec5_up; prec=5; round=Up; "833636.150970875" => "94.115000");
181    impl_test!(case_833636d150970875_prec5_halfup; prec=5; round=HalfUp; "833636.150970875" => "94.115");
182    impl_test!(case_833636d150970875_prec4_halfup; prec=4; round=HalfUp; "833636.150970875" => "94.12");
183    impl_test!(case_833636d150970875_prec20_up; prec=20; round=Up; "833636.150970875" => "94.115000");
184
185    #[cfg(property_tests)]
186    mod prop {
187        use super::*;
188        use proptest::*;
189        use num_traits::FromPrimitive;
190
191        proptest! {
192            #[test]
193            fn cbrt_of_cube_is_self(f: f64, prec in 15..50u64) {
194                // ignore non-normal numbers
195                prop_assume!(f.is_normal());
196
197                let n = BigDecimal::from_f64(f).unwrap().with_prec(prec);
198                let n_cubed = n.cube();
199                let x = n_cubed.cbrt();
200                prop_assert_eq!(x, n);
201            }
202        }
203    }
204}
bigdecimal/arithmetic/cbrt.rs

bigdecimal/arithmetic/
cbrt.rs
