[lmi-commits] [lmi] odd/expm1_log1p 42eb5d69 2/4: Use expeditious reimplementations of transcendentals

[Greg Chicares]

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.");