[lmi-commits] [lmi] odd/eraseme_error a3f38b5 03/10: Refactor instrumented reference implementation

[Greg Chicares]

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

    Refactor instrumented reference implementation
    
    This parallels the immediately preceding commit.
    
    Refactoring production and reference implementations separately can
    expose defects in the (identical) refactoring.
---
 zero.hpp | 31 ++++++++++++++++++++++++-------
 1 file changed, 24 insertions(+), 7 deletions(-)

diff --git a/zero.hpp b/zero.hpp
index 1836497..1a9f80b 100644
--- a/zero.hpp
+++ b/zero.hpp
@@ -544,7 +544,12 @@ double brent_zero
     if(tol < std::fabs(m) && 0.0 != fb)
         {
         // See if a bisection is forced.
-        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;
@@ -578,18 +583,30 @@ double brent_zero
                 }
             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;
         if(tol < std::fabs(d))