summaryrefslogtreecommitdiffstats
path: root/libdimension
diff options
context:
space:
mode:
Diffstat (limited to 'libdimension')
-rw-r--r--libdimension/polynomial.c13
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;
}