[Date Prev][Date Next][Thread Prev][Thread Next][Date Index][Thread Index]
[lmi-commits] [lmi] odd/expm1_log1p 42eb5d69 2/4: Use expeditious reimpl
From: |
Greg Chicares |
Subject: |
[lmi-commits] [lmi] odd/expm1_log1p 42eb5d69 2/4: Use expeditious reimplementations of transcendentals |
Date: |
Thu, 19 May 2022 06:37:38 -0400 (EDT) |
branch: odd/expm1_log1p
commit 42eb5d6972450c884064ca77836be3925c38e376
Author: Gregory W. Chicares <gchicares@sbcglobal.net>
Commit: Gregory W. Chicares <gchicares@sbcglobal.net>
Use expeditious reimplementations of transcendentals
With this change, all '.test' files created by 'make system_test' match.
Therefore, the observed differences between platforms seem to be
attributable solely to the implementations of expm1() and log1p().
---
math_functions.hpp | 30 ++++++++++++++++--------------
1 file changed, 16 insertions(+), 14 deletions(-)
diff --git a/math_functions.hpp b/math_functions.hpp
index f6718ede..c253bc27 100644
--- a/math_functions.hpp
+++ b/math_functions.hpp
@@ -24,6 +24,8 @@
#include "config.hpp"
+#include "bourn_cast.hpp"
+
#include <algorithm> // max(), min(), transform()
#include <cmath> // expm1(), log1p(), signbit()
#include <limits>
@@ -106,19 +108,19 @@ namespace lmi
#define LOGE2L 6.9314718055994530941723E-1L
#define LOG2EL 1.4426950408889634073599E0L
-long double expm1(long double x)
+inline double expm1(long double x)
{
if(std::fabs(x) < LOGE2L)
{
x *= LOG2EL;
__asm__("f2xm1" : "=t" (x) : "0" (x));
- return x;
+ return bourn_cast<double>(x);
}
else
- return std::exp(x) - 1.0;
+ return bourn_cast<double>(std::exp(x) - 1.0);
}
-long double log1p(long double x)
+inline double log1p(long double x)
{
if(std::fabs(x) < 1 - LOGE2L / 2)
{
@@ -126,10 +128,10 @@ long double log1p(long double x)
// __asm__("fldln2");
// __asm__("fxch");
__asm__("fyl2xp1" : "=t" (x) : "0" (x), "u" (y) : "st(1)");
- return x;
+ return bourn_cast<double>(x);
}
else
- return std::log(1.0 + x);
+ return bourn_cast<double>(std::log(1.0 + x));
}
} // namespace lmi
@@ -173,7 +175,7 @@ struct i_upper_n_over_n_from_i
// naively: (1+i)^(1/n) - 1
// substitute: (1+i)^n - 1 <-> std::expm1(std::log1p(i) * n)
- return std::expm1(std::log1p(i) / n);
+ return lmi::expm1(lmi::log1p(i) / n);
}
};
@@ -198,7 +200,7 @@ struct i_from_i_upper_n_over_n
{
// naively: (1+i)^n - 1
// substitute: (1+i)^n - 1 <-> std::expm1(std::log1p(i) * n)
- return std::expm1(std::log1p(i) * n);
+ return lmi::expm1(lmi::log1p(i) * n);
}
};
@@ -231,7 +233,7 @@ struct d_upper_n_from_i
// naively: n * (1 - (1+i)^(-1/n))
// substitute: (1+i)^n - 1 <-> std::expm1(std::log1p(i) * n)
- return -n * std::expm1(std::log1p(i) / -n);
+ return -n * lmi::expm1(lmi::log1p(i) / -n);
}
};
@@ -265,11 +267,11 @@ struct net_i_from_gross
// - fee *(1/n)
// )^n - 1
// substitute: (1+i)^n - 1 <-> std::expm1(std::log1p(i) * n)
- return std::expm1
+ return lmi::expm1
(
- n * std::log1p
- ( std::expm1(std::log1p(i) / n)
- - std::expm1(std::log1p(spread) / n)
+ n * lmi::log1p
+ ( lmi::expm1(lmi::log1p(i) / n)
+ - lmi::expm1(lmi::log1p(spread) / n)
- fee / n
)
);
@@ -330,7 +332,7 @@ struct coi_rate_from_q
{
// naively: 1 - (1-q)^(1/12)
// substitute: (1+i)^n - 1 <-> std::expm1(std::log1p(i) * n)
- T monthly_q = -std::expm1(std::log1p(-q) / 12);
+ T monthly_q = -lmi::expm1(lmi::log1p(-q) / 12);
if(T(1) == monthly_q)
{
throw std::logic_error("Monthly q equals unity.");