[lmi-commits] [lmi] odd/eraseme_long_double_irr 4ba46a03 3/3: Reimplement IRR in terms of double rather than long double

[Greg Chicares]

branch: odd/eraseme_long_double_irr
commit 4ba46a034ece569fd357d6fba48aa35bd195ffa3
Author: Gregory W. Chicares <gchicares@sbcglobal.net>
Commit: Gregory W. Chicares <gchicares@sbcglobal.net>

    Reimplement IRR in terms of double rather than long double
    
    Results of 'financial_test' unit test, before and after, comparing
    NPV of a stream of 100 elements at the calculated IRR rate (which
    would be exactly zero, in infinite precision), and also comparing
    IRR timings:
    
    before (long double)
      i686-w64-mingw32     -9.86988e-14
        iterator  form: 2.587e-02 s mean;      24182 us least of  39 runs
        container form: 2.469e-02 s mean;      24373 us least of  41 runs
      x86_64-w64-mingw32   -9.86988e-14
        iterator  form: 2.220e-02 s mean;      21482 us least of  46 runs
        container form: 2.178e-02 s mean;      21475 us least of  46 runs
      x86_64-pc-linux-gnu  -9.86988e-14
        iterator  form: 1.433e-02 s mean;      13447 us least of  70 runs
        container form: 1.373e-02 s mean;      13453 us least of  73 runs
    
    after (double)
      i686-w64-mingw32     -9.86988e-14
        iterator  form: 2.610e-02 s mean;      25739 us least of  39 runs
        container form: 2.628e-02 s mean;      25837 us least of  39 runs
      x86_64-w64-mingw32   -9.77707e-12
        iterator  form: 2.329e-02 s mean;      21414 us least of  43 runs
        container form: 2.202e-02 s mean;      21393 us least of  46 runs
      x86_64-pc-linux-gnu  -9.77707e-12
        iterator  form: 1.349e-02 s mean;      13415 us least of  75 runs
        container form: 1.346e-02 s mean;      13410 us least of  75 runs
    
    The original (long double) implementation seems preferable:
     - it's more accurate for x86_64, by two orders of magnitude; and
     - it's as fast for x86_64, and maybe faster for i686 (x87); and
     - the resulting NPV is identical across architectures.
---
 financial.hpp      | 27 ++++++++++++---------------
 financial_test.cpp | 37 +++++++++++++++++++++++--------------
 2 files changed, 35 insertions(+), 29 deletions(-)

diff --git a/financial.hpp b/financial.hpp
index 7b18939b..1969c79e 100644
--- a/financial.hpp
+++ b/financial.hpp
@@ -46,26 +46,23 @@ class calendar_date;
 // When no root is bracketed, -100% is always conservative enough;
 // but if a root is known to exceed the a priori upper bound, then
 // perhaps that upper bound could be returned instead.
-//
-// Is it really advantageous to use long double? Why not use a
-// template type argument instead?
 
 template<typename InputIterator>
