[lmi-commits] [lmi] odd/brent e2b9e33 8/8: Let Chandrupatla reject IQI when Brent would accept it

[Greg Chicares]

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

    Let Chandrupatla reject IQI when Brent would accept it
    
    Number of iterations for "genuine brent":
      145 preceding commit
      101 this commit
    versus 51 for Chandrupatla's method as given by Scherer.
    
    Hypothesis: This is a pure improvement that Brent would have made if
    it had occurred to him.
    
    On page 53 of AfMWD, Brent says "in practical tests our algorithm has
    never taken more than 3(k+1) function evaluations", where k+1 is the
    number required for pure bisection. One might speculatively hope that
    this refinement would reduce that factor from 3 to 2.
    
    Chandrupatla's algorithm, as given by Scherer at least, is certainly
    faster in this particular (tortured) example. However, it is unknown
    whether it has been rigorously tested: Chandrupatla's work is behind
    a paywall; Scherer's code for Brent's method was imperfect in many
    respects, so it would be unwise to trust his code for Chandrupatla's
    method; and Brent's method has been widely used for decades. To achieve
    the same performance as bisection, Brent's method would need to track
    three distinct points at all times, instead of setting a=c whenever
    fb < fc; but that would seem to represent a drastic change that could
    risk impairing the guarantees he took such great pains to establish.
---
 zero.hpp | 8 ++++++--
 1 file changed, 6 insertions(+), 2 deletions(-)

diff --git a/zero.hpp b/zero.hpp
index 43f0964..30436c2 100644
--- a/zero.hpp
+++ b/zero.hpp
@@ -602,9 +602,13 @@ int j = 0;
                 {
                 if(interpolate_inverse_quadratic == technique && !cond_c)
                     {
-                    std::cout << "X ";
+                    technique = interpolate_bisection1;
+                    d = e = m;
+                    }
+                else
+                    {
+                    d = p / q;
                     }
-                d = p / q;
                 }
             else
                 {