[lmi-commits] [lmi] master e59df26 14/23: Refactor

[Greg Chicares]

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

    Refactor
    
    Instead of solving for a rounded value of an unrounded objective
    function (which in practice actually does perform rounding anyway),
    solve for an unrounded value of a rounded objective function.
    
    This is crucial for binary64_midpoint(), which, given input bounds
    such as [0, 1e300], might iterate to 1e-100, eliminating half the
    possible floating-point values. Rounding that important 1e-100 iterate
    back to zero defeats that purpose.
    
    Practical necessity aside, this way is also much cleaner. There's no
    need to intrude rounding into the root-finding algorithm.
    
    Incidentally added code to show unrounded and rounded return values
    in the optional trace.
---
 zero.hpp      | 49 +++++++++++++++++--------------------------------
 zero_test.cpp | 45 ++++++++++++++++++++++++---------------------
 2 files changed, 41 insertions(+), 53 deletions(-)

diff --git a/zero.hpp b/zero.hpp
index 9222a46..133154a 100644
--- a/zero.hpp
+++ b/zero.hpp
@@ -32,7 +32,6 @@
 #include <cmath>                        // fabs(), isfinite(), isnan(), pow()
 #include <cstdint>                      // uint64_t
 #include <cstring>                      // memcpy()
-#include <functional>                   // function(), identity()
 #include <iomanip>                      // setprecision(), setw()
 #include <limits>
 #include <numeric>                      // midpoint()
@@ -192,24 +191,6 @@ inline double binary64_midpoint(double d0, double d1)
     return z;
 }
 
-namespace detail
-{
-using RoundT = std::function<double(double)>;
-} // namespace detail
-
-/// Workaround for clang--see:
-///    https://lists.nongnu.org/archive/html/lmi/2021-07/msg00001.html
-/// It is hoped that this can be replaced by std::identity soon.
-
-struct lmi_identity
-{
-    using is_transparent = void;
-
-    template<typename T>
-    constexpr T&& operator()(T&& t) const noexcept
-        {return std::forward<T>(t);}
-};
-
 /// Return a zero z of a function f within input bounds [a,b].
 ///
 /// Preconditions: bounds are distinct after rounding; and either
@@ -386,8 +367,6 @@ root_type lmi_root
     ,double          tolerance
     ,std::ostream&   os_trace  = null_stream()
     ,root_bias       bias      = bias_none
-//  ,detail::RoundT  round_dec = std::identity()
-    ,detail::RoundT  round_dec = lmi_identity()
     )
 {
     constexpr double epsilon {std::numeric_limits<double>::epsilon()};
@@ -405,9 +384,9 @@ root_type lmi_root
         ;
 
     // Declarations must precede lambda.
-    double a  {};
+    double a  {bound0};
     double fa {};
-    double b  {};
+    double b  {bound1};
     double fb {};
     double c  {};
     double fc {};
@@ -442,12 +421,10 @@ root_type lmi_root
 
     double t = tolerance;
 
-    a = round_dec(bound0);
-    b = round_dec(bound1);
-
     if(a == b)
         {
         recapitulate();
+        os_trace << " return value: " << a << " = a" << std::endl;
         return {a, improper_bounds, n_iter, n_eval};
         }
 
@@ -456,6 +433,7 @@ root_type lmi_root
     if(0.0 == fa) // Note 0.
         {
         recapitulate();
+        os_trace << " return value: " << a << " = a" << std::endl;
         return {a, root_is_valid, n_iter, n_eval};
         }
 
@@ -465,6 +443,7 @@ root_type lmi_root
     if(0.0 == fb) // Note 0 [bis].
         {
         recapitulate();
+        os_trace << " return value: " << b << " = b" << std::endl;
         return {b, root_is_valid, n_iter, n_eval};
         }
 
@@ -473,6 +452,7 @@ root_type lmi_root
     if(std::isnan(fa) || std::isnan(fb) || signum(fa) == signum(fb))
         {
         recapitulate();
+        os_trace << " return value: " << 0.0 << " = zero" << std::endl;
         return {0.0, root_not_bracketed, n_iter, n_eval};
         }
 
@@ -512,11 +492,13 @@ root_type lmi_root
                 )
                 {
                 recapitulate();
+                os_trace << " return value: " << b << " = b" << std::endl;
                 return {b, root_is_valid, n_iter, n_eval};
                 }
             else if(std::fabs(m) <= 2.0 * epsilon * std::fabs(c) + t)
                 {
                 recapitulate();
+                os_trace << " return value: " << c << " = c" << std::endl;
                 return {c, root_is_valid, n_iter, n_eval};
                 }
             else
