[lmi-commits] [lmi] master bcd3e721 2/8: Make relative-error function generally available

[Greg Chicares]

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

    Make relative-error function generally available
    
    Its name will become less terse in the next commit; this commit is
    intended to make the transplantation easy to see.
    
    Inaugurated the convention of writing "/*constexpr*/ where it can be
    replaced by "constexpr" when most of <cmath> becomes constexpr, in
    C++23 one hopes.
---
 math_functions.hpp      | 35 ++++++++++++++++++++++++++++++++++-
 math_functions_test.cpp | 44 ++++++++++++++++++++------------------------
 2 files changed, 54 insertions(+), 25 deletions(-)

diff --git a/math_functions.hpp b/math_functions.hpp
index fa09c375..6fbe59d5 100644
--- a/math_functions.hpp
+++ b/math_functions.hpp
@@ -25,7 +25,7 @@
 #include "config.hpp"
 
 #include <algorithm>                    // max(), min(), transform()
-#include <cmath>                        // signbit()
+#include <cmath>                        // fabs(), signbit()
 #include <limits>
 #include <numeric>                      // midpoint(), partial_sum()
 #include <stdexcept>
@@ -96,6 +96,39 @@ inline T outward_quotient(T numerator, T denominator)
     return (0 < numerator == 0 < denominator) ? x + y : x - y;
 }
 
+/// Absolute value of relative error.
+///
+/// Sometimes one of the arguments is a canonical expected value, and
+///   (observed-expected)/expected
+/// is wanted, with due regard to sign; in that case, code (o-e)/e
+/// directly. This function template, OTOH, is designed for the more
+/// general case where neither argument is favored; then a denominator
+/// must be chosen somehow, and the algebraic sign of the return value
+/// would generally not be useful.
+///
+/// Here the smaller of the absolute values of the arguments is chosen
+/// as the denominator. Alternatively, the greater, or some average,
+/// could have been chosen; but choosing the smaller gives a more
+/// conservative (i.e., larger) value.
+///
+/// Only floating point arguments are allowed, because no compelling
+/// use case for integer types is foreseen. Similarly, both arguments
+/// are constrained to be of the same type, because allowing different
+/// types seems unimportant.
+
+template<typename T>
+/*constexpr*/ T rel_err(T t, T u)
+{
+    static_assert(std::is_floating_point_v<T>);
+    constexpr T inf {std::numeric_limits<T>::infinity()};
+    auto const denominator {std::min(std::fabs(t), std::fabs(u))};
+    return
+          (0.0 == t && 0.0 == u) ? 0.0
+        : (0.0 == denominator)   ? inf
+        :                          std::fabs(t - u) / denominator
+        ;
+}
+
 /// Signed zeros, for comparison tests.
 
 enum signed_zero
diff --git a/math_functions_test.cpp b/math_functions_test.cpp
index d67ab439..a3d26b81 100644
--- a/math_functions_test.cpp
+++ b/math_functions_test.cpp
@@ -477,6 +477,24 @@ void test_compound_interest()
         );
 }
 
+void test_relative_error()
+{
+    constexpr double inf {std::numeric_limits<double>::infinity()};
+    constexpr double big {std::numeric_limits<double>::max()};
+
+    LMI_TEST_EQUAL(inf, rel_err(0.0, -2.0));
+    LMI_TEST_EQUAL(inf, rel_err(0.0, -1.0));
+    LMI_TEST_EQUAL(inf, rel_err(0.0, -0.5));
+    LMI_TEST_EQUAL(0.0, rel_err(0.0,  0.0));
+    LMI_TEST_EQUAL(inf, rel_err(0.0,  0.5));
+    LMI_TEST_EQUAL(inf, rel_err(0.0,  1.0));
+    LMI_TEST_EQUAL(inf, rel_err(0.0,  2.0));
+    LMI_TEST_EQUAL(0.0, rel_err(1.0,  1.0));
+    LMI_TEST_EQUAL(2.0, rel_err(1.0, -1.0));
+    LMI_TEST_EQUAL(big, rel_err(1.0,  big));
+    LMI_TEST_EQUAL(inf, rel_err(big, -big));
+}
+
 void test_signed_zero()
 {
     constexpr double inf  {std::numeric_limits<double>::infinity ()};
@@ -665,30 +683,6 @@ void test_expm1_log1p()
     double const i3 = lmi::expm1(-1.1512925464970228);
     LMI_TEST(materially_equal(-0.68377223398240425, i3, 1.0e-11));
 
-    constexpr double inf {std::numeric_limits<double>::infinity()};
-    constexpr double big {std::numeric_limits<double>::max()};
-    // Absolute value of relative error.
-    auto rel_err = [](double t, double u)
-        {
-        auto const denominator {std::min(std::fabs(t), std::fabs(u))};
-        return
-              (0.0 == t && 0.0 == u) ? 0.0
-            : (0.0 == denominator)   ? inf
-            :                          std::fabs(t - u) / denominator
-            ;
-        };
-    LMI_TEST_EQUAL(inf, rel_err(0.0, -2.0));
-    LMI_TEST_EQUAL(inf, rel_err(0.0, -1.0));
-    LMI_TEST_EQUAL(inf, rel_err(0.0, -0.5));
-    LMI_TEST_EQUAL(0.0, rel_err(0.0,  0.0));
-    LMI_TEST_EQUAL(inf, rel_err(0.0,  0.5));
-    LMI_TEST_EQUAL(inf, rel_err(0.0,  1.0));
-    LMI_TEST_EQUAL(inf, rel_err(0.0,  2.0));
-    LMI_TEST_EQUAL(0.0, rel_err(1.0,  1.0));
-    LMI_TEST_EQUAL(2.0, rel_err(1.0, -1.0));
-    LMI_TEST_EQUAL(big, rel_err(1.0,  big));
-    LMI_TEST_EQUAL(inf, rel_err(big, -big));
-
     // Test fdlibm vs. C RTL for many parameters.
     int    err_count0 {0};
     int    err_count1 {0};
@@ -838,6 +832,8 @@ int test_main(int, char*[])
 
     test_compound_interest();
 
+    test_relative_error();
+
     test_signed_zero();
 
     test_signum<bool         >(__FILE__, __LINE__);