From 0e71596def36b0c0331401728ee0451c5410001b Mon Sep 17 00:00:00 2001
From: Sam Hocevar <sam@hocevar.net>
Date: Mon, 12 Jan 2015 18:38:52 +0000
Subject: [PATCH] lolremez: greatly improve root search times by using simple
 regula falsi.

---
 tools/lolremez/solver.cpp | 37 ++++++++++++++++++++++++-------------
 1 file changed, 24 insertions(+), 13 deletions(-)

diff --git a/tools/lolremez/solver.cpp b/tools/lolremez/solver.cpp
index 21a22182..0753277f 100644
--- a/tools/lolremez/solver.cpp
+++ b/tools/lolremez/solver.cpp
@@ -179,28 +179,39 @@ void RemezSolver::FindZeroes()
      * place as the absolute error! */
     for (int i = 0; i < m_order + 1; i++)
     {
-        struct { real value, error; } left, right, mid;
+        struct { real value, error; } a, b, c;
 
-        left.value = m_control[i];
-        left.error = EvalEstimate(left.value) - EvalFunc(left.value);
-        right.value = m_control[i + 1];
-        right.error = EvalEstimate(right.value) - EvalFunc(right.value);
+        a.value = m_control[i];
+        a.error = EvalEstimate(a.value) - EvalFunc(a.value);
+        b.value = m_control[i + 1];
+        b.error = EvalEstimate(b.value) - EvalFunc(b.value);
 
         static real limit = ldexp((real)1, -500);
         static real zero = (real)0;
-        while (fabs(left.value - right.value) > limit)
+        while (fabs(a.value - b.value) > limit)
         {
-            mid.value = (left.value + right.value) / 2;
-            mid.error = EvalEstimate(mid.value) - EvalFunc(mid.value);
+            /* Interpolate linearly instead of taking the midpoint, this
+             * leads to far better convergence (6:1 speedups). */
+            real t = abs(b.error) / (abs(a.error) + abs(b.error));
+            real newc = b.value + t * (a.value - b.value);
+
+            /* If the third point didn't change since last iteration,
+             * we may be at an inflection point. Use the midpoint to get
+             * out of this situation. */
+            c.value = newc == c.value ? (a.value + b.value) / 2 : newc;
+            c.error = EvalEstimate(c.value) - EvalFunc(c.value);
+
+            if (c.error == zero)
+                break;
 
-            if ((left.error <= zero && mid.error <= zero)
-                 || (left.error >= zero && mid.error >= zero))
-                left = mid;
+            if ((a.error < zero && c.error < zero)
+                 || (a.error > zero && c.error > zero))
+                a = c;
             else
-                right = mid;
+                b = c;
         }
 
-        m_zeroes[i] = mid.value;
+        m_zeroes[i] = c.value;
     }
 
     printf(" -:- times for zeroes: estimate %f func %f weight %f\n",