vorot / roots Goto Github PK

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?

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.

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?

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

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?

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.

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

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!

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.

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

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.

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.

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.

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.

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