rpoly/
lib.rs

1#![doc = include_str!("../README.md")]
2
3mod rpoly_ak1;
4mod rpoly_types;
5
6use rpoly_ak1::{rpoly_ak1, LEADING_COEFFICIENT_ZERO_NUM, NOT_CONVERGENT_NUM};
7pub use rpoly_types::{
8    RpolyComplex, RpolyError,
9    RpolyError::{RpolyLeadingCoefficientZero, RpolyNotConvergent},
10    RpolyRoots,
11};
12
13/// Solve the polynomial given by:
14///
15/// a\[0\]*x^(n-1) + a\[1\]*x^(n-2) + ... + a\[n-2\]*x + a\[n-1\] = 0
16///
17/// where `n == a.len() == MDP1`.
18///
19/// # Errors
20///
21/// - [`RpolyLeadingCoefficientZero`]: If `a[0] == 0.0`.
22/// - [`RpolyNotConvergent`]: If the method fails to converge.
23///
24/// # Examples
25///
26/// ```
27/// # use rpoly::rpoly;
28/// #
29/// let roots = rpoly(&[1.0, -3.0, 2.0]).unwrap();  // x^2 - 3x + 2 = 0
30/// assert_eq!(roots.root_count(), 2);
31/// ```
32///
33pub fn rpoly<const MDP1: usize>(a: &[f64; MDP1]) -> Result<RpolyRoots<MDP1>, RpolyError> {
34    let mut degree = MDP1 - 1;
35    let mut zeror = [0.0; MDP1];
36    let mut zeroi = [0.0; MDP1];
37    rpoly_ak1(a, &mut degree, &mut zeror, &mut zeroi);
38    match degree {
39        LEADING_COEFFICIENT_ZERO_NUM => Err(RpolyLeadingCoefficientZero),
40        NOT_CONVERGENT_NUM => Err(RpolyNotConvergent),
41        _ => {
42            let mut roots = [RpolyComplex::default(); MDP1];
43            for i in 0..degree {
44                roots[i].re = zeror[i];
45                roots[i].im = zeroi[i];
46            }
47            Ok(RpolyRoots(roots))
48        }
49    }
50}
51
52#[test]
53fn test_rpoly() {
54    let roots = rpoly(&[
55        1.0,
56        -687712.0,
57        -3.12295507665e11,
58        7.503861976776e12,
59        3.803023224e15,
60        0.0,
61    ])
62    .unwrap();
63
64    assert!(roots.root_count() == 5);
65    for root in roots {
66        assert!(root.im == 0.0);
67        // dbg!(root);
68    }
69}
rpoly/lib.rs

rpoly/
lib.rs
