From 97a9efe5091a5c483f8081d60dccf6029be7ecc4 Mon Sep 17 00:00:00 2001 From: Tavian Barnes Date: Tue, 9 Nov 2010 16:36:08 -0500 Subject: 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. --- libdimension/polynomial.c | 13 +++++++++---- 1 file changed, 9 insertions(+), 4 deletions(-) (limited to 'libdimension/polynomial.c') 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; } -- cgit v1.2.3