rgsl/types/
polynomial.rs

1//
2// A rust binding for the GSL library by Guillaume Gomez (guillaume1.gomez@gmail.com)
3//
4
5/*!
6# General Polynomial Equations
7
8The roots of polynomial equations cannot be found analytically beyond the special cases of the quadratic, cubic and quartic equation. The algorithm
9described in this section uses an iterative method to find the approximate locations of roots of higher order polynomials.
10!*/
11
12use crate::Value;
13use ffi::FFI;
14
15ffi_wrapper!(
16    PolyComplexWorkspace,
17    *mut sys::gsl_poly_complex_workspace,
18    gsl_poly_complex_workspace_free
19);
20
21impl PolyComplexWorkspace {
22    /// This function allocates space for a gsl_poly_complex_workspace struct and a workspace suitable for solving a polynomial with n coefficients
23    /// using the routine gsl_poly_complex_solve.
24    ///
25    /// The function returns a pointer to the newly allocated gsl_poly_complex_workspace if no errors were detected, and a null pointer in the case
26    /// of error.
27    #[doc(alias = "gsl_poly_complex_workspace_alloc")]
28    pub fn new(n: usize) -> Option<Self> {
29        let tmp = unsafe { sys::gsl_poly_complex_workspace_alloc(n) };
30
31        if tmp.is_null() {
32            None
33        } else {
34            Some(Self::wrap(tmp))
35        }
36    }
37
38    /// This function computes the roots of the general polynomial P(x) = a_0 + a_1 x + a_2 x^2 + ... + a_{n-1} x^{n-1} using balanced-QR reduction
39    /// of the companion matrix. The parameter n specifies the length of the coefficient array. The coefficient of the highest order term must be
40    /// non-zero. The function requires a workspace w of the appropriate size. The n-1 roots are returned in the packed complex array z of length
41    /// 2(n-1), alternating real and imaginary parts.
42    ///
43    /// The function returns Success if all the roots are found. If the QR reduction does not converge, the error handler is invoked with an error
44    /// code of Failed. Note that due to finite precision, roots of higher multiplicity are returned as a cluster of simple roots with reduced
45    /// accuracy. The solution of polynomials with higher-order roots requires specialized algorithms that take the multiplicity structure into
46    /// account (see e.g. Z. Zeng, Algorithm 835, ACM Transactions on Mathematical Software, Volume 30, Issue 2 (2004), pp 218–236).
47    #[doc(alias = "gsl_poly_complex_solve")]
48    pub fn solve(&mut self, a: &[f64], z: &mut [f64]) -> Result<(), Value> {
49        let ret = unsafe {
50            sys::gsl_poly_complex_solve(
51                a.as_ptr(),
52                a.len() as _,
53                self.unwrap_unique(),
54                z.as_mut_ptr(),
55            )
56        };
57        result_handler!(ret, ())
58    }
59}