[lmi-commits] [lmi] odd/brent ec32d2b 5/8: Transplant Chandrupatla's condition into some Brent code

[Greg Chicares]

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

    Transplant Chandrupatla's condition into some Brent code
    
    [See previous commits on this branch for works cited.]
    
    Brent (page 51) says "it is natural to accept the [IQI point] if it lies
    between b and c, and up to three-quarters of the way from b to c".
    
    Scherer (page 90) restates that criterion, and continues "However, this
    does not include all possible cases where interpolation is applicable."
    Apparently he means that it doesn't reject all cases that it ought to.
    Then he says "Chandrupatla gave a more general discussion. The limiting
    condition is that the polynomial p(y) has a horizontal tangent at one of
    the boundaries", whereas Brent treated only the boundary at his 'c'.
    
    Scherer then derives Chandrupatla's condition by differentiating the
    Lagrange polynomial. It is easier to follow the reasoning here:
      https://github.com/SimpleArt/solver/wiki/Methods
    | ...first rescaling x and y so that the points lie in the unit square:
    |
    | x = (x2 - x1) / (x3 - x1)
    |
    | y = (f(x2) - f(x1)) / (f(x3) - f(x1))
    |
    | If x^2 < y and (1 - x)^2 < 1 - y, then inverse quadratic interpolation
    | is used.
    
    The graph following that quote makes everything clear. An equivalent
    graph is shown here:
      https://math.stackexchange.com/a/1552985
    
    Hypothesis: one of the "x^2 < y and (1 - x)^2 < 1 - y" criteria should
    be almost the same as Brent's, and the other should give a good reason
    for rejecting IQI when Brent would accept it.
---
 zero.hpp  | 10 +++++-----
 zero3.hpp |  6 +++---
 2 files changed, 8 insertions(+), 8 deletions(-)

diff --git a/zero.hpp b/zero.hpp
index 3a58f82..aefde9c 100644
--- a/zero.hpp
+++ b/zero.hpp
@@ -405,17 +405,17 @@ os_trace << " chandrupatla..." << std::endl;
             s = e;
             e = d;
 
+//          double xi  = (fb - fc) / (fa - fc);
+//          double phi = ( b -  c) / ( a -  c);
+// no...it's inverse interpolation
 //          double xi = (a - b) / (c - b);
 //          double phi = (fa - fb) / (fc - fb);
 // no...b is the latest, a is the previous, and c the oldest
 //          double xi = (b - a) / (c - a);
 //          double phi = (fb - fa) / (fc - fa);
-// no...try again...WAIT, ISN'T THIS THE CORRECT ONE?
-//          double xi  = ( b -  c) / ( a -  c);
-//          double phi = (fb - fc) / (fa - fc);
 // no...try again
-            double xi  = (fb - fc) / (fa - fc);
-            double phi = ( b -  c) / ( a -  c);
+            double xi  = ( b -  c) / ( a -  c);
+            double phi = (fb - fc) / (fa - fc);
 #if 0
             cond_c =
                    (phi * phi) < xi
diff --git a/zero3.hpp b/zero3.hpp
index 148d495..4ec6ce3 100644
--- a/zero3.hpp
+++ b/zero3.hpp
@@ -109,10 +109,10 @@ std::cout << ++j << " ch " << "IBLcQb"[technique] << ' ' 
<< t << ' ' << x << ' '
             {
             return xm;
             }
-// specimen values:
+// specimen observed values:
 // 0.0659374999999999961142 a
-// 0.0406250000000000013878 b this is one side of the bracket
-// 0.0912499999999999977796 c
+// 0.0406250000000000013878 b this is one end of the bracket
+// 0.0912499999999999977796 c so this is the other end of the bracket
 // 0.0447250871687336973292 true value
         double xi = (a - b) / (c - b);
         double phi = (fa - fb) / (fc - fb);