[lmi-commits] [lmi] master dce3edc 08/11: Show variable shifts in optional trace

[Greg Chicares]

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

    Show variable shifts in optional trace
    
    Operations 'j' and 'k' are shifts of {a,b,c} that maintain invariants:
        force_b_and_c_to_bracket_root    // j
        force_b_to_be_best_approximation // k
    and displaying them explicitly makes the trace easier to follow.
    
    For example, the fourth evaluation of the objective function (resulting
    from inverse quadratic interpolation) in this trace excerpt:
    
    #eval       a        fa       b         fb       c       fc
      3 L       5  0.609438 3.80896   0.337356     0.5 -1.69315
      4 Q 3.80896  0.337356 2.57754  -0.053165     0.5 -1.69315
      4 j 3.80896  0.337356 2.57754  -0.053165 3.80896 0.337356
      5 L 2.57754 -0.053165 2.74518 0.00984766 3.80896 0.337356
    
    caused b and c to be on the same side of a root (because fb and fc have
    the same sign), but Brent's method requires them to bracket a root--so
    c's old value is discarded and replaced by a's value, before choosing
    the fifth candidate for evaluation. (It might be clearer to suppress the
    redundant "4" in "4 j", but the benefit doesn't justify the complexity.)
    
    Brent's method may be viewed as containing a state machine (whose inputs
    are internal). The state consists of:
      three abscissae {a,b,c} and their ordinates {fa,fb,fc}
    and there are about ten transitions that affect that state. In the
    trace, each element of the state is a column, each transition is a row,
    and row labels like '4 j' show the evolution.
---
 zero.hpp      | 4 ++++
 zero_test.cpp | 8 ++++++++
 2 files changed, 12 insertions(+)

diff --git a/zero.hpp b/zero.hpp
index 84bb473..48367e1 100644
--- a/zero.hpp
+++ b/zero.hpp
@@ -354,6 +354,7 @@ root_type lmi_root
             fc = fa;
             d = e = b - a;
             impetus = force_b_and_c_to_bracket_root;
+            expatiate();
             }
         // If 'c' is a closer approximant than 'b', then swap them,
         // discarding the old value of 'a'.
@@ -362,6 +363,7 @@ root_type lmi_root
              a =  b;  b =  c;  c =  a;
             fa = fb; fb = fc; fc = fa;
             impetus = force_b_to_be_best_approximation;
+            expatiate();
             }
         double tol = 2.0 * epsilon * std::fabs(b) + t;
         double m = 0.5 * (c - b);
@@ -576,12 +578,14 @@ double brent_zero
   interpolate:
     c = a; fc = fa; d = e = b - a;
     impetus = force_b_and_c_to_bracket_root;
+    expatiate();
   extrapolate:
     if(std::fabs(fc) < std::fabs(fb))
         {
          a =  b;  b =  c;  c =  a;
         fa = fb; fb = fc; fc = fa;
         impetus = force_b_to_be_best_approximation;
+        expatiate();
         }
     tol = 2.0 * DBL_EPSILON * std::fabs(b) + t;
     m = 0.5 * (c - b);
diff --git a/zero_test.cpp b/zero_test.cpp
index dbc03a4..9f44676 100644
--- a/zero_test.cpp
+++ b/zero_test.cpp
@@ -489,16 +489,24 @@ void test_celebrated_equation()
     std::string const verified = 1 + R"--cut-here--(
 #eval            a           fa            b           fb            c         
  fc
   2 i -2.5600000000000001 -16.657216000000002 2.5600000000000001 
6.6572160000000018            0            0
+  2 j -2.5600000000000001 -16.657216000000002 2.5600000000000001 
6.6572160000000018 -2.5600000000000001 -16.657216000000002
   3 L 2.5600000000000001 6.6572160000000018 1.0980323260716793 
-5.8721945393772152 -2.5600000000000001 -16.657216000000002
+  3 j 2.5600000000000001 6.6572160000000018 1.0980323260716793 
-5.8721945393772152 2.5600000000000001 6.6572160000000018
   4 L 1.0980323260716793 -5.8721945393772152 1.783216881610604 
-2.8960493667789873 2.5600000000000001 6.6572160000000018
   5 Q 1.783216881610604 -2.8960493667789873 2.2478393639958036 
1.8621631139566732 2.5600000000000001 6.6572160000000018
+  5 j 1.783216881610604 -2.8960493667789873 2.2478393639958036 
1.8621631139566732 1.783216881610604 -2.8960493667789873
   6 L 2.2478393639958036 1.8621631139566732 2.0660057758331045 
-0.3135140955237814 1.783216881610604 -2.8960493667789873
+  6 j 2.2478393639958036 1.8621631139566732 2.0660057758331045 
-0.3135140955237814 2.2478393639958036 1.8621631139566732
   7 L 2.0660057758331045 -0.3135140955237814 2.0922079131171945 
-0.026123094109737011 2.2478393639958036 1.8621631139566732
   8 Q 2.0922079131171945 -0.026123094109737011 2.0945566700001779 
5.7910818359374616e-05 2.2478393639958036 1.8621631139566732
+  8 j 2.0922079131171945 -0.026123094109737011 2.0945566700001779 
5.7910818359374616e-05 2.0922079131171945 -0.026123094109737011
   9 L 2.0945566700001779 5.7910818359374616e-05 2.0945514746903111 
-7.6478343657981895e-08 2.0922079131171945 -0.026123094109737011
+  9 j 2.0945566700001779 5.7910818359374616e-05 2.0945514746903111 
-7.6478343657981895e-08 2.0945566700001779 5.7910818359374616e-05
  10 L 2.0945514746903111 -7.6478343657981895e-08 2.0945514815423065 
-2.2382096176443156e-13 2.0945566700001779 5.7910818359374616e-05
  11 Q 2.0945514815423065 -2.2382096176443156e-13 2.0945514815423265 
-8.8817841970012523e-16 2.0945566700001779 5.7910818359374616e-05
  12 Q 2.0945514815423265 -8.8817841970012523e-16 2.0945514815423274 
9.7699626167013776e-15 2.0945566700001779 5.7910818359374616e-05
+ 12 j 2.0945514815423265 -8.8817841970012523e-16 2.0945514815423274 
9.7699626167013776e-15 2.0945514815423265 -8.8817841970012523e-16
+ 12 k 2.0945514815423274 9.7699626167013776e-15 2.0945514815423265 
-8.8817841970012523e-16 2.0945514815423274 9.7699626167013776e-15
 )--cut-here--";
 
     LMI_TEST_EQUAL(verified, oss.str());