@@ -610,7 +592,6 @@ root_type lmi_root
             {
             b -= tol;
             }
-        b = round_dec(b);
 
         if(b == a) // Note 4.
             {
@@ -655,16 +636,19 @@ root_type decimal_root
     )
 {
     round_to<double> const round_dec {decimals, r_to_nearest};
+    auto fr = [&](double x) {return f(round_dec(x));};
 
-    return lmi_root
-        (f
-        ,bound0
-        ,bound1
+    auto z = lmi_root
+        (fr
+        ,round_dec(bound0)
+        ,round_dec(bound1)
         ,0.5 * std::pow(10.0, -decimals)
         ,os_trace
         ,bias
-        ,detail::RoundT(round_dec)
         );
+    z.root = round_dec(z.root);
+    os_trace << " return value: " << z.root << " (rounded)" << std::endl;
+    return z;
 }
 
 /// An instrumented translation of Brent's reference implementation.
@@ -852,6 +836,7 @@ double brent_zero
             {goto extrapolate;}
         }
     recapitulate();
+    os_trace << " return value: " << b << " = b" << std::endl;
     return b;
 }
 
diff --git a/zero_test.cpp b/zero_test.cpp
index 6e4f47e..2b8025c 100644
--- a/zero_test.cpp
+++ b/zero_test.cpp
@@ -705,13 +705,13 @@ void test_biases()
 ///    i686-w64-mingw32     x86_64-pc-linux-gnu
 ///   --------lmi-------     --------lmi-------     -----mathworks----
 ///   2.5600000000000001     2.5600000000000001     2.5600000000000001
-///   1.0980323260716796  +ϵ 1.0980323260716793     1.0980323260716793
-///   1.783216881610604   +ϵ 1.783216881610604   +ϵ 1.7832168816106038
-///   2.2478393639958036     2.2478393639958036     2.2478393639958036
+///   1.0980323260716793     1.0980323260716793     1.0980323260716793
+///   1.7832168816106038     1.7832168816106038     1.7832168816106038
+///   2.2478393639958032 -2ϵ 2.2478393639958032 -2ϵ 2.2478393639958036
 ///   2.0660057758331045     2.0660057758331045     2.0660057758331045
 ///   2.0922079131171945     2.0922079131171945     2.0922079131171945
-///   2.0945566700001774 -2ϵ 2.0945566700001779     2.0945566700001779
-///   2.0945514746903111     2.0945514746903111     2.0945514746903111
+///   2.0945566700001779     2.0945566700001779     2.0945566700001779
+///   2.0945514746903116 +2ϵ 2.0945514746903111     2.0945514746903111
 ///   2.0945514815423065     2.0945514815423065     2.0945514815423065
 ///   2.0945514815423265     2.0945514815423265     2.0945514815423265
 ///   2.0945514815423274     2.0945514815423274     2.0945514815423274
@@ -744,13 +744,13 @@ void test_celebrated_equation()
   0   2 j -2.5600000000000001 -16.657216000000002 2.5600000000000001 
6.6572160000000018 -2.5600000000000001 -16.657216000000002
   0   3 L 2.5600000000000001 6.6572160000000018 1.0980323260716793 
-5.8721945393772152 -2.5600000000000001 -16.657216000000002
   1   3 j 2.5600000000000001 6.6572160000000018 1.0980323260716793 
