summaryrefslogtreecommitdiffstats
diff options
context:
space:
mode:
-rw-r--r--libdimension/polynomial.c8
-rw-r--r--libdimension/tests/polynomial.c24
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);
+}