diff options
| author | Tavian Barnes <tavianator@gmail.com> | 2010-11-04 23:26:58 -0400 |
|---|---|---|
| committer | Tavian Barnes <tavianator@gmail.com> | 2010-11-04 23:28:51 -0400 |
| commit | 120c523adc70af185b743e955837e1fe9e9a6785 (patch) | |
| tree | abb32a3f9cce65479b03ec6181308d728ade6a8a | |
| parent | cf10131df6b18ebb31a341f00be47273eaf51da5 (diff) | |
| download | dimension-120c523adc70af185b743e955837e1fe9e9a6785.tar.xz | |
Be more lenient about the root bracketing in dmnsn_bisect_root().
| -rw-r--r-- | libdimension/polynomial.c | 15 |
1 files changed, 14 insertions, 1 deletions
diff --git a/libdimension/polynomial.c b/libdimension/polynomial.c index dc4cb6a..e242bed 100644 --- a/libdimension/polynomial.c +++ b/libdimension/polynomial.c @@ -231,7 +231,20 @@ dmnsn_bisect_root(const double poly[], size_t degree, double min, double max) mid = (min*evmax - max*evmin)/(evmax - evmin); evmid = dmnsn_evaluate_polynomial(poly, degree, mid); int sign = dmnsn_signbit(evmid); - if (sign == dmnsn_signbit(evmax)) { + + if (mid < min) { + /* This can happen due to numerical instability in the root bounding + algorithm, so behave like the normal secant method */ + max = min; + evmax = evmin; + min = mid; + evmin = evmid; + } else if (mid > max) { + min = max; + evmin = evmax; + max = mid; + evmax = evmid; + } else if (sign == dmnsn_signbit(evmax)) { max = mid; evmax = evmid; if (sign == lastsign) { |