-long double fv
+double fv
     (InputIterator first
     ,InputIterator last
-    ,long double i
+    ,double        i
     )
 {
     if(first == last)
         {
-        return 0.0L;
+        return 0.0;
         }
 
     // Symbol v, meaning 1/(1+i), is standard. A corresponding
     // symbol u, meaning (1+i), is not standard, but should be.
-    long double const u = 1.0L + i;
-    long double z = 0.0L;
+    double const u = 1.0 + i;
+    double z = 0.0;
     for(InputIterator j = first; j != last; ++j)
         {
         z += *j;
@@ -81,7 +78,7 @@ class irr_helper
     irr_helper
         (InputIterator first
         ,InputIterator last
-        ,long double   x
+        ,double        x
         ,int           decimals
         )
         :first_    {first}
@@ -90,12 +87,12 @@ class irr_helper
         ,decimals_ {decimals}
         {}
 
-    long double operator()(long double i)
+    double operator()(double i)
         {
         return fv(first_, last_, i) - x_;
         }
 
-    long double operator()()
+    double operator()()
         {
 // Uncomment to trace iterations:
 //      std::ofstream ofs_trace;
@@ -116,7 +113,7 @@ class irr_helper
                 {return z.root;}
             // Return -100% if NPVs of a priori bounds have same sign.
             case root_not_bracketed:
-                {return -1.0L;}
+                {return -1.0;}
             case improper_bounds:
                 {throw "IRR: improper bounds.";}
             }
@@ -126,15 +123,15 @@ class irr_helper
   private:
     InputIterator first_;
     InputIterator last_;
-    long double   x_;
+    double        x_;
     int           decimals_;
 };
 
 template<typename InputIterator>
-long double irr
+double irr
     (InputIterator first
     ,InputIterator last
-    ,long double x
+    ,double        x
     ,int decimals
     )
 {
diff --git a/financial_test.cpp b/financial_test.cpp
index 41c8440e..c71c7110 100644
--- a/financial_test.cpp
+++ b/financial_test.cpp
@@ -41,19 +41,19 @@
 /// to avoid division by zero.
 
 template<typename InputIterator>
-long double pv
+double pv
     (InputIterator first
     ,InputIterator last
-    ,long double i
+    ,double        i
     )
 {
     if(first == last)
         {
-        return 0.0L;
+        return 0.0;
         }
-    long double const v = 1.0L / (1.0L + i);
-    long double vn = 1.0L;
-    long double z = *first;
+    double const v = 1.0 / (1.0 + i);
+    double vn = 1.0;
+    double z = *first;
     InputIterator j = first;
     while(++j != last)
         {
@@ -113,9 +113,9 @@ int test_main(int, char*[])
     double pmts[3] = {100.0,  200.0,  300.0};
     double bfts[3] = {300.0, 1500.0, 5400.0};
 
-    LMI_TEST(materially_equal(104.0000L, fv(pmts + 0, pmts + 1, 0.04)));
-    LMI_TEST(materially_equal(316.1600L, fv(pmts + 0, pmts + 2, 0.04)));
-    LMI_TEST(materially_equal(640.8064L, fv(pmts + 0, pmts + 3, 0.04)));
+    LMI_TEST(materially_equal(104.0000, fv(pmts + 0, pmts + 1, 0.04)));
+    LMI_TEST(materially_equal(316.1600, fv(pmts + 0, pmts + 2, 0.04)));
+    LMI_TEST(materially_equal(640.8064, fv(pmts + 0, pmts + 3, 0.04)));
 
     // The next few tests compare floating-point quantities for exact
     // equality. Often that's inappropriate; however, the quantities
@@ -186,13 +186,22 @@ int test_main(int, char*[])
     // For any stream, NPV at the IRR rate should ideally be zero.
     std::vector<double> q{p};
     q.push_back(-b.back());
-    // This NPV is -9.777068044058979E-12 in a gnumeric spreadsheet,
-    // versus -9.86988e-014 with MinGW-w64 gcc-6.3.0
-    // versus  3.6593e-013  with i686-w64-mingw32 gcc-10.0
-    LMI_TEST(std::fabs(pv(q.begin(), q.end(), results.back())) <= 1e-12);
+    // On a 2.4GHz E5-2630 v3, on 2022-04-30, this NPV
+    std::cout << pv(q.begin(), q.end(), results.back()) << std::endl;
+    // is
+    //        -9.77707E-12 in gnumeric,
+    // versus, for a 'long double' implementation here:
+    //        -9.86988e-14 i686-w64-mingw32    gcc-10
+    //        -9.86988e-14 x86_64-w64-mingw32  gcc-10
+    //        -9.86988e-14 x86_64-pc-linux-gnu gcc-11.2.0
+    // versus, for a 'double' implementation here:
+    //        -9.86988e-14 i686-w64-mingw32    gcc-10
+    //        -9.77707e-12 x86_64-w64-mingw32  gcc-10
+    //        -9.77707e-12 x86_64-pc-linux-gnu gcc-11.2.0
+    LMI_TEST(std::fabs(pv(q.begin(), q.end(), results.back())) <= 1e-11);
 
     // Trivially, NPV at 0% interest is summation.
-    LMI_TEST(materially_equal(-4950.0L, pv(q.begin(), q.end(), 0.0)));
+    LMI_TEST(materially_equal(-4950.0, pv(q.begin(), q.end(), 0.0)));
 
     // Test const vectors.
     std::vector<double> const cp(p);