Problem calculating fourth-degree poly roots

[m2-farzan]

Unfortunately the library fails to evaluate roots of fourth-degree equations. Take this code for example:

use rustnomial::{Polynomial};

fn main() {
    let poly = Polynomial::<f64>::new(vec![1.0, 1.0, 1.0, 1.0, 1.0]);
    let roots = poly.roots();
    match roots {
        rustnomial::Roots::OnlyRealRoots(_) => {
            println!("roots() returned `Only Real Roots`");
        },
        _ => {
            println!("roots() returned some other type.");
        }
    }
}

Running this with rustnomial = "0.2.0" prints:

roots() returned `Only Real Roots`

while in fact x^4 + x^3 + x^2 + x + 1 has no real roots. The vector V in OnlyRealRoots(V) is empty too.

Some other times, the library returns OnlyRealRoots(V) with V containing only real roots, skipping complex ones.

Maybe I'm missing something here, but my understanding was that the library should have returned ManyComplexRoots in these cases.

BTW, thanks for the nice package. I'm using it in an educational project.

[philippeitis]

The issue in this case is that I'm using the roots library to find roots, and this library is limited to strictly real roots. I do implement the algebraic solutions for cubic polynomials and below - quartic polynomials are on my todo-list, but the equations are quite long and a pain to implement. As no closed-form solutions exist for the roots of polynomials with degree 5 or above, I do not return complex roots for these.

I've fixed OnlyRealRoots(V) having an empty V on the latest master commit. Complex root finding for quartic polynomials (and root finding if the polynomial can be reduced to a quartic polynomial with roots found by the roots library) can be implemented (https://en.wikipedia.org/wiki/Quartic_function). I'd be more than happy to accept a PR for this.

Implementing general root-finding is outside the scope of this library, and a good solution requires a lot of work and time to implement properly. If you find a library that manages to do these things, I'd be happy to use it.

[m2-farzan]

So I suggest adding a note in the README, clarifying that the root-finding functionality is limited to third/fourth degree polynomials. It'll save visitors some time.

About solving quartic functions, I just tried to implement it in Rust and it seems to have worked out. I also tried to port it to rustnomial but Rust is really Greek to me and I couldn't get a grasp on the existing code base, and I also had trouble trying to make my code work with generic template types instead of only f64. So I'll just copy-paste my code here. I understand that maybe it's totally worthless, but if you find it good, you could make the required changes and add it to the library. Again, maybe it's just easier to write from scratch...

#![allow(non_snake_case)]

use num::Complex;
use std::f64::consts::PI;

// Solve a x^2 + b x + c = 0
fn quadratic_roots(a: f64, b: f64, c: f64) -> [Complex<f64>; 2]
{
    let D = Complex::new(b.powf(2.0) - 4.0*a*c, 0.0);
    let x1 = (-b + D.sqrt()) / (2.0 * a);
    let x2 = (-b - D.sqrt()) / (2.0 * a);
    return [x1, x2];
}

// Solve a x^4 + b x^3 + c x^2 + d x + e = 0
fn quartic_roots(coefficients: [f64; 5]) -> [Complex<f64>; 4]
{
    let [a, b, c, d, e] = coefficients;

    let D0 = c.powf(2.0) - 3.0*b*d + 12.0*a*e;
    let D1 = 2.0*c.powf(3.0) - 9.0*b*c*d + 27.0*b.powf(2.0)*e + 27.0*a*d.powf(2.0) - 72.0*a*c*e;
    let _D = -1.0/27.0 * (D1.powf(2.0) - 4.0*D0.powf(3.0));

    if D0 == 0.0 && D1 == 0.0 // special case; three repeated roots
    {
        let x0 = quadratic_roots(6.0 * a, 3.0 * b, c)[0];
        let x1 = -b/a - 3.0*x0;
        return [x0, x0, x0, x1];

    }
    
    let p = (8.0*a*c - 3.0*b.powf(2.0)) / (8.0*a.powf(2.0));
    let q = (b.powf(3.0) - 4.0*a*b*c + 8.0*a.powf(2.0)*d) / (8.0*a.powf(3.0));
    
    let mut Q = ((
        D1 + Complex::new(D1.powf(2.0) - 4.0*D0.powf(3.0), 0.0).sqrt()
    ) / 2.0).powf(1.0 / 3.0);
    if Q == Complex::new(0.0, 0.0) // spacial case; use the other square root for D1^2-4*D0^3
    {
        Q = ((
            D1 - Complex::new(D1.powf(2.0) - 4.0*D0.powf(3.0), 0.0).sqrt()
        ) / 2.0).powf(1.0 / 3.0);
    }

    let mut S = 0.5 * (-2.0/3.0 * p + 1.0 /(3.0*a) * (Q + D0/Q)).sqrt();
    if S == Complex::new(0.0, 0.0) // spacial case; use another cube root for Q
    {
        Q = Q * Complex::from_polar(&1.0, &(2.0 * PI / 3.0));
        S = 0.5 * (-2.0/3.0 * p + 1.0 /(3.0*a) * (Q + D0/Q)).sqrt();
    }

    let x1 = -b/(4.0*a) - S + 0.5 * (-4.0*S.powf(2.0) - 2.0*p + q/S).sqrt();
    let x2 = -b/(4.0*a) - S - 0.5 * (-4.0*S.powf(2.0) - 2.0*p + q/S).sqrt();
    let x3 = -b/(4.0*a) + S + 0.5 * (-4.0*S.powf(2.0) - 2.0*p - q/S).sqrt();
    let x4 = -b/(4.0*a) + S - 0.5 * (-4.0*S.powf(2.0) - 2.0*p - q/S).sqrt();

    return [x1, x2, x3, x4];
}


fn main() {
    let X = quartic_roots([1.0, 1.0, 1.0, 1.0, 1.0]);
    for i in 0..4 {
        println!("{}", X[i]);
    }
}


#[test]
fn test_quartic_random() {
    let random_coefficients = [
        [0.548914, -0.917383, -0.693302, 0.634482, 0.153504],
        [-0.449297, 0.150007, 0.056093, 0.780327, -0.205118],
        [-0.628965, -0.280957, -0.736180, -0.432911, 0.339487],
        [-0.347198, 0.201530, 0.881887, 0.669201, -0.075353]
    ];

    for coefficients in random_coefficients {
        let roots = quartic_roots(coefficients);

        for x in roots {
            let [a, b, c, d, e] = coefficients;
            let y = a*x.powf(4.0) + b*x.powf(3.0) + c*x.powf(2.0) + d*x + e;
            assert!(y.norm() < 1e-9);
        }
    }
}

#[test]
fn test_quartic_special_case_factorable() {
    let roots = quartic_roots([1.0, 4.0, 6.0, 4.0, 1.0]);
    assert_eq!(roots, [Complex::new(-1.0, 0.0); 4]);
}