[lmi-commits] [lmi] odd/eraseme_error 2870f56 05/10: Reenumerate root-finding techniques

[Greg Chicares]

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

    Reenumerate root-finding techniques
    
    Single-character codes tersely elucidate traces. Options considered:
    
      X Y Z
      -----
      i 0 E_valuate  evaluate_bounds
      j 1 C_onstrain force_b_and_c_to_bracket_root
      k 2 B_est      force_b_to_be_best_approximation
      L L S_ecant    interpolate_linear
      Q Q I_QI       interpolate_inverse_quadratic
      0 d D_ither    dithering_near_root
      1 o O_OB       secant_out_of_bounds
      2 p P_arabola  parabola_not_single_valued
      3 g G_uarantee guarantee_linear_convergence
    
    Z: Letters alone could be chosen to suggest mnemonics--too precious.
    Y: Lowercase mnemonics for bisection remain obscure.
    X: Lowercase=initialization; capital=higher-order; digit=bisection
    
    Option X seems best. Symbols such as "-\|/_+" or unicode exotica
    were considered, but found ambiguous (which of "/\" is a secant?).
---
 zero.hpp      | 46 +++++++++++++++++++++++++++-------------------
 zero_test.cpp |  2 +-
 2 files changed, 28 insertions(+), 20 deletions(-)

diff --git a/zero.hpp b/zero.hpp
index f22db1c..a45284f 100644
--- a/zero.hpp
+++ b/zero.hpp
@@ -41,11 +41,15 @@ enum root_bias
     };
 
 enum root_technique
-    {interpolate_initialization
-    ,interpolate_bisection0
+    {evaluate_bounds
+    ,force_b_and_c_to_bracket_root
+    ,force_b_to_be_best_approximation
     ,interpolate_linear
     ,interpolate_inverse_quadratic
-    ,interpolate_bisection1 // bisection when quadratic rejected
+    ,dithering_near_root
+    ,secant_out_of_bounds
+    ,parabola_not_single_valued
+    ,guarantee_linear_convergence
     };
 
 enum root_validity
@@ -248,7 +252,7 @@ root_type lmi_root
     constexpr double epsilon   {std::numeric_limits<double>::epsilon()};
 
     int              n_iter    {0};
-    root_technique   technique {interpolate_initialization};
+    root_technique   technique {evaluate_bounds};
 
     os_trace
         << "#eval"
@@ -270,7 +274,7 @@ root_type lmi_root
         {
         os_trace
             <<        std::setw(3)  << n_iter
-            << ' '                  << "IBLQb"[technique]
+            << ' '                  << "ijkLQ0123"[technique]
             << ' ' << std::setw(12) << a << ' ' << std::setw(12) << fa
             << ' ' << std::setw(12) << b << ' ' << std::setw(12) << fb
             << ' ' << std::setw(12) << c << ' ' << std::setw(12) << fc
@@ -321,6 +325,7 @@ root_type lmi_root
             c = a;
             fc = fa;
             d = e = b - a;
+            technique = force_b_and_c_to_bracket_root;
             }
         // If 'c' is a closer approximant than 'b', then swap them,
         // discarding the old value of 'a'.
@@ -328,6 +333,7 @@ root_type lmi_root
             {
              a =  b;  b =  c;  c =  a;
             fa = fb; fb = fc; fc = fa;
+            technique = force_b_to_be_best_approximation;
             }
         double tol = 2.0 * epsilon * std::fabs(b) + t;
         double m = 0.5 * (c - b);
@@ -352,12 +358,12 @@ root_type lmi_root
             }
         if(std::fabs(e) < tol)
             {
-            technique = interpolate_bisection0;
+            technique = dithering_near_root;
             d = e = m;
             }
         else if(std::fabs(fa) <= std::fabs(fb))
             {
-            technique = interpolate_bisection0;
+            technique = secant_out_of_bounds;
             d = e = m;
             }
         else
@@ -395,7 +401,7 @@ root_type lmi_root
             //   "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)
                 {
-                technique = interpolate_bisection1;
+                technique = parabola_not_single_valued;
                 d = e = m;
                 }
             // AfMWD says on page 50:
@@ -405,7 +411,7 @@ root_type lmi_root
             //   "If |e| < δ or |p/q| ≥ ½|e| then we do a bisection"
             else if(std::fabs(0.5 * s * q) <= p)
                 {
-                technique = interpolate_bisection1;
+                technique = guarantee_linear_convergence;
                 d = e = m;
                 }
             else
@@ -501,7 +507,7 @@ double brent_zero
     // that f(a) and f(b) have different signs.
 
     int            n_iter    {0};
-    root_technique technique {interpolate_initialization};
+    root_technique technique {evaluate_bounds};
 
     os_trace
         << "#eval"
@@ -517,7 +523,7 @@ double brent_zero
         {
         os_trace
             <<        std::setw(3)  << n_iter
-            << ' '                  << "IBLQb"[technique]
+            << ' '                  << "ijkLQ0123"[technique]
             << ' ' << std::setw(12) << a << ' ' << std::setw(12) << fa
             << ' ' << std::setw(12) << b << ' ' << std::setw(12) << fb
             << ' ' << std::setw(12) << c << ' ' << std::setw(12) << fc
@@ -533,12 +539,14 @@ double brent_zero
     expatiate();
   interpolate:
     c = a; fc = fa; d = e = b - a;
+    technique = force_b_and_c_to_bracket_root;
   extrapolate:
     if(std::fabs(fc) < std::fabs(fb))
-       {
-        a =  b;  b =  c;  c =  a;
-       fa = fb; fb = fc; fc = fa;
-       }
+        {
+         a =  b;  b =  c;  c =  a;
+        fa = fb; fb = fc; fc = fa;
+        technique = force_b_to_be_best_approximation;
+        }
     tol = 2.0 * DBL_EPSILON * std::fabs(b) + t;
     m = 0.5 * (c - b);
     if(tol < std::fabs(m) && 0.0 != fb)
@@ -546,12 +554,12 @@ double brent_zero
         // See if a bisection is forced.
         if(std::fabs(e) < tol)
             {
-            technique = interpolate_bisection0;
+            technique = dithering_near_root;
             d = e = m;
             }
         else if(std::fabs(fa) <= std::fabs(fb))
             {
-            technique = interpolate_bisection0;
+            technique = secant_out_of_bounds;
             d = e = m;
             }
         else
@@ -590,7 +598,7 @@ double brent_zero
             //   "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)
                 {
-                technique = interpolate_bisection1;
+                technique = parabola_not_single_valued;
                 d = e = m;
                 }
             // AfMWD says on page 50:
@@ -600,7 +608,7 @@ double brent_zero
             //   "If |e| < δ or |p/q| ≥ ½|e| then we do a bisection"
             else if(std::fabs(0.5 * s * q) <= p)
                 {
-                technique = interpolate_bisection1;
+                technique = guarantee_linear_convergence;
                 d = e = m;
                 }
             else
diff --git a/zero_test.cpp b/zero_test.cpp
index 5bb93f8..603a850 100644
--- a/zero_test.cpp
+++ b/zero_test.cpp
@@ -450,7 +450,7 @@ void test_celebrated_equation()
     // "1 + " skips the newline:
     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 i -2.5600000000000001 -16.657216000000002 2.5600000000000001 
6.6572160000000018            0            0
   3 L 2.5600000000000001 6.6572160000000018 1.0980323260716793 
-5.8721945393772152 -2.5600000000000001 -16.657216000000002
   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