[lmi-commits] [lmi] odd/brent 0ad36f9 6/8: Display Brent and Chandrupatla conditions

[Greg Chicares]

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

    Display Brent and Chandrupatla conditions
    
    Brent's ALGOL unifies two conditions for rejecting an IQI iterate:
      if 2×p < 3×m×q - abs(tol×q) ∧ p < abs(0.5×s×q)
    of which the first is hypothesized to be almost the same as one of
    Chandrupatla's, the second being (Brent, page 50) "our main modification
    of Dekker's algorithm" which guarantees a reasonable rate of convergence
    even in extraordinarily ill-conditioned cases. The first is isolated as
    'cond_b', and Chandrupatla's is given as 'cond_c'; when Chandrupatla
    would reject an IQI iterate that Brent accepts, details are printed.
    
    This is the same f(x)=(x-1.7)^17 test as previously; it takes the same
    137 iterations whether the "#if" condition is 1 or 0, because when
    'cond_b' and 'cond_c' differ, it follows 'cond_b' as before--the only
    difference is the printing of details.
    
    (Puzzlingly, though, the number of decimal_root() iterations seems to
    have changed in this commit and the one preceding: going back two
      git switch --detach 3bd771343f5
    yields 129 iterations, while going back one
      git switch --detach ec32d2b4d2f
    yields 102. It would be a good idea to transplant Chandrupatla's
    criteria into 'brent_zero' as well, and test the outcome against the
          145 genuine brent
    result previously recorded and still reproduced now.)
    
    In this example at least, both ξ and φ are usually very close to one.
    At this time, the numerical stability of comparing ξ to φ^2 versus
    comparing their square roots is not considered.
    
    In all cases, it is the φ^2 < ξ criterion that fails, providing some
    support for the hypothesis that 'cond_b' corresponds to the other of
    Chandrupatla's two conditions.
    
    Let us consider the first rejection, on the third iteration:
      chandrupatla would reject the next iterate:
      1.98235290423467769294 in favor of
      0.999999999999844013665  based on these values:
    [reformatted]
      a, f(a): 2                           1.29140163000000323588e-09
      b, f(b): 1.99999999999968802733      1.29140162997717325777e-09
      c, f(c): 0                       -8272.40261886336338648
      0.999999999999844013665 xi
      1 phi
    Graphing points {a,b,c} makes rejection seem wise, as the slope of a
    line from b to a is over one million. Transformation to a unit parabola
    is quite straightforward, making it easy to validate 'xi' at least. It
    is also easy to verify that 'phi' is one, but that could be an accident.
---
 zero.hpp | 14 ++++++++------
 1 file changed, 8 insertions(+), 6 deletions(-)

diff --git a/zero.hpp b/zero.hpp
index aefde9c..4deb532 100644
--- a/zero.hpp
+++ b/zero.hpp
@@ -445,9 +445,10 @@ os_trace << " chandrupatla..." << std::endl;
                 d = p / q;
                 if(interpolate_inverse_quadratic == technique && !cond_c)
                     {
-#if 0
+#if 1
                     os_trace
-                        << "  chandrupatla would reject this:\n"
+                        << "\n"
+                        << "  chandrupatla would reject the next iterate:\n"
                         << "  " << b + p / q << " in favor of\n"
                         << "  " << b + m
                         << "  based on these values:\n"
@@ -463,10 +464,11 @@ os_trace << " chandrupatla..." << std::endl;
                         << "  " << phi << " phi\n"
                         << "  " << ((phi * phi) < xi) << " (phi * phi) < xi\n"
                         << "  " << (((1.0 - phi) * (1.0 - phi)) < (1.0 - xi)) 
<< " ((1.0 - phi) * (1.0 - phi)) < (1.0 - xi)\n"
-                        << std::endl;
-#endif // 0
-                    technique = interpolate_bisection1;
-                    d = e = m;
+                        << std::flush;
+#endif // 1
+// accept Chandrupatla's rejection:
+//                  technique = interpolate_bisection1;
+//                  d = e = m;
                     }
                 }
             else