bigdecimal/arithmetic/
exp.rs

1//! All the routines for calculating exp(x)
2//!
3
4use crate::*;
5use super::*;
6
7
8/// Calculate e^n
9pub(crate) fn impl_exp(n: BigDecimalRef, ctx: &Context) -> BigDecimal {
10    use arithmetic::division::scaled_uint_division_into;
11    use arithmetic::inverse::impl_inverse_uint_scale;
12
13    if n.is_zero() {
14        return BigDecimal::one();
15    }
16
17    // n ~= x * 2^k
18    let (x, k) = factor_two_to_k_scale(WithScale { value: n.digits, scale: n.scale });
19
20    let target_precision = ctx.precision().get();
21
22    // always at least 3 u64's worth of bits
23    let target_precision_bits = digit_to_bit_count(target_precision).max(64 * 3) + k as u64;
24
25    let mut num = x.clone();
26    let mut den = WithScale { value: 1u8.into(), scale: 0 };
27
28    // sum = 1 + x
29    let mut sum = WithScale { value: 1u8.into(), scale: 0 };
30    addition::addassign_scaled_biguint(&mut sum, num.as_ref());
31
32    // we should have form `1.xxxxx` so only one digit is integer part
33    if true {
    {
        match (&sum.count_int_digits(), &1) {
            (left_val, right_val) => {
                if !(*left_val == *right_val) {
                    let kind = ::core::panicking::AssertKind::Eq;
                    ::core::panicking::assert_failed(kind, &*left_val,
                        &*right_val, ::core::option::Option::None);
                }
            }
        }
    };
};debug_assert_eq!(sum.count_int_digits(), 1);
34
35    let mut delta: WithScale<BigUint> = Default::default();
36
37    // assuming linear convergence, should break after N; we loop
38    // through 2*N for safety
39    let stop = (target_precision * 2).max(10);
40
41    for i in 2..stop {
42        // each loop iteration:
43        //    num = x^i
44        //    den = factorial(i)
45        //    delta = num / den
46        //    sum += delta
47        num.mulassign_scaled_biguint(&x);
48        den.value *= i;
49        remove_trailing_zeros(&mut den, &[1]);
50        scaled_uint_division_into(&mut delta, &num, &den, target_precision_bits - i);
51        sum.addassign_scaled_biguint(&delta);
52
53        // we have converged if number of leading zeros in delta is
54        // larger than the target precision
55        let leading_zero_count = delta.count_int_digits().neg();
56        if leading_zero_count > target_precision as i64 {
57            break;
58        }
59    }
60
61    // reuse 'delta' as scratchpad
62    let mut tmp = delta.value;
63
64    // at this point: sum = exp(n / 2^k)
65    //
66    // we now square it 'k' times to get final result
67    arithmetic::pow::pow_2_k_scaled_biguint(
68        &mut sum, &mut tmp, k, target_precision_bits
69    );
70
71    let result = BigDecimal::from(sum);
72    if n.sign == Sign::Minus {
73        return ctx.invert(&result);
74    } else {
75        return ctx.round_decimal(result);
76    }
77}
78
79/// Factor scaled biguint by 2^k
80///
81/// Returns pair of sacled-BigUint in range [0.0, 0.5] and 'k', the
82/// number of times to square the BigUint to return the number to
83/// the original value.
84///
85fn factor_two_to_k_scale(n: WithScale<&BigUint>) -> (WithScale<BigUint>, u16) {
86    let log2_n = n.value.bits() as f64;
87    let log2_s = (n.scale as f64) * LOG2_10;
88    let k = 1.0 + log2_n - log2_s;
89    if k <= 0.5 {
90        let r = n.value.clone();
91        return ((r, n.scale).into(), 0);
92    }
93
94    let k = k.ceil() as u64;
95    let mut x = BigUint::from(5u8).pow(k);
96    x *= n.value;
97
98    let mut result = WithScale {
99        value: x,
100        scale: n.scale + k as i64,
101    };
102
103    // strip trailing zeros at different speeds (last one must be '1')
104    remove_trailing_zeros(&mut result, &[8, 1]);
105
106    (result, k.as_())
107}
108
109/// Remove trailing zeros by divmoding n by given powers of ten
110///
111/// Multiple powers may be chosen to speed up removal (should end with '1'
112/// to remove all zero, i.e. mod-10)
113///
114fn remove_trailing_zeros<'a>(
115    n: &mut WithScale<BigUint>,
116    powers: impl IntoIterator<Item=&'a u8>,
117) {
118    for &i in powers.into_iter() {
119        if true {
    if !(i < 20) { ::core::panicking::panic("assertion failed: i < 20") };
};debug_assert!(i < 20);
120
121        let s = 10u64.pow(i as u32);
122        while (&n.value % s).is_zero() {
123            n.value /= s;
124            n.scale -= i as i64;
125        }
126    }
127}
128
129#[cfg(test)]
130mod test {
131    use super::*;
132
133    include!("exp.tests.rs");
134}
bigdecimal/arithmetic/exp.rs

bigdecimal/arithmetic/
exp.rs
