[lmi-commits] [lmi] master dd3abb11 18/22: Demonstrate a defect

[Greg Chicares]

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

    Demonstrate a defect
---
 Makefile.am             |  2 ++
 math_functions.hpp      |  9 +++++++++
 math_functions_test.cpp | 43 +++++++++++++++++++++++++++++++++++++++++++
 objects.make            |  2 ++
 4 files changed, 56 insertions(+)

diff --git a/Makefile.am b/Makefile.am
index 3f512af9..ea4982f2 100644
--- a/Makefile.am
+++ b/Makefile.am
@@ -876,6 +876,8 @@ materially_equal_test_LDADD = \
   libtest_common.la
 
 math_functions_test_LDADD = \
+  fdlibm_expm1.c \
+  fdlibm_log1p.c \
   libtest_common.la
 
 mc_enum_test_SOURCES = \
diff --git a/math_functions.hpp b/math_functions.hpp
index 19c3e5e8..b32bd3c2 100644
--- a/math_functions.hpp
+++ b/math_functions.hpp
@@ -32,6 +32,15 @@
 #include <type_traits>                  // /is_.*/
 #include <vector>
 
+extern "C" double fdlibm_expm1(double);
+extern "C" double fdlibm_log1p(double);
+
+namespace lmi
+{
+inline double expm1(double z) {return fdlibm_expm1(z);}
+inline double log1p(double z) {return fdlibm_log1p(z);}
+} // namespace lmi
+
 // TODO ?? Write functions here for other refactorable uses of
 // std::pow() throughout lmi, to facilitate reuse and unit testing.
 
diff --git a/math_functions_test.cpp b/math_functions_test.cpp
index c2795f1e..2690607d 100644
--- a/math_functions_test.cpp
+++ b/math_functions_test.cpp
@@ -147,6 +147,47 @@ struct i_upper_n_over_n_from_i_T
 };
 } // Unnamed namespace.
 
+/// Test fdlibm expm1() and log1p().
+///
+/// Testing for exact floating-point equality seems to be a patent
+/// mistake, but should work for the fdlibm implementations on
+/// any implementation where double is binary64.
+
+void test_expm1_log1p()
+{
+    std::cout.setf(std::ios_base::fixed, std::ios_base::floatfield);
+    std::cout.precision(23);
+
+    double const x = lmi::expm1(1.01);
+    double const y = lmi::log1p(0.01);
+    double const z = lmi::expm1(lmi::log1p(0.04) / 12);
+
+#if 0
+    // digits      1 23456789012345678901234
+    LMI_TEST_EQUAL(1.74560101501691655734305, x);
+    // digits          123456789012345678901
+    LMI_TEST_EQUAL(0.00995033085316808334209, y);
+    // digits          123456789012345678901
+    LMI_TEST_EQUAL(0.00327373978219886374239, z);
+#else // not 0
+    std::cout << "Those tests failed, so compare..." << std::endl;
+
+    double const a = std::expm1(1.01);
+    double const b = std::log1p(0.01);
+    double const c = std::expm1(std::log1p(0.04) / 12);
+
+    std::cout
+        << "  " << x << " lmi::expm1(1.01)\n"
+        << "  " << a << " std::expm1(1.01)\n"
+        << "  " << y << " lmi::log1p(0.01)\n"
+        << "  " << b << " std::log1p(0.01)\n"
+        << "  " << z << " lmi::expm1(lmi::log1p(0.04) / 12)\n"
+        << "  " << c << " std::expm1(std::log1p(0.04) / 12)\n"
+        << std::endl
+        ;
+#endif // not 0
+}
+
 /// This function isn't a unit test per se. Its purpose is to show
 /// how a sample calculation is affected by
 ///   exponential versus power method,
@@ -564,6 +605,8 @@ int test_main(int, char*[])
 
     assay_speed();
 
+    test_expm1_log1p();
+
     sample_results();
 
     return 0;
diff --git a/objects.make b/objects.make
index 85d38976..e102e0ac 100644
--- a/objects.make
+++ b/objects.make
@@ -791,6 +791,8 @@ materially_equal_test$(EXEEXT): \
 
 math_functions_test$(EXEEXT): \
   $(common_test_objects) \
+  fdlibm_expm1.o \
+  fdlibm_log1p.o \
   math_functions_test.o \
   timer.o \