Fix a bug in cubic polynomial solver, and add more tests.

1 files changed, 160 insertions, 32 deletions

diff --git a/libdimension/tests/polynomial.c b/libdimension/tests/polynomial.c
index 335fd27..b734cc2 100644
--- a/libdimension/tests/polynomial.c
+++ b/libdimension/tests/polynomial.c
@@ -24,59 +24,187 @@
 
 #include "tests.h"
 #include "../polynomial.c"
+#include <stdarg.h>
 
-/* poly[] = 2*(x + 1)*(x - 1.2345)*(x - 2.3456)*(x - 5)*(x - 100) */
-static const double poly[6] = {
-  [5] =  2.0,
-  [4] = -215.1602,
-  [3] =  1540.4520864,
-  [2] = -2430.5727856,
-  [1] = -1292.541872,
-  [0] =  2895.6432,
-};
+#define DMNSN_CLOSE_ENOUGH 1.0e-6
 
-static double roots[5];
-static size_t nroots;
+static void
+dmnsn_assert_roots(const double poly[], size_t degree, size_t nroots_ex, ...)
+{
+  double roots[degree - 1];
+  size_t nroots = dmnsn_polynomial_solve(poly, degree, roots);
+  ck_assert_int_eq(nroots, nroots_ex);
+
+  va_list ap;
+  va_start(ap, nroots_ex);
+  for (size_t i = 0; i < nroots; ++i) {
+    double root_ex = va_arg(ap, double);
+    bool found = false;
+    for (size_t j = 0; j < nroots; ++j) {
+      if (fabs(root_ex - roots[j]) >= dmnsn_epsilon) {
+        continue;
+      }
+      found = true;
+
+      double evroot = dmnsn_polynomial_evaluate(poly, degree, roots[j]);
+      ck_assert(fabs(evroot) < DMNSN_CLOSE_ENOUGH);
+
+      double evmin = dmnsn_polynomial_evaluate(poly, degree, roots[j] - dmnsn_epsilon);
+      double evmax = dmnsn_polynomial_evaluate(poly, degree, roots[j] + dmnsn_epsilon);
+      ck_assert(fabs(evroot) <= fabs(evmin) && fabs(evroot) <= fabs(evmax));
+
+      break;
+    }
+
+    if (!found) {
+      for (size_t j = 0; j < nroots; ++j) {
+        fprintf(stderr, "roots[%zu] == %.17g\n", j, roots[j]);
+      }
+      ck_abort_msg("Expected root %.17g not found", root_ex);
+    }
+  }
+  va_end(ap);
+}
 
-#define DMNSN_CLOSE_ENOUGH 1.0e-6
 
-DMNSN_TEST_SETUP(polynomial)
+DMNSN_TEST(linear, no_positive_roots)
 {
-  nroots = dmnsn_polynomial_solve(poly, 5, roots);
+  /* poly[] = x + 1 */
+  static const double poly[] = {
+    [1] = 1.0,
+    [0] = 1.0,
+  };
+  dmnsn_assert_roots(poly, 1, 0);
 }
 
-DMNSN_TEST(polynomial, finds_positive_roots)
+DMNSN_TEST(linear, one_root)
 {
-  ck_assert_int_eq(nroots, 4);
+  /* poly[] = x - 1 */
+  static const double poly[] = {
+    [1] =  1.0,
+    [0] = -1.0,
+  };
+  dmnsn_assert_roots(poly, 1, 1, 1.0);
 }
 
-DMNSN_TEST(polynomial, local_min_roots)
+
+DMNSN_TEST(quadratic, no_roots)
 {
-  for (size_t i = 0; i < nroots; ++i) {
-    double evmin = dmnsn_polynomial_evaluate(poly, 5, roots[i] - dmnsn_epsilon);
-    double ev    = dmnsn_polynomial_evaluate(poly, 5, roots[i]);
-    double evmax = dmnsn_polynomial_evaluate(poly, 5, roots[i] + dmnsn_epsilon);
-    ck_assert(fabs(ev) < fabs(evmin) && fabs(ev) < fabs(evmax));
-  }
+  /* poly[] = x^2 + 1 */
+  static const double poly[] = {
+    [2] = 1.0,
+    [1] = 0.0,
+    [0] = 1.0,
+  };
+  dmnsn_assert_roots(poly, 2, 0);
 }
 
