find_roots_quartic returns incorrect roots

[ryan-williams]

Similar to #23, but in this case find_roots_quartic returns substantially incorrect values:

/*
[dependencies]
roots = "0.0.8"
*/

use roots::find_roots_quartic;

fn main() {
    let a4 = 0.000000030743755847066437;
    let a3 = 0.000000003666731306801131;
    let a2 = 1.0001928389119579;
    let a1 = 0.000011499702220469921;
    let a0 = -0.6976068572771268;
    let f = |x: f64| {
        let x2 = x * x;
        a4*x2*x2 + a3*x2*x + a2*x2 + a1*x + a0
    };
    let roots = find_roots_quartic(a4, a3, a2, a1, a0);
    for root in roots.as_ref() {
        println!("  x: {}, f(x): {}", root, f(*root));
    }
}

Output (open in Rust Explorer):

x: -5.106021787341158, f(x): 25.37884091853862
x: 5.106010194296296, f(x): 25.37884091853862

Note that plugging the roots back into f yields 25.37884091853862 (supposed to be 0).

Wolfram Alpha gives correct values:

• $x$ ≈ -0.835153846196954
• $x$ ≈ 0.835142346155438

Looking at the coefficients, we can see that $a_4$, $a_3$, and $a_1$ are negligibly small; it's basically a parabola defined by $a_2 ≈ 1$ and $a_0 ≈ 0.69$, so the roots are $≈\sqrt(0.69) ≈ 0.83$. The $±5.1$ returned by find_roots_quartic is way off.

Use case

I'm using find_roots_quartic to find intersections between ellipses (cf. runsascoded/shapes), which reduces to solving a quartic equation. The coefficients/quartic above relate to intersecting a unit circle (centered at the origin) with this ellipse:

{
    c: { x: -1.100285308561806, y: -1.1500279763995946e-5 },
    r: { x: 1.000263820108834, y: 1.0000709021402923 }
}

center at $≈(-1.1, 0)$, radii $≈(1, 1)$:

The quartic above is actually solving for $y$-coordinates, in terms of this picture; $±5.1$ is pretty nonsensical for these shapes!

[vorot]

This sort of errors may be caused by the underflow/overflow when calculating the discriminant, similar to #31.

A solution could be to increase the precision of variables used for that.

Anyway, an unit test for this case is welcome.