Stability fix for dmnsn_bisect_root().

When one of the bounds is close to a different root, make sure the result
is more accurate than that bound.  Otherwise we find the wrong root, and
eventually hang.  This could be seen with a 1920x1080 render of demo.pov,
for example.

1 files changed, 9 insertions, 4 deletions

diff --git a/libdimension/polynomial.c b/libdimension/polynomial.c
index e242bed..b68d440 100644
--- a/libdimension/polynomial.c
+++ b/libdimension/polynomial.c
@@ -219,14 +219,15 @@ dmnsn_root_bound(double poly[], size_t degree)
 static inline double
 dmnsn_bisect_root(const double poly[], size_t degree, double min, double max)
 {
+  if (max - min < dmnsn_epsilon)
+    return min;
+
   double evmin = dmnsn_evaluate_polynomial(poly, degree, min);
   double evmax = dmnsn_evaluate_polynomial(poly, degree, max);
+  double initial_evmin = evmin, initial_evmax = evmax;
   double mid = 0.0, evmid;
   int lastsign = -1;
 
-  if (max - min < dmnsn_epsilon)
-    return min;
-
   do {
     mid = (min*evmax - max*evmin)/(evmax - evmin);
     evmid = dmnsn_evaluate_polynomial(poly, degree, mid);
@@ -261,7 +262,11 @@ dmnsn_bisect_root(const double poly[], size_t degree, double min, double max)
     }
 
     lastsign = sign;
-  } while (fabs(evmid) >= fabs(mid)*dmnsn_epsilon);
+  } while (fabs(evmid) >= fabs(mid)*dmnsn_epsilon
+           /* These conditions improve stability when one of the bounds is
+              close to a different root than we are trying to find */
+           || fabs(evmid) > fabs(initial_evmin)
+           || fabs(evmid) > fabs(initial_evmax));
 
   return mid;
 }