vorot/roots

Library of well known algorithms for numerical root finding.

roots's Issues

More Informative Error

Firstly, amazing work on the crate!

I was wondering whether you would consider having more informative errors for when the search fails. For example, in the case of non-convergence, returning what the last estimate or last bracket was.

I realize that the different algorithms work differently and thus may store different underlying information, but I thought this could be particularly useful if one wanted to chain two root finding methods in case the first one fails.

Immutable Convergency

I've noticed that the root finding algorithms all take a mutable reference to the Convergency trait, which seems odd. In particular, I was hoping to have a const Convergency set for a module, but this is impossible to achieve because it must be mutable.

Are there any use case for a mutable convergency test other than in DebugConvergency? I would have thought that if DebugConvergency is the only reason why this is used, would it not be better to make it internally mutable through a RwLock (or similar)? As that is for debugging, the extra overhead should not be too much of a concern.

In fact, would it be better if the functions to use an associated type instead of &mut dyn Convergency? This would allow for mutable cases as well as immutable cases. I've created a small gist that shows what I mean.

`find_roots_quartic` does not find roots in a specific case.

I created a test case for find_roots_quartic where it should return two roots according to wolfram, but returns none.

Is this something I should expect from numerical errors? If yes, what would be a better method to solve such an equation?

Easier integration with rust_decimal OR fallible functions

Hello, thanks for the library. It works amazing. Perhaps my only source of pain is the combination of the following two facts:

1. there is no integration with rust_decimal::Decimal; so i am forced to use conversion functions e.g. Decimal::to_f64()
2. however, conversion functions are fallible and return Option<f64>;

Besides the clutter, my concern is that i'm usually forced to unwrap() (risking panic) or unwrap_or(meaningless) when instead I'd rather have a version of the root finding functions which accept a fallible evaluation function where None signals "user error, stop iterating".

Do you think it'd be possible to support such a use case?

Getting the number of roots

Hi,

Thanks for writing this library! I'm currently trying to get the number of roots to a quartic equation with this:

let v = roots::find_roots_quartic(a4, a3, a2, a1, a0);
let mut num_roots = 0;
match v {
    roots::Roots::Four(roots) => {
        num_roots = 4;
    },
    roots::Roots::Three(roots) => {
        num_roots = 3;
    },
    roots::Roots::Two(roots) => {
        num_roots = 2;
    },
    roots::Roots::One(roots) => {
        num_roots = 1;
    },
    _ => {num_roots = 0;}
}

Is there a better way to do this? Thanks!

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

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.

find_roots_eigen Not Working for Polynomial of Degree 4 or More

I was transforming a python code to rust, the python library numpy seems to only use eigen for their roots function. I wanted to use the same. In my test I want the roots x^5 -2.5x^4 +5.0x^3 -5.0x^2 +2.5x -0.5, so :

roots::find_roots_eigen(vec![-2.5,5.0,-5.0,2.5,-0.5]);

The function gives me 0.890896422653356 while I was expecting 0.5.

Also, in src/numerical/eigen.rs there is a test with fake value to pass.

let roots = dbg!(find_roots_eigen(vec![
            -3.75f64 / -14.0625f64,
            29.75f64 / -14.0625f64,
            4.0f64 / -14.0625f64,
            -16.0f64 / -14.0625f64,
        ]));

// (According to Wolfram Alpha, roots must be -1.1016116464173349f64, 0.9682783130840016f64)
// This means that this function cannot handle such small discriminant.
// But at least it finds the right number of them.
assert_eq!(roots[0], 0.9990584398692597f64);
assert_eq!(roots[1], 0.12511486301943303f64);

If you know it does not work, why not throw an error or something? Maybe there is a need to rework the algorithm?

Precision of square eq solver can be improved for a->0

HI! Instead of using (-b+-sqrt(d))/2/a you can use this formulae to get rid of 0/0 division in case a->0:

    sqD = sqrt(d);
    x1 = b > 0 ? (2 * c / (-b - sqD)) : (2 * c / (-b + sqD));
    x2 = -b / a - x1;

Thanx.

Very small values of `a` produce incorrect results in `find_roots_cubic`

