[lmi-commits] [lmi] odd/expm1_log1p 9eea8419 1/4: Expeditiously reimplement expm1() and log1p()

[Greg Chicares]

branch: odd/expm1_log1p
commit 9eea8419fa76c7759f0548eeeaebb3a04063c9d5
Author: Gregory W. Chicares <gchicares@sbcglobal.net>
Commit: Gregory W. Chicares <gchicares@sbcglobal.net>

    Expeditiously reimplement expm1() and log1p()
    
    Different results are observed on pc-linux-gnu vs msw, so implement
    these two functions in a casual way in order to determine whether
    they're the cause. Perhaps they are:
    
    x86_64, sse, gcc, posix
      1.74560101501691655734305 std::expm1(1.01)
      0.00995033085316808334209 std::log1p(0.01)
      0.00327373978219886374239 std::expm1(std::log1p(0.04) / 12)
      1.74560101501691651831177 lmi::expm1(1.01)
      0.00995033085316808305410 lmi::log1p(0.01)
      0.00327373978219886392640 lmi::expm1(lmi::log1p(0.04) / 12)
    
    x86_64, sse, gcc, msw
      1.74560101501691633529845 std::expm1(1.01)
      0.00995033085316808334209 std::log1p(0.01)
      0.00327373978219886417607 std::expm1(std::log1p(0.04) / 12)
      1.74560101501691651831177 lmi::expm1(1.01)
      0.00995033085316808305410 lmi::log1p(0.01)
      0.00327373978219886392640 lmi::expm1(lmi::log1p(0.04) / 12)
---
 math_functions.hpp      | 34 ++++++++++++++++++++++++++++++++++
 math_functions_test.cpp | 10 ++++++++++
 2 files changed, 44 insertions(+)

diff --git a/math_functions.hpp b/math_functions.hpp
index 31e8530f..f6718ede 100644
--- a/math_functions.hpp
+++ b/math_functions.hpp
@@ -99,6 +99,40 @@ T signum(T t)
     return (0 == t) ? 0 : std::signbit(t) ? -1 : 1;
 }
 
+namespace lmi
+{
+// Cf. lmi commit 16afb361c29 .
+
+#define LOGE2L  6.9314718055994530941723E-1L
+#define LOG2EL  1.4426950408889634073599E0L
+
+long double expm1(long double x)
+{
+  if(std::fabs(x) < LOGE2L)
+      {
+      x *= LOG2EL;
+      __asm__("f2xm1" : "=t" (x) : "0" (x));
+      return x;
+      }
+  else
+    return std::exp(x) - 1.0;
+}
+
+long double log1p(long double x)
+{
+  if(std::fabs(x) < 1 - LOGE2L / 2)
+      {
+      long double y = LOGE2L;
+//    __asm__("fldln2");
+//    __asm__("fxch");
+      __asm__("fyl2xp1" : "=t" (x) : "0" (x), "u" (y) : "st(1)");
+      return x;
+      }
+  else
+    return std::log(1.0 + x);
+}
+} // namespace lmi
+
 // Actuarial functions.
 //
 // Some inputs are nonsense, like interest rates less than 100%.
diff --git a/math_functions_test.cpp b/math_functions_test.cpp
index c2795f1e..11f08c3d 100644
--- a/math_functions_test.cpp
+++ b/math_functions_test.cpp
@@ -196,6 +196,16 @@ void sample_results()
         << "  double prec, std::pow\n"
         ;
 
+    std::cout
+        << "  " << std::expm1(1.01) << " std::expm1(1.01)\n"
+        << "  " << std::log1p(0.01) << " std::log1p(0.01)\n"
+        << "  " << std::expm1(std::log1p(0.04) / 12) << " 
std::expm1(std::log1p(0.04) / 12)\n"
+        << "  " << lmi::expm1(1.01) << " lmi::expm1(1.01)\n"
+        << "  " << lmi::log1p(0.01) << " lmi::log1p(0.01)\n"
+        << "  " << lmi::expm1(lmi::log1p(0.04) / 12) << " 
lmi::expm1(lmi::log1p(0.04) / 12)\n"
+        << std::endl
+        ;
+
     fenv_initialize();
 }