From 54f9b5e58befe75c22976cc48c381679a5127251 Mon Sep 17 00:00:00 2001
From: Tavian Barnes <tavianator@gmail.com>
Date: Fri, 29 Oct 2010 00:40:12 -0400
Subject: Only return non-zero roots of polynomials.

---
 libdimension/dimension/polynomial.h |  4 ++--
 libdimension/polynomial.c           | 29 ++++++++++++++---------------
 2 files changed, 16 insertions(+), 17 deletions(-)

diff --git a/libdimension/dimension/polynomial.h b/libdimension/dimension/polynomial.h
index 52a44dd..f7d4517 100644
--- a/libdimension/dimension/polynomial.h
+++ b/libdimension/dimension/polynomial.h
@@ -41,8 +41,8 @@ dmnsn_evaluate_polynomial(double poly[], size_t degree, double x)
   return ret;
 }
 
-/* Stores the non-negative roots of poly[] in x[], and returns the number of
-   such roots that were stored */
+/* Stores the positive roots of poly[] in x[], and returns the number of such
+   roots that were stored */
 size_t dmnsn_solve_polynomial(double poly[], size_t degree, double x[]);
 
 /* Helper function to print a polynomial */
diff --git a/libdimension/polynomial.c b/libdimension/polynomial.c
index 6e8ce4e..5736a4d 100644
--- a/libdimension/polynomial.c
+++ b/libdimension/polynomial.c
@@ -27,7 +27,7 @@ static inline size_t
 dmnsn_solve_linear(double poly[2], double x[1])
 {
   x[0] = -poly[0]/poly[1];
-  return (x[0] >= 0.0) ? 1 : 0;
+  return (x[0] > 0.0) ? 1 : 0;
 }
 
 static inline size_t
@@ -38,29 +38,26 @@ dmnsn_solve_quadratic(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] > 0.0) ? ((x[1] > 0.0) ? 2 : 1) : 0;
   } else {
     return 0;
   }
 }
 
-static inline size_t
-dmnsn_zero_roots(double poly[], size_t *degree, double x[])
+/* Eliminate trivial zero roots from poly[] */
+static inline void
+dmnsn_eliminate_zero_roots(double poly[], size_t *degree)
 {
-  size_t i = 0;
-  while (i <= *degree && poly[i] == 0.0) {
-    x[i] = 0.0;
-    ++i;
-  }
+  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 <= *degree; ++j) {
+    for (size_t j = 0; j + i <= deg; ++j) {
       poly[j] = poly[j + i];
     }
   }
 
   *degree -= i;
-  return i;
 }
 
 /* Get the real degree of a polynomial, ignoring leading zeros */
@@ -97,6 +94,8 @@ dmnsn_bisect_root(double poly[], size_t degree, double min, double max)
       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 {
@@ -257,12 +256,12 @@ dmnsn_solve_polynomial(double poly[], size_t degree, double x[])
     p[i] = poly[i];
   }
 
-  size_t i = 0; /* Index into x[] */
-
-  /* Eliminate simple zero roots */
-  i += dmnsn_zero_roots(p, &degree, x + i);
   /* 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[] */
-- 
cgit v1.2.3