Numerical fixes for polynomial.c.

3 files changed, 115 insertions, 95 deletions

diff --git a/bench/libdimension/polynomial.c b/bench/libdimension/polynomial.c
index bcf9281..eaf7efe 100644
--- a/bench/libdimension/polynomial.c
+++ b/bench/libdimension/polynomial.c
@@ -27,29 +27,29 @@ main(void)
 #define NPOLY 5
   double p[NPOLY][NPOLY + 1], x[NPOLY];
 
-  // p[0][] = x - 0.5;
+  /* p[0][] = x - 0.5; */
   p[0][1] =  1.0;
   p[0][0] = -0.5;
 
-  // p[1][] = (x + 0.5)*(x - 0.5)
+  /* p[1][] = (x + 0.5)*(x - 0.5) */
   p[1][2] =  1.0;
   p[1][1] =  0.0;
   p[1][0] = -0.25;
 
-  // p[2][] = (x + 1)*(x - 1.2345)*(x - 100)
+  /* p[2][] = (x + 1)*(x - 1.2345)*(x - 100) */
   p[2][3] =  1.0;
   p[2][2] = -100.2345;
   p[2][1] =  22.2155;
   p[2][0] =  123.45;
 
-  // p[3][] = (x + 1)*(x - 1.2345)*(x - 5)*(x - 100)
+  /* p[3][] = (x + 1)*(x - 1.2345)*(x - 5)*(x - 100) */
   p[3][4] =  1.0;
   p[3][3] = -105.2345;
   p[3][2] =  523.388;
   p[3][1] =  12.3725;
   p[3][0] = -617.25;
 
-  // p[4][] = (x + 1)*(x - 1.2345)*(x - 2.3456)*(x - 5)*(x - 100)
+  /* p[4][] = (x + 1)*(x - 1.2345)*(x - 2.3456)*(x - 5)*(x - 100) */
   p[4][5] =  1.0;
   p[4][4] = -107.5801;
   p[4][3] =  770.2260432;
diff --git a/libdimension/polynomial.c b/libdimension/polynomial.c
index 14c33d0..6e2aed2 100644
--- a/libdimension/polynomial.c
+++ b/libdimension/polynomial.c
@@ -27,7 +27,7 @@ static inline size_t
 dmnsn_solve_linear(const double poly[2], double x[1])
 {
   x[0] = -poly[0]/poly[1];
-  return (x[0] > 0.0) ? 1 : 0;
+  return (x[0] >= dmnsn_epsilon) ? 1 : 0;
 }
 
 static inline size_t
@@ -38,34 +38,18 @@ dmnsn_solve_quadratic(const double poly[3], double x[2])
     double s = sqrt(disc);
     x[0] = (-poly[1] + s)/(2.0*poly[2]);
     x[1] = (-poly[1] - s)/(2.0*poly[2]);
-    return (x[0] > 0.0) ? ((x[1] > 0.0) ? 2 : 1) : 0;
+    return (x[0] >= dmnsn_epsilon) ? ((x[1] >= dmnsn_epsilon) ? 2 : 1) : 0;
   } else {
     return 0;
   }
 }
 
