`find_roots_quartic` misses double-roots, is unstable when multiplying all coefficients by a constant

[ryan-williams]

Code

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

use roots::{find_roots_quartic, Roots};

/// Synthesize coefficients from 4 real roots (and a `scale` parameter, equal to `a_4`), call `find_roots_quartic`
fn get_roots(r0: f64, r1: f64, r2: f64, r3: f64, scale: f64) -> Roots<f64> {
    let unscaled_coeffs = [
        1.,
        -(r0 + r1 + r2 + r3),
        r0*r1 + r0*r2 + r0*r3 + r1*r2 + r1*r3 + r2*r3,
        -(r0*r1*r2 + r0*r1*r3 + r0*r2*r3 + r1*r2*r3),
        r0*r1*r2*r3,
    ];
    let [ a4, a3, a2, a1, a0 ] = unscaled_coeffs.map(|c| c * scale);
    let f = |a: f64, n: usize| {
        let x_str = match n {
            0 => "",
            1 => "x",
            2 => "x²",
            3 => "x³",
            4 => "x⁴",
            _ => panic!("Unsupported exponent: {}", n),
        };
        let coeff_str = if a.abs() == 1. { "".to_string() } else { format!("{}", a.abs()) };
        let prefix = if n == 4 { 
            if a < 0. { "-" } else { "" }
        } else { 
            if a < 0. { " - " } else { " + " }
        };
        format!("{}{}{}", prefix, coeff_str, x_str)
    };
    println!("  Scaled by {}:", scale);
    println!("    {}{}{}{}{}", f(a4, 4), f(a3, 3), f(a2, 2), f(a1, 1), f(a0, 0));
    println!("    coeffs: {:?}", [ a4, a3, a2, a1, a0 ]);
    find_roots_quartic(a4, a3, a2, a1, a0)
}

fn scaled_roots<const N: usize>(r0: f64, r1: f64, r2: f64, r3: f64, scales: [f64; N]) {
    println!("Expected roots: {} {} {} {}:", r0, r1, r2, r3);
    for scale in scales {
        let roots = get_roots(r0, r1, r2, r3, scale);
        println!("    Computed roots: {:?}", roots);
    }
}

fn main() {
    // Compute the same roots, but with coefficients scaled by [ 0.1, 1, 10 ]
    scaled_roots(0.1, 0.1, 1., 1., [ 0.1, 1., 10. ]);
}

1. Start with known roots, generate coefficients
2. Scale coefficients by 0.1, 1., or 10.
   • "scale" becomes a4 (the coefficient of x⁴)
   • shouldn't effect the computed root values
3. Run find_roots_quartic, print Roots result

Open in "Rust Explorer"

Output:

Expected roots: 0.1 0.1 1 1:
  Scaled by 0.1:
    0.1x⁴ - 0.22000000000000003x³ + 0.14100000000000001x² - 0.022000000000000006x + 0.0010000000000000002
    coeffs: [0.1, -0.22000000000000003, 0.14100000000000001, -0.022000000000000006, 0.0010000000000000002]
    Computed roots: Two([0.9999999882195979, 1.0000000117804024])  ❌ misses "0.1" double-root
  Scaled by 1:
    x⁴ - 2.2x³ + 1.4100000000000001x² - 0.22000000000000003x + 0.010000000000000002
    coeffs: [1.0, -2.2, 1.4100000000000001, -0.22000000000000003, 0.010000000000000002]
    Computed roots: Four([0.09999999283067545, 0.1000000071693245, 0.9999999928306755, 1.0000000071693247])  ✅ close enough
  Scaled by 10:
    10x⁴ - 22x³ + 14.100000000000001x² - 2.2x + 0.10000000000000002
    coeffs: [10.0, -22.0, 14.100000000000001, -2.2, 0.10000000000000002]
    Computed roots: No([])  ❌

• First iteration ("Scaled by 0.1") misses the double-root at 0.1 ❌
• Second iteration ("Scaled by 1") returns four roots, with ≈1e-7 relative accuracy ✅
• Third iteration ("Scaled by 10") doesn't return any roots ❌

Possibly related to #23.