[lmi-commits] [lmi] master 028b454 23/23: Find any root with only 64 function evaluations

[Greg Chicares]

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

    Find any root with only 64 function evaluations
    
    Statistics:
    
      evaluations  max   mean  sample-std-dev    commit        date
          7212      63   18.6       6.94       d6bd8029e  20210718T1630Z
          7501      75   19.3       7.36       d2dcbf860  20210723T2039Z
          7331      64   18.9       7.12       dc63e62cd  20210724T2119Z
          7332      65   18.9       7.13          this        today
    
    Picking 64 as the sprauchling limit has the charming consequence that
    no solve can take more than 128 evaluations, though it cost one
    evaluation.
    
    'Mang a' the ups an' doons o' life,
      There's lots o' things gae wrang,
    To disappoint us in our aims,
      And put us aff our gang;
    But could we see the goal we seek,
      That's yont the hill we speel,
    We'd think our failure to get there
      Was maybe just as weel.
    
    We've routh o' disappointments as
      We sprauchle up the brae,
    We dinna get the things we want,
      We lose the things we hae;
    An' yet 'mid a' that comes an' gangs,
      As time rows aff the reel,
    We're forced to this conclusion, that
      It's maybe just as weel.
    
    We deem events that loom afar
      Will prove to be a curse,
    Yet when they happen--o'd, hoo sune
      We think them the reverse;
    In fact, we never can foresee
      What may prove wae or weel,
    An' tho' we mourn owre blasted hopes,
      It's maybe just as weel.
    
    --Robert McLean Calder
---
 calendar_date.cpp |  1 +
 financial.hpp     |  1 +
 gpt_specamt.cpp   |  1 +
 ihs_avsolve.cpp   |  1 +
 solve.cpp         |  1 +
 zero.hpp          | 26 ++++++++++++++++++++++----
 zero_test.cpp     | 19 +++++++++++++------
 7 files changed, 40 insertions(+), 10 deletions(-)

diff --git a/calendar_date.cpp b/calendar_date.cpp
index 957880b..19a0e13 100644
--- a/calendar_date.cpp
+++ b/calendar_date.cpp
@@ -750,6 +750,7 @@ class birthdate_limit
             , 366 + a_priori_maximum_
             ,bias_
             ,0
+            ,64
             );
         LMI_ASSERT(root_is_valid == z.validity);
         int j = bourn_cast<int>(z.root);
diff --git a/financial.hpp b/financial.hpp
index 2788c75..bea9b39 100644
--- a/financial.hpp
+++ b/financial.hpp
@@ -107,6 +107,7 @@ class irr_helper
             ,1000.0     // Assumed upper bound.
             ,bias_lower // Return the final bound with the lower FV.
             ,decimals_
+            ,64
 //          ,os_trace
             );
         switch(z.validity)
diff --git a/gpt_specamt.cpp b/gpt_specamt.cpp
index 4760d64..d3727d8 100644
--- a/gpt_specamt.cpp
+++ b/gpt_specamt.cpp
@@ -165,6 +165,7 @@ currency gpt_specamt::CalculateSpecAmt
         ,999999999.99
         ,bias_higher
         ,z.round_min_specamt.decimals()
+        ,64
         );
 
     // Because it is implausible that the upper bound is too low,
diff --git a/ihs_avsolve.cpp b/ihs_avsolve.cpp
index 0ab7c85..6ec6b7d 100644
--- a/ihs_avsolve.cpp
+++ b/ihs_avsolve.cpp
@@ -486,6 +486,7 @@ currency AccountValue::Solve
         ,upper_bound
         ,bias
         ,decimals
+        ,64
         ,os_trace
         );
     currency const solution_cents = round_minutiae().c(solution.root);
diff --git a/solve.cpp b/solve.cpp
index 34e55ab..741e6dc 100644
--- a/solve.cpp
+++ b/solve.cpp
@@ -328,6 +328,7 @@ currency AccountValue::Solve()
         ,UpperBound
         ,Bias
         ,Decimals
+        ,64
         ,status()
         );
     currency const solution_cents = round_to_cents.c(solution.root);
diff --git a/zero.hpp b/zero.hpp
index d42690f..3d33e45 100644
--- a/zero.hpp
+++ b/zero.hpp
@@ -30,6 +30,7 @@
 #include "ssize_lmi.hpp"
 
 #include <cfloat>                       // DBL_EPSILON, DECIMAL_DIG
+#include <climits>                      // INT_MAX
 #include <cmath>                        // fabs(), isfinite(), isnan(), pow()
 #include <cstdint>                      // uint64_t
 #include <cstring>                      // memcpy()
@@ -64,6 +65,7 @@ enum root_impetus : char
     ,force_b_to_be_best_approximation = 'k'
     ,interpolate_linear               = 'L'
     ,interpolate_inverse_quadratic    = 'Q'
+    ,interpolate_guaranteed_64_evals  = 'G'
     ,dithering_near_root              = '0'
     ,secant_out_of_bounds             = '1'
     ,parabola_not_single_valued       = '2'
@@ -276,6 +278,13 @@ inline double binary64_midpoint(double d0, double d1)
 /// enforce it, and also handles the special case where both ordinates
 /// are zero.
 ///
+/// For Brent's method, the worst-case number of iterations is the
+/// square of the number required by naive bisection, so it may take
+/// an unreasonable amount of time for ill-conditioned problems. The
+/// optional 'sprauchling_limit' argument specifies the maximum number
+/// of evaluations to allow before switching to binary64 bisection,
+/// which is guaranteed to converge in 64 further evaluations.
+///
 /// Notes referred to in the source code
 ///
 /// Note 0. If one of the bounds is a zero, it is returned as soon as