-5.8721945393772152 2.5600000000000001 6.6572160000000018
-  1   4 L 1.0980323260716793 -5.8721945393772152 1.783216881610604 
-2.8960493667789873 2.5600000000000001 6.6572160000000018
-  2   5 Q 1.783216881610604 -2.8960493667789873 2.2478393639958036 
1.8621631139566732 2.5600000000000001 6.6572160000000018
-  3   5 j 1.783216881610604 -2.8960493667789873 2.2478393639958036 
1.8621631139566732 1.783216881610604 -2.8960493667789873
-  3   6 L 2.2478393639958036 1.8621631139566732 2.0660057758331045 
-0.3135140955237814 1.783216881610604 -2.8960493667789873
-  4   6 j 2.2478393639958036 1.8621631139566732 2.0660057758331045 
-0.3135140955237814 2.2478393639958036 1.8621631139566732
-  4   7 L 2.0660057758331045 -0.3135140955237814 2.0922079131171945 
-0.026123094109737011 2.2478393639958036 1.8621631139566732
-  5   8 Q 2.0922079131171945 -0.026123094109737011 2.0945566700001779 
5.7910818359374616e-05 2.2478393639958036 1.8621631139566732
+  1   4 L 1.0980323260716793 -5.8721945393772152 1.7832168816106038 
-2.8960493667789873 2.5600000000000001 6.6572160000000018
+  2   5 Q 1.7832168816106038 -2.8960493667789873 2.2478393639958032 
1.862163113956667 2.5600000000000001 6.6572160000000018
+  3   5 j 1.7832168816106038 -2.8960493667789873 2.2478393639958032 
1.862163113956667 1.7832168816106038 -2.8960493667789873
+  3   6 L 2.2478393639958032 1.862163113956667 2.0660057758331045 
-0.3135140955237814 1.7832168816106038 -2.8960493667789873
+  4   6 j 2.2478393639958032 1.862163113956667 2.0660057758331045 
-0.3135140955237814 2.2478393639958032 1.862163113956667
+  4   7 L 2.0660057758331045 -0.3135140955237814 2.0922079131171945 
-0.026123094109737011 2.2478393639958032 1.862163113956667
+  5   8 Q 2.0922079131171945 -0.026123094109737011 2.0945566700001779 
5.7910818359374616e-05 2.2478393639958032 1.862163113956667
   6   8 j 2.0922079131171945 -0.026123094109737011 2.0945566700001779 
5.7910818359374616e-05 2.0922079131171945 -0.026123094109737011
   6   9 L 2.0945566700001779 5.7910818359374616e-05 2.0945514746903111 
-7.6478343657981895e-08 2.0922079131171945 -0.026123094109737011
   7   9 j 2.0945566700001779 5.7910818359374616e-05 2.0945514746903111 
-7.6478343657981895e-08 2.0945566700001779 5.7910818359374616e-05
@@ -762,6 +762,8 @@ void test_celebrated_equation()
 10 iterations, 12 evaluations; final interval:
  b +2.09455148154232650981 fb -8.88178419700125232339e-16
  c +2.09455148154232739799 fc +9.76996261670137755573e-15
+ return value: +2.09455148154232650981 = b
+ return value: +2.09455148154232650981 (rounded)
 )--cut-here--";
 
     LMI_TEST_EQUAL(verified, oss.str());
@@ -824,20 +826,20 @@ 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__, 167);
+    test_a_decimal_function(f01, root_01, -1.0, 4.0, 17, __LINE__, 158);
     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__, 155);
+    test_a_decimal_function(f01, root_01, -1.0, 4.0, 16, __LINE__, 152);
     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__, 143);
+    test_a_decimal_function(f01, root_01, -1.0, 4.0, 15, __LINE__, 138);
     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__, 128);
+    test_a_decimal_function(f01, root_01, -1.0, 4.0, 14, __LINE__, 135);
     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__, 117);
+    test_a_decimal_function(f01, root_01, -1.0, 4.0, 12, __LINE__, 107);
     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);};
     auto root_02 = 1.7;
-    test_a_decimal_function(f02, root_02, 0.0, 2.0, 17     , __LINE__, 145);
+    test_a_decimal_function(f02, root_02, 0.0, 2.0, 17     , __LINE__, 148);
     test_a_function        (f02, root_02, 0.0, 2.0, 1.0e-15, __LINE__);
 
     auto f03 = [](double x) {return std::cos(x) - 0.999;};
@@ -929,7 +931,8 @@ 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(168 <= r.n_eval && r.n_eval <= 172); // weak
+    LMI_TEST_RELATION(163,<=,r.n_eval); // weak
+    LMI_TEST_RELATION(r.n_eval,<=,166); // weak
 
     d = brent_zero(eq_2_1, -100.0, 100.0, 0.5);
     zeta = -100.0;
@@ -942,7 +945,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(20, r.n_eval); // weak
+    LMI_TEST_EQUAL(18, r.n_eval); // weak
     // Number of evaluations required:
     //   23 for brent_zero() [above]
     //   20 for decimal_root()
@@ -968,7 +971,7 @@ void test_hodgepodge()
     LMI_TEST(  53 == max_n_eval_bolzano(-100.0, 100.0, 0.0, -100.0));
     LMI_TEST(2807 == max_n_eval_brent  (-100.0, 100.0, 0.0, -100.0));
     LMI_TEST(r.n_eval <= 2807);
-    LMI_TEST_EQUAL(66, r.n_eval); // weak
+    LMI_TEST_EQUAL(67, r.n_eval); // weak
 
     r = decimal_root(signum_offset, -1.0, 1.0, bias_none, 13);
     LMI_TEST(root_is_valid == r.validity);