fastnum/decimal/dec/math/
nth_root.rs1use 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}