Fix the polynomial solver when the bounds are exact.

2 files changed, 31 insertions, 1 deletions

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) {
diff --git a/libdimension/tests/polynomial.c b/libdimension/tests/polynomial.c
index 208edf9..335fd27 100644
--- a/libdimension/tests/polynomial.c
+++ b/libdimension/tests/polynomial.c
@@ -23,6 +23,7 @@
  */
 
 #include "tests.h"
+#include "../polynomial.c"
 
 /* poly[] = 2*(x + 1)*(x - 1.2345)*(x - 2.3456)*(x - 5)*(x - 100) */
 static const double poly[6] = {
@@ -66,3 +67,26 @@ DMNSN_TEST(polynomial, accurate_roots)
     ck_assert(fabs(ev) < DMNSN_CLOSE_ENOUGH);
   }
 }
+
+/* repeated_root[] = (x - 1)^6 */
+static const double repeated_root[7] = {
+  [6] =  1.0,
+  [5] = -6.0,
+  [4] =  15.0,
+  [3] = -20.0,
+  [2] =  15.0,
+  [1] = -6.0,
+  [0] =  1.0,
+};
+
+DMNSN_TEST(stability, equal_bounds)
+{
+  double root = dmnsn_bisect_root(repeated_root, 6, 1.0, 1.0);
+  ck_assert_msg(root == 1.0, "root == %.17g", root);
+}
+
+DMNSN_TEST(stability, equal_values_at_bounds)
+{
+  double root = dmnsn_bisect_root(repeated_root, 6, 0.9, 1.1);
+  ck_assert_msg(fabs(root - 1.0) < dmnsn_epsilon, "root == %.17g", root);
+}