summaryrefslogtreecommitdiffstats
diff options
context:
space:
mode:
authorTavian Barnes <tavianator@gmail.com>2010-11-04 23:26:58 -0400
committerTavian Barnes <tavianator@gmail.com>2010-11-04 23:28:51 -0400
commit120c523adc70af185b743e955837e1fe9e9a6785 (patch)
treeabb32a3f9cce65479b03ec6181308d728ade6a8a
parentcf10131df6b18ebb31a341f00be47273eaf51da5 (diff)
downloaddimension-120c523adc70af185b743e955837e1fe9e9a6785.tar.xz
Be more lenient about the root bracketing in dmnsn_bisect_root().
-rw-r--r--libdimension/polynomial.c15
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) {