From 8d8d44434bb876b8bf77ce2e77c2946cff0432c7 Mon Sep 17 00:00:00 2001
From: Tavian Barnes <tavianator@gmail.com>
Date: Thu, 28 Oct 2010 13:43:37 -0400
Subject: Use the false position method to bisect isolating intervals.

---
 libdimension/polynomial.c | 50 +++++++++++++++++++----------------------------
 1 file changed, 20 insertions(+), 30 deletions(-)

(limited to 'libdimension/polynomial.c')

diff --git a/libdimension/polynomial.c b/libdimension/polynomial.c
index e968832..d78361c 100644
--- a/libdimension/polynomial.c
+++ b/libdimension/polynomial.c
@@ -76,47 +76,38 @@ dmnsn_real_degree(double poly[], size_t degree)
   return 0;
 }
 
-/* Improve a root with Newton's method */
-static inline double
-dmnsn_improve_root(double poly[], size_t degree, double x)
-{
-  double dx;
-  do {
-    /* Calculate the value of the polynomial and its derivative at once */
-    double p = poly[degree], dp = 0.0;
-    for (ssize_t i = degree - 1; i >= 0; --i) {
-      dp = dp*x + p;
-      p  = p*x  + poly[i];
-    }
-
-    dx = p/dp;
-    x -= dx;
-  } while (fabs(dx/x) > dmnsn_epsilon);
-
-  return x;
-}
-
-/* Use the method of bisection to find a root in a range that contains exactly
-   one root, counting multiplicity */
+/* Use the false position method to find a root in a range that contains exactly
+   one root */
 static inline double
 dmnsn_bisect_root(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;
 
-  while (max - min > dmnsn_epsilon) {
-    double mid = (min + max)/2.0;
-    double evmid = dmnsn_evaluate_polynomial(poly, degree, mid);
-    if (dmnsn_signbit(evmid) == dmnsn_signbit(evmax)) {
+  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) {
+        evmin /= 2.0;
+      }
     } else {
       min = mid;
       evmin = evmid;
+      if (sign == lastsign) {
+        evmax /= 2.0;
+      }
     }
-  }
 
-  return (min + max)/2.0;
+    lastsign = sign;
+  } while (fabs(evmid) > fabs(mid)*dmnsn_epsilon);
+
+  return mid;
 }
 
 /* Use synthetic division to eliminate the root `r' from poly[] */
@@ -158,7 +149,7 @@ dmnsn_descartes_rule(double poly[], size_t degree)
 }
 
 /* Get the (n k) binomial coefficient */
-static double
+static inline double
 dmnsn_binom(size_t n, size_t k)
 {
   if (k > n - k) {
@@ -270,7 +261,6 @@ dmnsn_solve_polynomial(double poly[], size_t degree, double x[])
 
       /* Bisect within the found range */
       double r = dmnsn_bisect_root(p, degree, ranges[j][0], ranges[j][1]);
-      r = dmnsn_improve_root(p, degree, r);
 
       /* Use synthetic division to eliminate the root `r' */
       dmnsn_eliminate_root(p, &degree, r);
-- 
cgit v1.2.3