Use Rouillier and Zimmerman's version of the Uspensky algorithm.

3 files changed, 153 insertions, 191 deletions

diff --git a/libdimension/dimension/geometry.h b/libdimension/dimension/geometry.h
index d2ad362..1248e01 100644
--- a/libdimension/dimension/geometry.h
+++ b/libdimension/dimension/geometry.h
@@ -124,12 +124,17 @@ dmnsn_degrees(double radians)
   return radians*45.0/atan(1.0);
 }
 
-/** Return the sign bit of a scalar. */
+/** Return the sign of a scalar. */
 DMNSN_INLINE int
-dmnsn_signbit(double n)
+dmnsn_sign(double n)
 {
-  /* Guarantee a 1 or 0 return, to allow testing two signs for equality */
-  return signbit(n) ? 1 : 0;
+  if (n > 0.0) {
+    return 1;
+  } else if (n < 0.0) {
+    return -1;
+  } else {
+    return 0;
+  }
 }
 
 /* Shorthand for vector/matrix construction */
diff --git a/libdimension/polynomial.c b/libdimension/polynomial.c
index 56f087d..86c5ebb 100644
--- a/libdimension/polynomial.c
+++ b/libdimension/polynomial.c
@@ -20,7 +20,8 @@
 
 /**
  * @file
- * Polynomials.
+ * Real root isolation algorithm based on work by Vincent, Uspensky, Collins and
+ * Akritas, Johnson, Krandick, and Rouillier and Zimmerman.
  */
 
 #include "dimension-internal.h"
@@ -41,7 +42,7 @@ dmnsn_real_degree(const double poly[], size_t degree)
 
 /** Divide each coefficient by the leading coefficient. */
 static inline void