@@ -351,8 +360,9 @@ root_type lmi_root
     ,double          bound0
     ,double          bound1
     ,double          tolerance
-    ,std::ostream&   os_trace  = null_stream()
-    ,root_bias       bias      = bias_none
+    ,int             sprauchling_limit = INT_MAX
+    ,std::ostream&   os_trace          = null_stream()
+    ,root_bias       bias              = bias_none
     )
 {
     int              n_iter  {0};
@@ -490,7 +500,13 @@ root_type lmi_root
                 ; // Do nothing.
                 }
             }
-        if(std::fabs(e) < tol)
+        if(sprauchling_limit < n_eval)
+            {
+            impetus = interpolate_guaranteed_64_evals;
+            n = binary64_midpoint(b, c); // "next" iterate
+            d = e = n - b;
+            }
+        else if(std::fabs(e) < tol)
             {
             impetus = dithering_near_root;
             d = e = n - b;
@@ -604,7 +620,8 @@ root_type decimal_root
     ,double          bound1
     ,root_bias       bias
     ,int             decimals
-    ,std::ostream&   os_trace = null_stream()
+    ,int             sprauchling_limit = INT_MAX
+    ,std::ostream&   os_trace          = null_stream()
     )
 {
     round_to<double> const round_dec {decimals, r_to_nearest};
@@ -631,6 +648,7 @@ root_type decimal_root
         ,round_dec(bound0)
         ,round_dec(bound1)
         ,0.5 * std::pow(10.0, -decimals)
+        ,sprauchling_limit
         ,os_trace
         ,bias
         );
diff --git a/zero_test.cpp b/zero_test.cpp
index d1fa18f..415095a 100644
--- a/zero_test.cpp
+++ b/zero_test.cpp
@@ -147,7 +147,7 @@ void test_a_function
 
     std::ostringstream os1;
     os1.precision(DECIMAL_DIG);
-    root_type r = lmi_root(f, bound0, bound1, tol, os1);
+    root_type r = lmi_root(f, bound0, bound1, tol, INT_MAX, os1);
     INVOKE_LMI_TEST(root_is_valid == r.validity, file, line);
     error = r.root - exact_root;
     INVOKE_LMI_TEST_RELATION(std::fabs(error),<=,maximum_error,file,line);
@@ -346,7 +346,7 @@ void test_fundamentals()
     // Same, with expatiation.
 
     std::ostringstream oss;
-    r = decimal_root(e_function, 0.5, 5.0, bias_none, 9, oss);
+    r = decimal_root(e_function, 0.5, 5.0, bias_none, 9, INT_MAX, oss);
     std::cout << oss.str() << std::endl;
 
     // Test use with function object.
@@ -753,7 +753,7 @@ void test_celebrated_equation()
     auto f = [](double x) {return x * x * x - 2.0 * x - 5.0;};
     std::ostringstream oss;
     oss.precision(17);
-    root_type r = decimal_root(f, -2.56, 2.56, bias_none, 21, oss);
+    root_type r = decimal_root(f, -2.56, 2.56, bias_none, 21, INT_MAX, oss);
     LMI_TEST(root_is_valid == r.validity);
     // This constant is from the cited blog; lmi yields this,
     // which agrees to sixteen significant digits:
@@ -827,7 +827,7 @@ void test_wikipedia_example()
 {
     auto f = [](double x) {return (x + 3.0) * (x - 1.0) * (x - 1.0);};
     std::ostringstream oss;
-    root_type r = decimal_root(f, -4.0, 4.0 / 3.0, bias_none, 15, oss);
+    root_type r = decimal_root(f, -4.0, 4.0 / 3.0, bias_none, 15, INT_MAX, 
oss);
     LMI_TEST(root_is_valid == r.validity);
     LMI_TEST(std::fabs(-3.0 - r.root) <= 1.0e-15);
     // Display this to investigate further:
@@ -965,7 +965,7 @@ void test_hodgepodge()
     // rather than a theoretical maximum. Perhaps they'll always
     // succeed, because floating-point behavior is determinate;
     // but small variations betoken no catastrophe.
-    LMI_TEST_RELATION(159,<=,r.n_eval); // weak
+    LMI_TEST_RELATION(156,<=,r.n_eval); // weak
     LMI_TEST_RELATION(r.n_eval,<=,166); // weak
 
     d = brent_zero(eq_2_1, -100.0, 100.0, 0.5);
@@ -1029,13 +1029,20 @@ void test_hodgepodge()
     LMI_TEST_EQUAL(55, r.n_eval); // weak
 
     std::ostringstream oss;
-    r = lmi_root(signum_offset, -1.0e300, 1.0e300, 5.0e-19, oss);
+    r = lmi_root(signum_offset, -1.0e300, 1.0e300, 5.0e-19, INT_MAX, oss);
     LMI_TEST(root_is_valid == r.validity);
     LMI_TEST(materially_equal(-1.0 / 3.0, r.root));
     LMI_TEST(r.n_eval <= 3023);
     LMI_TEST_EQUAL(1052, r.n_eval); // weak
     // Display this to investigate further:
 //  std::cout << oss.str() << std::endl;
+
+    // Find a root of this irksome function in 64 evaluations,
+    // to maximal precision, in an enormous interval.
+    r = lmi_root(signum_offset, -1.0e300, 1.0e300, 5.0e-19, 0);
+    LMI_TEST(root_is_valid == r.validity);
+    LMI_TEST(materially_equal(-1.0 / 3.0, r.root));
+    LMI_TEST_RELATION(64,<=,r.n_eval);
 }
 
 void test_former_rounding_problem()