[lmi-commits] [lmi] master dc63e62 16/23: Cache evaluations for rounded decimal solves

[Greg Chicares]

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

    Cache evaluations for rounded decimal solves
    
    Avoiding superfluous evaluations (of iterates that round to the same
    value) regains most of the performance recently lost:
    
      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          this        today
---
 zero.hpp      | 23 ++++++++++++++++++++++-
 zero_test.cpp | 15 ++++++++-------
 2 files changed, 30 insertions(+), 8 deletions(-)

diff --git a/zero.hpp b/zero.hpp
index 6349598..1b89b0f 100644
--- a/zero.hpp
+++ b/zero.hpp
@@ -27,6 +27,7 @@
 #include "math_functions.hpp"           // signum()
 #include "null_stream.hpp"
 #include "round_to.hpp"
+#include "ssize_lmi.hpp"
 
 #include <cfloat>                       // DBL_EPSILON, DECIMAL_DIG
 #include <cmath>                        // fabs(), isfinite(), isnan(), pow()
@@ -36,6 +37,7 @@
 #include <limits>
 #include <numeric>                      // midpoint()
 #include <ostream>
+#include <unordered_map>
 #include <utility>                      // forward()
 
 enum root_bias
@@ -608,7 +610,23 @@ root_type decimal_root
     )
 {
     round_to<double> const round_dec {decimals, r_to_nearest};
-    auto fr = [&](double x) {return f(round_dec(x));};
+
+    std::unordered_map<double,double> m;
+
+    auto fr = [&](double x) // f(), rounded
+        {
+        double const r = round_dec(x);
+        auto const i = m.find(r);
+        if(m.end() != i)
+            {
+            os_trace << "Superfluous evaluation avoided" << std::endl;
+            return i->second;
+            }
+        else
+            {
+            return m[r] = static_cast<double>(f(r));
+            }
+        };
 
     auto z = lmi_root
         (fr
@@ -619,6 +637,9 @@ root_type decimal_root
         ,bias
         );
     z.root = round_dec(z.root);
+    os_trace << " function evaluations: " << z.n_eval;
+    z.n_eval = lmi::ssize(m);
+    os_trace << " " << z.n_eval << " nominal, actual" << std::endl;
     os_trace << " return value: " << z.root << " (rounded)" << std::endl;
     return z;
 }
diff --git a/zero_test.cpp b/zero_test.cpp
index 2b8025c..28696cc 100644
--- a/zero_test.cpp
+++ b/zero_test.cpp
@@ -763,6 +763,7 @@ void test_celebrated_equation()
  b +2.09455148154232650981 fb -8.88178419700125232339e-16
  c +2.09455148154232739799 fc +9.76996261670137755573e-15
  return value: +2.09455148154232650981 = b
+ function evaluations: +12 +12 nominal, actual
  return value: +2.09455148154232650981 (rounded)
 )--cut-here--";
 
@@ -826,15 +827,15 @@ void test_various_functions()
 //  test_a_function        (f01, root_01, -1.0, 4.0, 0.5 * 1.0e-19, __LINE__);
 //  test_a_decimal_function(f01, root_01, -1.0, 4.0, 18, __LINE__, 168);
 //  test_a_function        (f01, root_01, -1.0, 4.0, 0.5 * 1.0e-18, __LINE__);
-    test_a_decimal_function(f01, root_01, -1.0, 4.0, 17, __LINE__, 158);
+    test_a_decimal_function(f01, root_01, -1.0, 4.0, 17, __LINE__, 149);
     test_a_function        (f01, root_01, -1.0, 4.0, 0.5 * 1.0e-17, __LINE__);
-    test_a_decimal_function(f01, root_01, -1.0, 4.0, 16, __LINE__, 152);
+    test_a_decimal_function(f01, root_01, -1.0, 4.0, 16, __LINE__, 140);
     test_a_function        (f01, root_01, -1.0, 4.0, 0.5 * 1.0e-16, __LINE__);
-    test_a_decimal_function(f01, root_01, -1.0, 4.0, 15, __LINE__, 138);
+    test_a_decimal_function(f01, root_01, -1.0, 4.0, 15, __LINE__, 127);
     test_a_function        (f01, root_01, -1.0, 4.0, 0.5 * 1.0e-15, __LINE__);
-    test_a_decimal_function(f01, root_01, -1.0, 4.0, 14, __LINE__, 135);
+    test_a_decimal_function(f01, root_01, -1.0, 4.0, 14, __LINE__, 125);
     test_a_function        (f01, root_01, -1.0, 4.0, 0.5 * 1.0e-14, __LINE__);
-    test_a_decimal_function(f01, root_01, -1.0, 4.0, 12, __LINE__, 107);
+    test_a_decimal_function(f01, root_01, -1.0, 4.0, 12, __LINE__,  93);
     test_a_function        (f01, root_01, -1.0, 4.0, 0.5 * 1.0e-12, __LINE__);
 
     auto f02 = [](double x) {return std::pow(x - 1.7, 17.0);};
@@ -931,7 +932,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(163,<=,r.n_eval); // weak
+    LMI_TEST_RELATION(159,<=,r.n_eval); // weak
     LMI_TEST_RELATION(r.n_eval,<=,166); // weak
 
     d = brent_zero(eq_2_1, -100.0, 100.0, 0.5);
@@ -945,7 +946,7 @@ void test_hodgepodge()
     LMI_TEST(10 == max_n_eval_bolzano(-100.0, 100.0, 0.5, -100.0));
     LMI_TEST(98 == max_n_eval_brent  (-100.0, 100.0, 0.5, -100.0));
     LMI_TEST(r.n_eval <= 98);
-    LMI_TEST_EQUAL(18, r.n_eval); // weak
+    LMI_TEST_EQUAL(11, r.n_eval); // weak
     // Number of evaluations required:
     //   23 for brent_zero() [above]
     //   20 for decimal_root()