fastnum/decimal/dec/math/
nth_root.rs

1use crate::decimal::{
2    dec::{
3        convert::to_f64,
4        intrinsics::Intrinsics,
5        math::{add::add, div::div, mul::mul},
6        parse::{from_f64, from_u32},
7    },
8    utils::types,
9    Decimal,
10};
11
12type D<const N: usize> = Decimal<N>;
13
14#[inline]
15pub(crate) const fn nth_root<const N: usize>(d: D<N>, n: u32) -> D<N> {
16    if d.is_nan() {
17        return d.op_invalid();
18    }
19
20    if d.is_zero() || d.is_one() {
21        return d;
22    }
23
24    if d.is_negative() {
25        return d.signaling_nan();
26    }
27
28    if d.is_infinite() {
29        return d;
30    }
31
32    nth_root_newton(d, n)
33}
34
35#[inline]
36const fn nth_root_newton<const N: usize>(d: D<N>, n: u32) -> D<N> {
37    let approx_f64 = to_f64(d);
38    let guess = types::f64::sqrt(approx_f64);
39
40    let mut result = from_f64(guess).compound(&d);
41
42    let mut result_next;
43
44    let n_minus_one = from_u32(n - 1);
45    let one_div_n = D::ONE.div(from_u32(n));
46    let mut x_n;
47    let mut j;
48    let mut i = 1;
49
50    while result.is_ok() && i < Intrinsics::<N>::SERIES_MAX_ITERATIONS {
51        x_n = result;
52        j = n - 2;
53
54        while j > 0 {
55            x_n = mul(x_n, result);
56            j -= 1;
57        }
58
59        result_next = mul(one_div_n, add(mul(n_minus_one, result), div(d, x_n)));
60
61        if result.eq(&result_next) {
62            break;
63        }
64
65        result = result_next;
66        i += 1;
67    }
68
69    result
70}