[lmi-commits] [lmi] master 77cc4d65 11/11: Prefer nonstd::power() to std::pow() in a particular case

[Greg Chicares]

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

    Prefer nonstd::power() to std::pow() in a particular case
    
    The exponent in this case is always an integer in {1, 3, 6, 12}.
    For x86_64, pc-linux-gnu and MinGW-w64 yielded different results
    when std::pow() was used; nonstd::power() makes results identical.
    
    The number of multiplications is low, so there's little opportunity
    for error to accumulate. And, because the exponent is lower than 15,
    right-to-left binary exponentiation uses the optimal addition chain,
    so the number of multiplications is the lowest possible.
    
    Although 'i' may be small (so that 'v' is close to unity), 'p' is
    likely to be large, and therefore 'vp' is generally outside the range
    where using
      vp^N <=> 1 + expm1(log1p(vp - 1) * N)
    would be significantly more accurate.
---
 commutation_functions.cpp | 12 ++++++------
 1 file changed, 6 insertions(+), 6 deletions(-)

diff --git a/commutation_functions.cpp b/commutation_functions.cpp
index a456833a..4614280f 100644
--- a/commutation_functions.cpp
+++ b/commutation_functions.cpp
@@ -26,9 +26,9 @@
 #include "assert_lmi.hpp"
 #include "et_vector.hpp"                // [VECTORIZE]
 #include "ssize_lmi.hpp"
+#include "stl_extensions.hpp"           // nonstd::power()
 
 #include <algorithm>                    // rotate_copy() [VECTORIZE]
-#include <cmath>                        // pow()
 #include <functional>                   // multiplies    [VECTORIZE]
 #include <numeric>                      // partial_sum()
 
@@ -161,14 +161,14 @@ ULCommFns::ULCommFns
         // Present value of $1 one month hence.
         double vp = v * p;
         // Present value of $1 twelve months hence.
-        double vp12 = std::pow(vp, 12);
-        double vpn  = std::pow(vp, periods_per_year);
+        double vp12 = nonstd::power(vp, 12);
+        double vpn  = nonstd::power(vp, periods_per_year);
         // Twelve times a'' upper 12 (Eckley equations 28 and 31),
         // determined analytically using the geometric series theorem.
 //      double aa = 1.0;
 //      // Eckley equation (31).
-//      double sa = (1.0 - vp12) / (1.0 - std::pow(vp, 6));
-//      double qa = (1.0 - vp12) / (1.0 - std::pow(vp, 3));
+//      double sa = (1.0 - vp12) / (1.0 - nonstd::power(vp, 6));
+//      double qa = (1.0 - vp12) / (1.0 - nonstd::power(vp, 3));
 //      // Eckley equation (28).
 //      double ma = (1.0 - vp12) / (1.0 - vp);
         // The prefix k indicates the processing mode, which is
@@ -176,7 +176,7 @@ ULCommFns::ULCommFns
         double ka = 1.0;
         if(1.0 != vp)
             {
-            ka = (1.0 - vp12) / (1.0 - std::pow(vp, months_per_period));
+            ka = (1.0 - vp12) / (1.0 - nonstd::power(vp, months_per_period));
             }
         kd[j] = ka * ad[j];
         kc[j] = ka * ad[j] * v * q;