-/* Eliminate trivial zero roots from poly[] */
-static inline void
-dmnsn_eliminate_zero_roots(double poly[], size_t *degree)
-{
-  size_t i, deg = *degree;
-  for (i = 0; i <= deg && poly[i] == 0.0; ++i);
-
-  if (i > 0) {
-    for (size_t j = 0; j + i <= deg; ++j) {
-      poly[j] = poly[j + i];
-    }
-  }
-
-  *degree -= i;
-}
-
 /* Get the real degree of a polynomial, ignoring leading zeros */
 static inline size_t
 dmnsn_real_degree(const double poly[], size_t degree)
 {
   for (ssize_t i = degree; i >= 0; ++i) {
-    if (poly[i] != 0.0) {
+    if (fabs(poly[i]) >= dmnsn_epsilon) {
       return i;
     }
   }
@@ -73,56 +57,20 @@ dmnsn_real_degree(const double poly[], size_t degree)
   return 0;
 }
 
-/* Use the false position method to find a root in a range that contains exactly
-   one root */
-static inline double
-dmnsn_bisect_root(const double poly[], size_t degree, double min, double max)
+/* Eliminate trivial zero roots from poly[] */
+static inline void
+dmnsn_eliminate_zero_roots(double poly[], size_t *degree)
 {
-  double evmin = dmnsn_evaluate_polynomial(poly, degree, min);
-  double evmax = dmnsn_evaluate_polynomial(poly, degree, max);
-  double mid = 0.0, evmid;
-  int lastsign = -1;
-
-  if (max - min <= dmnsn_epsilon)
-    return min;
+  size_t i, deg = *degree;
+  for (i = 0; i <= deg && fabs(poly[i]) < dmnsn_epsilon; ++i);
 
-  do {
-    mid = (min*evmax - max*evmin)/(evmax - evmin);
-    evmid = dmnsn_evaluate_polynomial(poly, degree, mid);
-    int sign = dmnsn_signbit(evmid);
-    if (sign == dmnsn_signbit(evmax)) {
-      max = mid;
-      evmax = evmid;
-      if (sign == lastsign) {
-        /* Don't allow the algorithm to keep the same endpoint for three
-           iterations in a row; this ensures superlinear convergence */
-        evmin /= 2.0;
-      }
-    } else {
-      min = mid;
-      evmin = evmid;
-      if (sign == lastsign) {
-        evmax /= 2.0;
-      }
+  if (i > 0) {
+    for (size_t j = 0; j + i <= deg; ++j) {
+      poly[j] = poly[j + i];
     }
-
-    lastsign = sign;
-  } while (fabs(evmid) > fabs(mid)*dmnsn_epsilon);
-
-  return mid;
-}
-
-/* Use synthetic division to eliminate the root `r' from poly[] */
-static inline void
-dmnsn_eliminate_root(double poly[], size_t *degree, double r)
-{
-  double rem = poly[*degree];
-  for (ssize_t i = *degree - 1; i >= 0; --i) {
-    double temp = poly[i];
-    poly[i] = rem;
-    rem = temp + r*rem;
   }
-  --*degree;
+
+  *degree -= i;
 }
 
 /* Returns the number of sign changes between coefficients of `poly' */
@@ -132,7 +80,7 @@ dmnsn_descartes_rule(const double poly[], size_t degree)
   int lastsign = 0;
   size_t i;
   for (i = 0; i <= degree; ++i) {
-    if (poly[i] != 0.0) {
+    if (fabs(poly[i]) >= dmnsn_epsilon) {
       lastsign = dmnsn_signbit(poly[i]);
       break;
     }
@@ -141,7 +89,7 @@ dmnsn_descartes_rule(const double poly[], size_t degree)
   size_t changes = 0;
   for (++i; i <= degree; ++i) {
     int sign = dmnsn_signbit(poly[i]);
-    if (poly[i] != 0.0 && sign != lastsign) {
+    if (fabs(poly[i]) >= dmnsn_epsilon && sign != lastsign) {
       lastsign = sign;
       ++changes;
     }
@@ -199,11 +147,29 @@ dmnsn_uspensky_bounds(const double poly[], size_t degree, double bounds[][2],
     bounds[0][1] = INFINITY;
     return 1;
   } else {
+    /* Number of roots found so far */
     size_t n = 0;
 
+    /* First divide poly[] by (x - 1) to test for a root at x = 1 */
+    double pdiv1[degree], rem = poly[degree];
+    for (ssize_t i = degree - 1; i >= 0; --i) {
+      pdiv1[i] = rem;
+      rem += poly[i];
+    }
+    if (fabs(rem) < dmnsn_epsilon) {
+      bounds[n][0] = 1.0;
+      bounds[n][1] = 1.0;
+      ++n;
+      if (n == max_roots)
+        return n;
+
+      --degree;
+      poly = pdiv1;
+    }
+
     /* a[] is the expanded form of poly(x + 1), b[] is the expanded form of
        (x + 1)^degree * poly(1/(x + 1)) */
-    double a[degree], b[degree];
+    double a[degree + 1], b[degree + 1];
     for (size_t i = 0; i <= degree; ++i) {
       a[i] = poly[i];
       b[i] = poly[degree - i];
@@ -214,18 +180,10 @@ dmnsn_uspensky_bounds(const double poly[], size_t degree, double bounds[][2],
       }
     }
 
-    /* Test for a root at 1.0 */
-    if (a[0] == 0.0) {
-      bounds[n][0] = 1.0;
-      bounds[n][1] = 1.0;
-      ++n;
-      if (n == max_roots)
-        return n;
-    }
-
     /* Recursively test for roots in b[] */
     size_t nb = dmnsn_uspensky_bounds(b, degree, bounds + n, max_roots - n);
     for (size_t i = n; i < n + nb; ++i) {
+      /* Transform the found roots of b[] into roots of poly[] */
       double temp = bounds[i][0];
       bounds[i][0] = 1.0/(bounds[i][1] + 1.0);
       bounds[i][1] = 1.0/(temp + 1.0);
@@ -237,6 +195,7 @@ dmnsn_uspensky_bounds(const double poly[], size_t degree, double bounds[][2],
     /* Recursively test for roots in a[] */
     size_t na = dmnsn_uspensky_bounds(a, degree, bounds + n, max_roots - n);
     for (size_t i = n; i < n + na; ++i) {
+      /* Transform the found roots of a[] into roots of poly[] */
       bounds[i][0] += 1.0;
       bounds[i][1] += 1.0;
     }
@@ -246,35 +205,96 @@ dmnsn_uspensky_bounds(const double poly[], size_t degree, double bounds[][2],
   }
 }
 
-/* Modified Uspensky's algorithm */
+/* Calculate a finite upper bound for the roots of poly[] */
+static inline double
+dmnsn_root_bound(double poly[], size_t degree)
+{
+  double bound = 0.0;
+  for (size_t i = 0; i < degree; ++i) {
+    bound = dmnsn_max(bound, fabs(poly[i]));
+  }
+  bound /= fabs(poly[degree]);
+  bound += 1.0;
+  return bound;
+}
+
+/* Use the false position method to find a root in a range that contains exactly
+   one root */
+static inline double
+dmnsn_bisect_root(const double poly[], size_t degree, double min, double max)
+{
+  double evmin = dmnsn_evaluate_polynomial(poly, degree, min);
+  double evmax = dmnsn_evaluate_polynomial(poly, degree, max);
+  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);
+    int sign = dmnsn_signbit(evmid);
+    if (sign == dmnsn_signbit(evmax)) {
+      max = mid;
+      evmax = evmid;
+      if (sign == lastsign) {
+        /* Don't allow the algorithm to keep the same endpoint for three
+           iterations in a row; this ensures superlinear convergence */
+        evmin /= 2.0;
+      }
+    } else {
+      min = mid;
+      evmin = evmid;
+      if (sign == lastsign) {
+        evmax /= 2.0;
+      }
+    }
+
+    lastsign = sign;
+  } while (fabs(evmid) >= fabs(mid)*dmnsn_epsilon);
+
+  return mid;
+}
+
+/* Use synthetic division to eliminate the root `r' from poly[] */
+static inline void
+dmnsn_eliminate_root(double poly[], size_t *degree, double r)
+{
+  double rem = poly[*degree];
+  for (ssize_t i = *degree - 1; i >= 0; --i) {
+    double temp = poly[i];
+    poly[i] = rem;
+    rem = temp + r*rem;
+  }
+  --*degree;
+}
+
+/* Uspensky's algorithm */
 size_t
 dmnsn_solve_polynomial(const double poly[], size_t degree, double x[])
 {
   /* Copy the polynomial so we can be destructive */
   double p[degree + 1];
-  for (ssize_t i = degree; i >= 0; --i) {
+  for (size_t i = 0; i <= degree; ++i) {
     p[i] = poly[i];
   }
 
+  /* Index into x[] */
+  size_t i = 0;
+
   /* Account for leading zero coefficients */
   degree = dmnsn_real_degree(p, degree);
   /* Eliminate simple zero roots */
   dmnsn_eliminate_zero_roots(p, &degree);
 
-  size_t i = 0; /* Index into x[] */
-
   if (degree >= 3) {
-    /* Find isolating intervals for degree - 2 roots of p[] */
+    /* Find isolating intervals for (degree - 2) roots of p[] */
     double ranges[degree - 2][2];
     size_t n = dmnsn_uspensky_bounds(p, degree, ranges, degree - 2);
 
-    /* Calculate an finite upper bound for the roots of p[] */
-    double absmax = 0.0;
-    for (size_t j = 0; j < degree; ++j) {
-      absmax = dmnsn_max(absmax, fabs(p[j]));
-    }
-    absmax /= fabs(p[degree]);
-    absmax += 1.0;
+    /* Calculate a finite upper bound for the roots of p[] */
+    double absmax = dmnsn_root_bound(p, degree);
 
     for (size_t j = 0; j < n; ++j) {
       /* Replace large or infinite upper bounds with a finite one */
diff --git a/tests/libdimension/polynomial.c b/tests/libdimension/polynomial.c
index 230a1a7..8c6b4da 100644
--- a/tests/libdimension/polynomial.c
+++ b/tests/libdimension/polynomial.c
@@ -29,7 +29,7 @@ int
 main()
 {
   double poly[6], x[5];
-  // poly[] = (x + 1)*(x - 1.2345)*(x - 2.3456)*(x - 5)*(x - 100)
+  /* poly[] = (x + 1)*(x - 1.2345)*(x - 2.3456)*(x - 5)*(x - 100) */
   poly[5] =  1.0;
   poly[4] = -107.5801;
   poly[3] =  770.2260432;