-dmnsn_normalize_polynomial(double poly[], size_t degree)
+dmnsn_polynomial_normalize(double poly[], size_t degree)
 {
   for (size_t i = 0; i < degree; ++i) {
     poly[i] /= poly[degree];
@@ -50,188 +51,157 @@ dmnsn_normalize_polynomial(double poly[], size_t degree)
 }
 
 /** Eliminate trivial zero roots from \p poly[]. */
-static inline size_t
-dmnsn_eliminate_zero_roots(double poly[], size_t degree)
+static inline void
+dmnsn_eliminate_zero_roots(double **poly, size_t *degree)
 {
   size_t i;
-  for (i = 0; i <= degree && fabs(poly[i]) < dmnsn_epsilon; ++i);
-
-  if (i > 0) {
-    for (size_t j = 0; j + i <= degree; ++j) {
-      poly[j] = poly[j + i];
+  for (i = 0; i <= *degree; ++i) {
+    if (fabs((*poly)[i]) >= dmnsn_epsilon) {
+      break;
     }
   }
 
-  return degree - i;
+  *poly += i;
+  *degree -= i;
 }
 
-/** Returns the number of sign changes between coefficients of \p poly[]. */
-static inline size_t
-dmnsn_descartes_rule(const double poly[], size_t degree)
+/** Calculate a finite upper bound on the roots of a normalized polynomial. */
+static inline double
+dmnsn_root_bound(const double poly[], size_t degree)
 {
-  size_t changes = 0;
-  int lastsign = dmnsn_signbit(poly[degree]);
-  for (size_t i = degree; i-- > 0;) {
-    int sign = dmnsn_signbit(poly[i]);
-    if (fabs(poly[i]) >= dmnsn_epsilon && sign != lastsign) {
-      lastsign = sign;
-      ++changes;
-    }
+  double bound = fabs(poly[0]);
+  for (size_t i = 1; i < degree; ++i) {
+    bound = dmnsn_max(bound, fabs(poly[i]));
   }
-
-  return changes;
+  bound += 1.0;
+  return bound;
 }
 
-/** How many levels of Pascal's triangle to precompute. */
-#define DMNSN_NBINOM 11
-
-/** Pre-computed values of Pascal's triangle. */
-static const double dmnsn_pascals_triangle[DMNSN_NBINOM][DMNSN_NBINOM] = {
-  { 1.0 },
-  { 1.0, 1.0 },
-  { 1.0, 2.0, 1.0 },
-  { 1.0, 3.0, 3.0, 1.0 },
-  { 1.0, 4.0, 6.0, 4.0, 1.0 },
-  { 1.0, 5.0, 10.0, 10.0, 5.0, 1.0 },
-  { 1.0, 6.0, 15.0, 20.0, 15.0, 6.0, 1.0 },
-  { 1.0, 7.0, 21.0, 35.0, 35.0, 21.0, 7.0, 1.0 },
-  { 1.0, 8.0, 28.0, 56.0, 70.0, 56.0, 28.0, 8.0, 1.0 },
-  { 1.0, 9.0, 36.0, 84.0, 126.0, 126.0, 84.0, 36.0, 9.0, 1.0 },
-  { 1.0, 10.0, 45.0, 120.0, 210.0, 252.0, 210.0, 120.0, 45.0, 10.0, 1.0 }
-};
-
-/** Compute the (n k) binomial coefficient. */
-static inline double
-dmnsn_binom(size_t n, size_t k)
+/** Copy a polynomial. */
+static inline void
+dmnsn_polynomial_copy(double dest[], const double src[], size_t degree)
 {
-  if (k > n - k) {
-    k = n - k;
+  for (size_t i = 0; i <= degree; ++i) {
+    dest[i] = src[i];
   }
+}
 
-  double ret = 1.0;
-  for (size_t i = 0; i < k; ++i) {
-    ret *= n - i;
-    ret /= i + 1;
+/** Transform a polynomial by P'(x) = P(x + 1). */
+static inline void
+dmnsn_polynomial_translate(double poly[], size_t degree)
+{
+  for (size_t i = 0; i <= degree; ++i) {
+    for (size_t j = degree - i; j <= degree - 1; ++j) {
+      poly[j] += poly[j + 1];
+    }
   }
+}
 
-  return ret;
+/** Transform a polynomial by P'(x) = P(c*x). */
+static inline void
+dmnsn_polynomial_scale(double poly[], size_t degree, double c)
+{
+  double factor = c;
+  for (size_t i = 1; i <= degree; ++i) {
+    poly[i] *= factor;
+    factor *= c;
+  }
 }
 
-/** Find ranges that contain a single root, with Uspensky's algorithm. */
+/** Returns the result of Descartes' rule on x^degree * poly(1/(x + 1)). */
 static size_t
-dmnsn_uspensky_bounds(const double poly[], size_t degree, double bounds[][2],
-                      size_t max_roots)
+dmnsn_descartes_bound(const double poly[], size_t degree)
 {
-  size_t signchanges = dmnsn_descartes_rule(poly, degree);
-  if (signchanges == 0) {
-    return 0;
-  } else if (signchanges == 1) {
-    bounds[0][0] = +0.0;
-    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 (size_t i = degree; i-- > 0;) {
-      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;
-    }
+  /* Copy the polynomial so we can be destructive */
+  double p[degree + 1];
+  dmnsn_polynomial_copy(p, poly, degree);
 
-    /* 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 + 1], b[degree + 1];
-    if (dmnsn_likely(degree < DMNSN_NBINOM)) {
-      /* Use precomputed binomial coefficients if possible */
-      for (size_t i = 0; i <= degree; ++i) {
-        a[i] = poly[i];
-        b[i] = poly[degree - i];
-        for (size_t j = i + 1; j <= degree; ++j) {
-          double binom = dmnsn_pascals_triangle[j][i];
-          a[i] += binom*poly[j];
-          b[i] += binom*poly[degree - j];
-        }
-      }
-    } else {
-      for (size_t i = 0; i <= degree; ++i) {
-        a[i] = poly[i];
-        b[i] = poly[degree - i];
-        for (size_t j = i + 1; j <= degree; ++j) {
-          double binom = dmnsn_binom(j, i);
-          a[i] += binom*poly[j];
-          b[i] += binom*poly[degree - j];
-        }
-      }
+  /* Calculate poly(1/(1/x + 1)) which avoids reversal */
+  for (size_t i = 1; i <= degree; ++i) {
+    for (size_t j = i; j >= 1; --j) {
+      p[j] += p[j - 1];
     }
+  }
 
-    /* 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);
+  /* Find the number of sign changes in p[] */
+  size_t changes = 0;
+  int lastsign = dmnsn_sign(p[0]);
+  for (size_t i = 1; changes <= 1 && i <= degree; ++i) {
+    int sign = dmnsn_sign(p[i]);
+    if (sign != 0 && sign != lastsign) {
+      ++changes;
+      lastsign = sign;
     }
-    n += nb;
-    if (n == max_roots)
-      return n;
+  }
+
+  return changes;
+}
 
-    /* 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;
+/** Depth-first search of possible isolating intervals. */
+static size_t
+dmnsn_root_bounds_recursive(double poly[], size_t degree, double *c, double *k,
+                            double bounds[][2], size_t nbounds)
+{
+  size_t s = dmnsn_descartes_bound(poly, degree);
+  if (s >= 2) {
+    *c *= 2.0;
+    *k /= 2.0;
+    double currc = *c, currk = *k;
+
+    /* Get the left child */
+    dmnsn_polynomial_scale(poly, degree, 1.0/2.0);
+
+    size_t n = dmnsn_root_bounds_recursive(poly, degree, c, k, bounds, nbounds);
+    if (n == nbounds) {
+      return n;
     }
-    n += na;
 
+    /* Get the right child from the last tested polynomial */
+    dmnsn_polynomial_translate(poly, degree);
+    dmnsn_polynomial_scale(poly, degree, currk/(*k));
+
+    *c = currc + 1.0;
+    *k = currk;
+
+    bounds += n;
+    nbounds -= n;
+    n += dmnsn_root_bounds_recursive(poly, degree, c, k, bounds, nbounds);
     return n;
+  } else if (s == 1) {
+    bounds[0][0] = (*c)*(*k);
+    bounds[0][1] = (*c + 1.0)*(*k);
+    return 1;
+  } else {
+    return 0;
   }
 }
 
-/** Calculate a finite upper bound for the roots of the normalized polynomial
-    \p poly[]. */
-static inline double
-dmnsn_root_bound(const double poly[], size_t degree)
+/** Find ranges that contain a single root. */
+static size_t
+dmnsn_root_bounds(const double poly[], size_t degree, double bounds[][2],
+                  size_t nbounds)
 {
-  double bound = 0.0;
-  for (size_t i = 0; i < degree; ++i) {
-    bound = dmnsn_max(bound, fabs(poly[i]));
-  }
-  bound += 1.0;
-  return bound;
-}
+  /* Copy the polynomial so we can be destructive */
+  double p[degree + 1];
+  dmnsn_polynomial_copy(p, poly, degree);
 
-/** Improve a root with Newton's method. */
-static inline double
-dmnsn_improve_root(const double poly[], size_t degree, double x)
-{
-  double p;
-  do {
-    /* Calculate the value of the polynomial and its derivative at once */
-    p = poly[degree];
-    double dp = 0.0;
-    for (size_t i = degree; i-- > 0;) {
-      dp = dp*x + p;
-      p  = p*x  + poly[i];
-    }
+  /* Scale the roots to within (0, 1] */
+  double bound = dmnsn_root_bound(p, degree);
+  dmnsn_polynomial_scale(p, degree, bound);
 
-    x -= p/dp;
-  } while (fabs(p) >= fabs(x)*dmnsn_epsilon);
+  /* Bounding intervals are of the form (c*k, (c + 1)*k) */
+  double c = 0.0, k = 1.0;
 
-  return x;
+  /* Isolate the roots */
+  size_t n = dmnsn_root_bounds_recursive(p, degree, &c, &k, bounds, nbounds);
+
+  /* Scale the roots back to within (0, bound] */
+  for (size_t i = 0; i < n; ++i) {
+    bounds[i][0] *= bound;
+    bounds[i][1] *= bound;
+  }
+
+  return n;
 }
 
 /** Use the false position method to find a root in a range that contains
@@ -239,21 +209,16 @@ dmnsn_improve_root(const double poly[], size_t degree, double x)
 static inline double
 dmnsn_bisect_root(const double poly[], size_t degree, double min, double max)
 {
-  if (max - min < dmnsn_epsilon) {
-    /* Bounds are equal, use Newton's method */
-    return dmnsn_improve_root(poly, degree, min);
-  }
-
   double evmin = dmnsn_polynomial_evaluate(poly, degree, min);
   double evmax = dmnsn_polynomial_evaluate(poly, degree, max);
-  double initial_evmin = fabs(evmin), initial_evmax = fabs(evmax);
-  double mid = 0.0, evmid;
-  int lastsign = -1;
+  double evinitial = dmnsn_min(fabs(evmin), fabs(evmax));
+  double mid, evmid;
+  int lastsign = 0;
 
   do {
     mid = (min*evmax - max*evmin)/(evmax - evmin);
     evmid = dmnsn_polynomial_evaluate(poly, degree, mid);
-    int sign = dmnsn_signbit(evmid);
+    int sign = dmnsn_sign(evmid);
 
     if (mid < min) {
       /* This can happen due to numerical instability in the root bounding
@@ -267,7 +232,7 @@ dmnsn_bisect_root(const double poly[], size_t degree, double min, double max)
       evmin = evmax;
       max = mid;
       evmax = evmid;
-    } else if (sign == dmnsn_signbit(evmax)) {
+    } else if (sign == dmnsn_sign(evmax)) {
       max = mid;
       evmax = evmid;
       if (sign == lastsign) {
@@ -284,11 +249,11 @@ dmnsn_bisect_root(const double poly[], size_t degree, double min, double max)
     }
 
     lastsign = sign;
-  } while (fabs(evmid) >= fabs(mid)*dmnsn_epsilon
-           /* These conditions improve stability when one of the bounds is
-              close to a different root than we are trying to find */
-           || fabs(evmid) > initial_evmin
-           || fabs(evmid) > initial_evmax);
+  } while ((fabs(evmid) >= fabs(mid)*dmnsn_epsilon
+            /* This condition improves stability when one of the bounds is
+               close to a different root than we are trying to find */
+            || fabs(evmid) > evinitial)
+           && max - min >= fabs(mid)*dmnsn_epsilon);
 
   return mid;
 }
@@ -348,7 +313,7 @@ dmnsn_solve_cubic(double poly[4], double x[3])
   double q = poly[0] - poly[2]*(9.0*poly[1] - 2.0*b2)/27.0;
 
   double disc = 4.0*p*p*p + 27.0*q*q;
-  double bdiv3 = poly[2]/3;
+  double bdiv3 = poly[2]/3.0;
 
   if (disc < 0.0) {
     /* Three real roots -- this implies p < 0 */
@@ -392,15 +357,13 @@ dmnsn_solve_cubic(double poly[4], double x[3])
     return 0;
 }
 
-/* Uspensky's algorithm */
+/* Solve a polynomial */
 DMNSN_HOT size_t
 dmnsn_polynomial_solve(const double poly[], size_t degree, double x[])
 {
   /* Copy the polynomial so we can be destructive */
-  double p[degree + 1];
-  for (size_t i = 0; i <= degree; ++i) {
-    p[i] = poly[i];
-  }
+  double copy[degree + 1], *p = copy;
+  dmnsn_polynomial_copy(p, poly, degree);
 
   /* Index into x[] */
   size_t i = 0;
@@ -408,23 +371,17 @@ dmnsn_polynomial_solve(const double poly[], size_t degree, double x[])
   /* Account for leading zero coefficients */
   degree = dmnsn_real_degree(p, degree);
   /* Normalize the leading coefficient to 1.0 */
-  dmnsn_normalize_polynomial(p, degree);
+  dmnsn_polynomial_normalize(p, degree);
   /* Eliminate simple zero roots */
-  degree = dmnsn_eliminate_zero_roots(p, degree);
+  dmnsn_eliminate_zero_roots(&p, &degree);
 
   static const size_t max_algebraic = 3;
   if (degree > max_algebraic) {
     /* Find isolating intervals for (degree - max_algebraic) roots of p[] */
     double ranges[degree - max_algebraic][2];
-    size_t n = dmnsn_uspensky_bounds(p, degree, ranges, degree - max_algebraic);
-
-    /* Calculate a finite upper bound for the roots of p[] */
-    double absmax = dmnsn_root_bound(p, degree);
+    size_t n = dmnsn_root_bounds(p, degree, ranges, degree - max_algebraic);
 
     for (size_t j = 0; j < n; ++j) {
-      /* Replace large or infinite upper bounds with a finite one */
-      ranges[j][1] = dmnsn_min(ranges[j][1], absmax);
-
       /* Bisect within the found range */
       double r = dmnsn_bisect_root(p, degree, ranges[j][0], ranges[j][1]);
 
diff --git a/libdimension/tests/polynomial.c b/libdimension/tests/polynomial.c
index d66c975..b0d3b77 100644
--- a/libdimension/tests/polynomial.c
+++ b/libdimension/tests/polynomial.c
@@ -32,13 +32,13 @@ main(void)
   dmnsn_die_on_warnings(true);
 
   double poly[6], x[5];
-  /* 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;
-  poly[2] = -1215.2863928;
-  poly[1] = -646.270936;
-  poly[0] =  1447.8216;
+  /* poly[] = 2*(x + 1)*(x - 1.2345)*(x - 2.3456)*(x - 5)*(x - 100) */
+  poly[5] =  2.0;
+  poly[4] = -215.1602;
+  poly[3] =  1540.4520864;
+  poly[2] = -2430.5727856;
+  poly[1] = -1292.541872;
+  poly[0] =  2895.6432;
 
   size_t n = dmnsn_polynomial_solve(poly, 5, x);
   if (n != 4) {
@@ -52,7 +52,7 @@ main(void)
     double evmin = dmnsn_polynomial_evaluate(poly, 5, x[i] - dmnsn_epsilon);
     double ev    = dmnsn_polynomial_evaluate(poly, 5, x[i]);
     double evmax = dmnsn_polynomial_evaluate(poly, 5, x[i] + dmnsn_epsilon);
-    if (fabs(evmin) < ev || fabs(evmax) < ev) {
+    if (fabs(evmin) < fabs(ev) || fabs(evmax) < fabs(ev)) {
       fprintf(stderr, "--- Root %.15g is inaccurate! ---\n", x[i]);
       return EXIT_FAILURE;
     }