From a0a11f34d75eef9971a870fd51dd6d08f7708311 Mon Sep 17 00:00:00 2001 From: Tavian Barnes Date: Fri, 21 Dec 2012 19:30:30 -0500 Subject: Fix the polynomial solver when the bounds are exact. --- libdimension/polynomial.c | 8 +++++++- 1 file changed, 7 insertions(+), 1 deletion(-) (limited to 'libdimension/polynomial.c') diff --git a/libdimension/polynomial.c b/libdimension/polynomial.c index a08318b..9491861 100644 --- a/libdimension/polynomial.c +++ b/libdimension/polynomial.c @@ -214,6 +214,12 @@ dmnsn_bisect_root(const double poly[], size_t degree, double min, double max) { double evmin = dmnsn_polynomial_evaluate(poly, degree, min); double evmax = dmnsn_polynomial_evaluate(poly, degree, max); + + /* Handle equal bounds, and equal values at the bounds. */ + if (dmnsn_unlikely(fabs(evmax - evmin) < dmnsn_epsilon)) { + return (min + max)/2.0; + } + double evinitial = dmnsn_min(fabs(evmin), fabs(evmax)); double mid, evmid; int lastsign = 0; @@ -425,7 +431,7 @@ dmnsn_polynomial_print(FILE *file, const double poly[], size_t degree) if (i < degree) { fprintf(file, (poly[i] >= 0.0) ? " + " : " - "); } - fprintf(file, "%.15g", fabs(poly[i])); + fprintf(file, "%.17g", fabs(poly[i])); if (i >= 2) { fprintf(file, "*x^%zu", i); } else if (i == 1) { -- cgit v1.2.3