[Date Prev][Date Next][Thread Prev][Thread Next][Date Index][Thread Index]
[lmi-commits] [lmi] odd/expm1_log1p 9eea8419 1/4: Expeditiously reimplem
From: |
Greg Chicares |
Subject: |
[lmi-commits] [lmi] odd/expm1_log1p 9eea8419 1/4: Expeditiously reimplement expm1() and log1p() |
Date: |
Thu, 19 May 2022 06:37:38 -0400 (EDT) |
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();
}