From 120c523adc70af185b743e955837e1fe9e9a6785 Mon Sep 17 00:00:00 2001
From: Tavian Barnes <tavianator@gmail.com>
Date: Thu, 4 Nov 2010 23:26:58 -0400
Subject: Be more lenient about the root bracketing in dmnsn_bisect_root().

---
 libdimension/polynomial.c | 15 ++++++++++++++-
 1 file changed, 14 insertions(+), 1 deletion(-)

(limited to 'libdimension/polynomial.c')

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) {
-- 
cgit v1.2.3