vorot / roots

Library of well known algorithms for numerical root finding.
BSD 2-Clause "Simplified" License
50 stars 18 forks source link

Potential bug in find_roots_quartic #16

Closed jingnanshi closed 2 years ago

jingnanshi commented 2 years ago

The logic for determining whether there are no roots for the quartic equation seems wrong. Specifically, on the line below: https://github.com/vorot/roots/blob/ad6ec955452429bc4e38a4fb9d25936b66b4b11f/src/analytical/quartic.rs#L166

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.

jingnanshi commented 2 years ago

PR merged. Closing this issue.