-DMNSN_TEST(polynomial, accurate_roots)
+DMNSN_TEST(quadratic, no_positive_roots)
 {
-  for (size_t i = 0; i < nroots; ++i) {
-    double ev = dmnsn_polynomial_evaluate(poly, 5, roots[i]);
-    ck_assert(fabs(ev) < DMNSN_CLOSE_ENOUGH);
-  }
+  /* poly[] = (x + 1)^2 */
+  static const double poly[] = {
+    [2] = 1.0,
+    [1] = 2.0,
+    [0] = 1.0,
+  };
+  dmnsn_assert_roots(poly, 2, 0);
+}
+
+DMNSN_TEST(quadratic, one_positive_root)
+{
+  /* poly[] = (x + 1)*(x - 1) */
+  static const double poly[] = {
+    [2] =  1.0,
+    [1] =  0.0,
+    [0] = -1.0,
+  };
+  dmnsn_assert_roots(poly, 2, 1, 1.0);
+}
+
+DMNSN_TEST(quadratic, two_roots)
+{
+  /* poly[] = (x - 1.2345)*(x - 2.3456) */
+  static const double poly[] = {
+    [2] =  1.0,
+    [1] = -3.5801,
+    [0] =  2.8956432,
+  };
+  dmnsn_assert_roots(poly, 2, 2, 1.2345, 2.3456);
+}
+
+
+DMNSN_TEST(cubic, no_positive_roots)
+{
+  /* poly[] = x^3 + 1 */
+  static const double poly[] = {
+    [3] = 1.0,
+    [2] = 0.0,
+    [1] = 0.0,
+    [0] = 1.0,
+  };
+  dmnsn_assert_roots(poly, 3, 0);
+}
+
+DMNSN_TEST(cubic, one_root)
+{
+  /* poly[] = x^3 - 1 */
+  static const double poly[] = {
+    [3] =  1.0,
+    [2] =  0.0,
+    [1] =  0.0,
+    [0] = -1.0,
+  };
+  dmnsn_assert_roots(poly, 3, 1, 1.0);
+}
+
+DMNSN_TEST(cubic, two_roots)
+{
+  /* poly[] = (x - 1)*(x - 4)^2 **/
+  static const double poly[] = {
+    [3] =   1.0,
+    [2] =  -9.0,
+    [1] =  24.0,
+    [0] = -16.0,
+  };
+  dmnsn_assert_roots(poly, 3, 2, 1.0, 4.0);
+}
+
+DMNSN_TEST(cubic, three_roots)
+{
+  /* poly[] = (x - 1.2345)*(x - 2.3456)*(x - 100) */
+  static const double poly[] = {
+    [3] =    1.0,
+    [2] = -103.5801,
+    [1] =  360.9056432,
+    [0] = -289.56432,
+  };
+  dmnsn_assert_roots(poly, 3, 3, 1.2345, 2.3456, 100.0);
+}
+
+
+DMNSN_TEST(quintic, four_roots)
+{
+  /* poly[] = 2*(x + 1)*(x - 1.2345)*(x - 2.3456)*(x - 5)*(x - 100) */
+  static const double poly[] = {
+    [5] =     2.0,
+    [4] =  -215.1602,
+    [3] =  1540.4520864,
+    [2] = -2430.5727856,
+    [1] = -1292.541872,
+    [0] =  2895.6432,
+  };
+  dmnsn_assert_roots(poly, 5, 4, 1.2345, 2.3456, 5.0, 100.0);
 }
 
 /* repeated_root[] = (x - 1)^6 */
 static const double repeated_root[7] = {
-  [6] =  1.0,
-  [5] = -6.0,
+  [6] =   1.0,
+  [5] =  -6.0,
   [4] =  15.0,
   [3] = -20.0,
   [2] =  15.0,
-  [1] = -6.0,
-  [0] =  1.0,
+  [1] =  -6.0,
+  [0] =   1.0,
 };
 
 DMNSN_TEST(stability, equal_bounds)