[lmi-commits] [lmi] odd/eraseme_error 1dbe0b4 02/10: Refactor

[Greg Chicares]

branch: odd/eraseme_error
commit 1dbe0b44e5bfb1aa1f33695a25add6d5cf1317fa
Author: Gregory W. Chicares <gchicares@sbcglobal.net>
Commit: Gregory W. Chicares <gchicares@sbcglobal.net>

    Refactor
    
    Brent's method uses bisection in four different circumstances, which
    will usefully be distinguished soon in the optional trace.
    
    This is a pure refactoring, modulo NaNs, for which
      x < y
    and
      y <= x
    are not equivalent.
---
 zero.hpp | 31 ++++++++++++++++++++++++-------
 1 file changed, 24 insertions(+), 7 deletions(-)

diff --git a/zero.hpp b/zero.hpp
index da1ae65..1836497 100644
--- a/zero.hpp
+++ b/zero.hpp
@@ -350,7 +350,12 @@ root_type lmi_root
                 ; // Do nothing.
                 }
             }
-        if(std::fabs(e) < tol || std::fabs(fa) <= std::fabs(fb))
+        if(std::fabs(e) < tol)
+            {
+            technique = interpolate_bisection0;
+            d = e = m;
+            }
+        else if(std::fabs(fa) <= std::fabs(fb))
             {
             technique = interpolate_bisection0;
             d = e = m;
@@ -383,18 +388,30 @@ root_type lmi_root
                 }
             s = e;
             e = d;
-            if
-                (   2.0 * p < 3.0 * m * q - std::fabs(tol * q)
-                &&  p < std::fabs(0.5 * s * q)
-                )
+            // Use the criteria in Brent's ALGOL, which differ
+            // slightly from their descriptions in his text.
+            //
+            // AfMWD says on page 51:
+            //   "we reject i [i.e., b + p/q] if 2|p| ≥ 3|mq|"
+            if(3.0 * m * q - std::fabs(tol * q) <= 2.0 * p)
                 {
-                d = p / q;
+                technique = interpolate_bisection1;
+                d = e = m;
                 }
-            else
+            // AfMWD says on page 50:
+            //   "Let e be the value of p/q at the step before the
+            //   last one."
+            // (That value is 's', both above and in the ALGOL.)
+            //   "If |e| < δ or |p/q| ≥ ½|e| then we do a bisection"
+            else if(std::fabs(0.5 * s * q) <= p)
                 {
                 technique = interpolate_bisection1;
                 d = e = m;
                 }
+            else
+                {
+                d = p / q;
+                }
             }
         a = b;
         fa = fb;