Hi, it looks like when presented with very small values for a - that is, for curves that are almost but note quite quadratic - find_roots_cubic produces wildly incorrect results. The following test case illustrates the issue:

    use roots::*;

    let a = -0.000000000000000040410628481035;
    let b = 0.0126298310280606;
    let c = -0.100896606408756;
    let d = 0.0689539597036461;

    let x = 0.754710877053; // Expected root

    assert!((a*x*x*x + b*x*x + c*x + d).abs() < 0.001);

    let roots = find_roots_cubic(a, b, c, d);
    let roots = match roots {
        Roots::No(_)    => vec![],
        Roots::One(r)   => r.to_vec(),
        Roots::Two(r)   => r.to_vec(),
        Roots::Three(r) => r.to_vec(),
        Roots::Four(r)  => r.to_vec()
    };

    println!("{:?}", roots);
    assert!(roots.into_iter().any(|r| (r-x).abs() < 0.01));

For me, this outputs [312537357195212.44] as the result and fails the assertion. There appears to be a threshold around a=-0.0000000002 where the results become correct again.

I've been working around this by treating small values of a as indicating that the curve is quadratic, which works for my use case (solving for points on bezier curves, which are quite commonly very close to being quadratic) but I'm not sure is good enough for a general solution where the other coefficients might also have very small values.

`find_roots_quadratic` susceptible to overflow/underflow

A newly published paper explicitly calls out the find_roots_quadratic() function of this library as being susceptible to underflow/overflow, saying:

4.5 Rust “roots” module
The roots module of the Rust language offers the function roots::find_roots_quadratic() to solve quadratic equations.
According to its source code, the algorithm used handles separately the linear case. The discriminant is computed with the
naive formula without any extra precision, and no provision is made to avoid underflows and overflows.

In addition, it lists 18 cases where that function gives the wrong number of solutions:

Might be worth taking a look!
Link to paper: https://cnrs.hal.science/hal-04116310
Link to github: https://github.com/goualard-f/QuadraticEquation.jl

Disable debug messages

Call roots::find_roots_quartic outputs a lot of dbg messages on console.

Potential bug in find_roots_quartic

The logic for determining whether there are no roots for the quartic equation seems wrong. Specifically, on the line below:

roots/src/analytical/quartic.rs

Line 166 in ad6ec95

let no_roots = pp > F::zero() || dd > F::zero();

This condition check for pp is only valid if the determinant is greater than zero. So the correct check should be this:

 let no_roots = discriminant > F::zero() && (pp > F::zero() || dd > F::zero()); 

A test case that can demonstrate this bug is the following:

a0 = {f64} 1.5971379368787519
a1 = {f64} -5.7774307699949397
a2 = {f64} 7.9323705711747614
a3 = {f64} -4.8721513473605924
a4 = {f64} 1.1248467624839498

where the correct real roots are 1.22591 and 1.25727 (wolfram alpha link) but the correct implementation will give no roots.

Let me know what you think and I can submit a PR. Thank you for your great work.

roots::numerical::SearchError is private

The numerical algorithms return SearchError as part of the Result; however the crate does not expose SearchError. As a result, I can call the numerical algorithm functions, but if I wish to pass the result to another function, I cannot actually implement a function with the return type Result<f64, SearchError>.

Let me know if you would like a PR.

Using generics instead of trait objects?

Hello!

First of all, nice job implementing all of these algorithms in Rust :)

Now on to my issue: I started using your library and noticed that it uses trait objects quite heavily for passing functions and convergencies. Now, given that a function has to be evaluated rather often during computation, it might be useful for performance reasons to switch to generics for this. Have you considered that? I wonder what the implications would be.

Also it might be more ergonomic for users, as it means you do not have to pass in references to closures and trait objects, but can just pass them by value. For functions, you still can use a reference, as they implement Fn transitively. So I'd argue it is strictly more general to do it like this. You may even implement Convergency on references to keep this change backwards-compatible.

I would be happy to provide a pull request along those lines.

FeatureSuggestion: FnMut for closures

Hi,

I have found myself a few times using roots (nice crate! thx!) in situations where:

• Each step of the iteration involves running agains a large structure
• I can re-evaluate my target function by doing minor modificaitons to my large structure

However, since e..g roots::find_root_brent requires a Fn I cannot pass a closure which captures a &mut ref to my struct and modify it in each iteration - instead I have to clone it every time to get a Fn closure instead of FnMut.

Would it be a good idea to accept FnMut closures for all the root functions or am I missing something?

Publish new version to crates.io

Version 0.0.7 still has the issue fixed by #19. Can you publish a new version to crates.io?

Timeframe for publishing 0.0.5

Thanks for merging my PR. Is there a timeframe for the next crates.io release?

`find_roots_quartic` returns incorrect roots

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!