[lmi-commits] [lmi] master e41a14f 4/8: Improve instrumentation

[Greg Chicares]

branch: master
commit e41a14f0b21849d2745736f771fb7f423470fe03
Author: Gregory W. Chicares <gchicares@sbcglobal.net>
Commit: Gregory W. Chicares <gchicares@sbcglobal.net>

    Improve instrumentation
    
    When a trace is to be printed, it is useful to know what interpolation
    technique was used at each step.
    
    Brent resorts to bisection in two cases, described on pages 50 and 51 of
    _Algorithms for Minimization without Derivatives_ (ISBN 0-13-022335-2):
     - "Our main modification of Dekker's algorithm ensures that a bisection
       is done at least once in every 2 log₂(|b-c|/δ) consecutive steps."
     - "When inverse quadratic interpolation is used the interpolating
       parabola cannot be a good approximation to f unless it is single-
       valued between (f,f(b)) and (c,f(c)). Thus, it is natural to accept
       the point i if it lies between b and c, and up to three-quarters of
       the way from b to c: consider the limiting case where the
       interpolating parabola has a vertical tangent at c and f(b) = -f(c).
       Hence, we reject i if 2|p| ≥ 3|mq|."
    The trace distinguishes two cases as 'B' and 'b' respectively, which
    may correspond to the two cases above (the ALGOL original is obscure).
---
 zero.hpp | 44 ++++++++++++++++++++++++++++----------------
 1 file changed, 28 insertions(+), 16 deletions(-)

diff --git a/zero.hpp b/zero.hpp
index 28d7845..46f31e3 100644
--- a/zero.hpp
+++ b/zero.hpp
@@ -34,32 +34,42 @@
 #include <string>
 #include <utility>
 
-enum root_validity
-    {root_is_valid
-    ,root_not_bracketed
-    };
-
 enum root_bias
     {bias_none   // Return root z with f(z) closest to 0.0 .
     ,bias_lower  // Require f(z) <= 0.0 .
     ,bias_higher // Require  0.0 <= f(z).
     };
 
+enum interpolation_technique
+    {interpolate_initialization
+    ,interpolate_bisection0
+    ,interpolate_linear
+    ,interpolate_inverse_quadratic
+    ,interpolate_bisection1 // bisection when quadratic rejected
+    };
+
+enum root_validity
+    {root_is_valid
+    ,root_not_bracketed
+    };
+
 typedef std::pair<double,root_validity> root_type;
 
 namespace detail
 {
 inline void expatiate
-    (int         & number_of_iterations
-    ,std::ostream& os_trace
-    ,double        x
-    ,double        fx
+    (int                    & number_of_iterations
+    ,std::ostream           & os_trace
+    ,interpolation_technique  technique
+    ,double                   x
+    ,double                   fx
     )
 {
     if(os_trace.good())
         {
         os_trace
             << "iteration " << number_of_iterations++
+            << " "          << "IBLQb"[technique]
             << " iterand "  << x
             << " value "    << fx
             << std::endl
@@ -252,7 +262,8 @@ root_type decimal_root
 
     static double const epsilon = std::numeric_limits<double>::epsilon();
 
-    int number_of_iterations = 0;
+    int number_of_iterations          {0};
+    interpolation_technique technique {interpolate_initialization};
 
     double t = 0.5 * std::pow(10.0, -decimals);
 
@@ -262,14 +273,14 @@ root_type decimal_root
     double b = round_(bound1);
 
     double fa = static_cast<double>(f(a));
-    detail::expatiate(number_of_iterations, os_trace, a, fa);
+    detail::expatiate(number_of_iterations, os_trace, technique, a, fa);
     if(0.0 == fa) // Note 0.
         {
         return std::make_pair(a, root_is_valid);
         }
 
     double fb = static_cast<double>(f(b));
-    detail::expatiate(number_of_iterations, os_trace, b, fb);
+    detail::expatiate(number_of_iterations, os_trace, technique, b, fb);
     if(0.0 == fb) // Note 0 [bis].
         {
         return std::make_pair(b, root_is_valid);
@@ -322,7 +333,7 @@ root_type decimal_root
             }
         if(std::fabs(e) < tol || std::fabs(fa) <= std::fabs(fb))
             {
-            // Bisection.
+            technique = interpolate_bisection0;
             d = e = m;
             }
         else
@@ -331,13 +342,13 @@ root_type decimal_root
             double s = fb / fa; // Note 3.
             if(a == c)
                 {
-                // Linear interpolation.
+                technique = interpolate_linear;
                 p = 2.0 * m * s;
                 q = 1.0 - s;
                 }
             else
                 {
-                // Inverse quadratic interpolation.
+                technique = interpolate_inverse_quadratic;
                 q = fa / fc;
                 double r = fb / fc;
                 p = s * (2.0 * m * q * (q - r) - (b - a) * (r - 1.0));
@@ -362,6 +373,7 @@ root_type decimal_root
                 }
             else
                 {
+                technique = interpolate_bisection1;
                 d = e = m;
                 }
             }
@@ -392,7 +404,7 @@ root_type decimal_root
         else
             {
             fb = static_cast<double>(f(b));
-            detail::expatiate(number_of_iterations, os_trace, b, fb);
+            detail::expatiate(number_of_iterations, os_trace, technique, b, 
fb);
             }
         }
 }