[lmi-commits] [lmi] odd/eraseme_error ee0dc1c 2/9: Fix the problem in a slapdash manner

[Greg Chicares]

branch: odd/eraseme_error
commit ee0dc1c7f63b4da98bc32924231ad023525550e9
Author: Gregory W. Chicares <gchicares@sbcglobal.net>
Commit: Gregory W. Chicares <gchicares@sbcglobal.net>

    Fix the problem in a slapdash manner
    
    At an interest rate of i, standardly defining v=1/(1+i),
    calculate future value of a stream of payments P1...Pn as
      Pn * (1+i)^1 + Pn-1 * (1+i)^2 ... + P1 * (1+i)^n
    instead of
      (Pn * v^n + Pn-1 * v^(n-1) ... + P1 * v^1) / v^(n+1)
    pointedly avoiding the NaNs that result from using the former
    formula at i = -100% --> v = 1/(1-1).
    
    The unit test demonstrates a worthwhile speed improvement:
    
      Speed test: vector of irrs, length 100, 5 decimals
      [before this commit]
    i686    4.902e-03 s mean;       4521 us least of 100 runs
    x86_64  3.128e-03 s mean;       3116 us least of 100 runs
      [after this commit]
    i686    1.488e-03 s mean;       1113 us least of 100 runs
    x86_64  1.491e-03 s mean;       1476 us least of 100 runs
---
 financial.hpp      | 50 +++++++++++++++++++++++++++++++++++++++++++++-----
 financial_test.cpp | 44 ++++++++++++++++++++++++++++++++++++--------
 2 files changed, 81 insertions(+), 13 deletions(-)

diff --git a/financial.hpp b/financial.hpp
index 70b8938..5d5bc89 100644
--- a/financial.hpp
+++ b/financial.hpp
@@ -24,8 +24,10 @@
 
 #include "config.hpp"
 
+#include "alert.hpp"
 #include "assert_lmi.hpp"
 #include "bourn_cast.hpp"
+#include "materially_equal.hpp"
 #include "mc_enum_type_enums.hpp"       // mcenum_mode
 #include "miscellany.hpp"               // ios_out_app_binary()
 #include "zero.hpp"
@@ -72,6 +74,39 @@ long double fv
 }
 
 template<typename InputIterator>
+long double future_value
+    (InputIterator first
+    ,InputIterator last
+    ,long double i
+    )
+{
+    if(first == last)
+        {
+        return 0.0L;
+        }
+    // v = 1/(1+i) is standard; u = 1+i should be, too
+    // For i=-100%, this deliberately returns zero!
+    long double const u = 1.0L + i;
+    long double un = 1.0L;
+    long double z = 0.0L;
+    for(InputIterator j = last; j != first;)
+        {
+        --j;
+        un *= u;
+        z += *j * un;
+        }
+    warning().precision(21);
+#if 0
+    auto y = fv(first, last, i);
+    if(1.0e-10 < std::fabs(z - y) || !materially_equal(z, y))
+        warning() << "UNEQUAL:\n  " << z << "\n  " << y << LMI_FLUSH;
+//  LMI_ASSERT(std::fabs(z - y) < 1.0e-13);
+    LMI_ASSERT(materially_equal(z, y));
+#endif // 0
+    return z;
+}
+
+template<typename InputIterator>
 long double npv
     (InputIterator first
     ,InputIterator last
@@ -112,19 +147,19 @@ class irr_helper
 
     long double operator()(long double i)
         {
-        return fv(first_, last_, i) - x_;
+        return future_value(first_, last_, i) - x_;
         }
 
     long double operator()()
         {
     std::ofstream ofs_trace;
-    ofs_trace.open("trace_irr.txt", ios_out_app_binary());
+////ofs_trace.open("trace_irr.txt", ios_out_app_binary());
     std::ostream os_trace(ofs_trace.rdbuf());
         root_type const z = decimal_root
             (*this
             ,-1.0       // A priori lower bound.
             ,1000.0     // Assumed upper bound.
-            ,bias_lower // Return the final bound with the lower fv.
+            ,bias_lower // Return the final bound with the lower FV.
             ,decimals_
             ,os_trace
             );
@@ -134,7 +169,12 @@ 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.0L;}
+                {
+                os_trace << "ROOT NOT BRACKETED\n" << std::endl;
+//              return z.root;
+                return 1234.5678;
+                }
             case improper_bounds:
                 {throw "IRR: improper bounds.";}
             }
@@ -166,7 +206,7 @@ long double irr
     ,int decimals
     )
 {
-    return irr_helper<InputIterator>(first, last, 0.0L, decimals)();
+    return irr_helper<InputIterator>(first, last - 1, -*(last - 1), 
decimals)();
 }
 
 template
diff --git a/financial_test.cpp b/financial_test.cpp
index 0a5630d..e4a1a63 100644
--- a/financial_test.cpp
+++ b/financial_test.cpp
@@ -39,28 +39,56 @@ int test_main(int, char*[])
     double pmts[3] = {100.0,  200.0,  300.0};
     double bfts[3] = {300.0, 1500.0, 5400.0};
 
+    std::cout << fv(pmts + 0, pmts + 1, 0.04) << std::endl;
+    std::cout << future_value(pmts + 0, pmts + 1, 0.04) << std::endl;
+
+    std::cout << fv(pmts + 0, pmts + 2, 0.04) << std::endl;
+    std::cout << future_value(pmts + 0, pmts + 2, 0.04) << std::endl;
+
+    std::cout << fv(pmts + 0, pmts + 3, 0.04) << std::endl;
+    std::cout << future_value(pmts + 0, pmts + 3, 0.04) << std::endl;
+
     // The next few tests compare floating-point quantities for exact
     // equality. Often that's inappropriate; however, the quantities
     // are integer-valued and the algorithm is designed to round them
     // exactly.
 
     irr_helper<double*> xxx(pmts, pmts + 1, bfts[0], 5);
-    LMI_TEST(2.0 == xxx());
+    LMI_TEST_EQUAL(2.0, xxx());
 
-    LMI_TEST(( 2.0 == irr_helper<double*>(pmts, pmts + 1, bfts[0], 5)()));
+    LMI_TEST_EQUAL( 2.0, irr_helper<double*>(pmts, pmts + 1, bfts[0], 5)());
 
-    LMI_TEST(( 2.0 == irr_helper<double*>(pmts, pmts + 3, bfts[2], 5)()));
+    LMI_TEST_EQUAL( 2.0, irr_helper<double*>(pmts, pmts + 3, bfts[2], 5)());
 
-    LMI_TEST((-1.0 == irr_helper<double*>(pmts, pmts + 3, 0.0    , 5)()));
+    LMI_TEST_EQUAL(-1.0, irr_helper<double*>(pmts, pmts + 3, 0.0    , 5)());
 
     // Test with arrays.
     double cash_flows[4] = {pmts[0], pmts[1], pmts[2], -bfts[2]};
-    LMI_TEST(2.0 == irr(cash_flows, 4 + cash_flows, 5));
+    LMI_TEST_EQUAL(2.0, irr(cash_flows, 4 + cash_flows, 5));
+    std::cout
+        << cash_flows[0] << "\n"
+        << cash_flows[1] << "\n"
+        << cash_flows[2] << "\n"
+        << cash_flows[3] << "\n"
+        << std::endl
+        ;
 
+    std::cout << fv(cash_flows + 0, cash_flows + 4, 0.04) << std::endl;
+    std::cout << future_value(cash_flows + 0, cash_flows + 4,  0.04) << 
std::endl;
+    std::cout << "fv of 'cash_flows' at i = -1, 0, 1:" << std::endl;
+    std::cout << future_value(cash_flows + 0, cash_flows + 4, -1.0 ) << 
std::endl;
+    std::cout << future_value(cash_flows + 0, cash_flows + 4,  0.0 ) << 
std::endl;
+    std::cout << future_value(cash_flows + 0, cash_flows + 4,  1.0 ) << 
std::endl;
+//return 0;
     // Test with vectors.
     std::vector<double> v(cash_flows, 4 + cash_flows);
-    LMI_TEST(2.0 == irr(v.begin(), v.end(), 0.0, 5));
-    LMI_TEST(2.0 == irr(v.begin(), v.end(), 5));
+    // The next test is begging to fail. It was designed so that 200%
+    // is a root; but appending a zero at the end of the stream makes
+    // -100% a root, and arguably a preferable one. Reason: any stream
+    // of payments accumulates to a future value of zero at a rate of
+    // -100%.
+    LMI_TEST_EQUAL(-1.0, irr(v.begin(), v.end(), 0.0, 5));
+    LMI_TEST_EQUAL(2.0, irr(v.begin(), v.end(), 5));
 
     std::vector<double> p; // Payments.
     std::vector<double> b; // Benefits.
@@ -184,7 +212,7 @@ int test_main(int, char*[])
     // but, given that we must return only one, wouldn't some
     // value like 0% or -100% be more suitable?
     irr(p0, b0, r0, p0.size(), p0.size(), decimals);
-    LMI_TEST_EQUAL(r0[3], 1000);
+    LMI_TEST_EQUAL(r0[3], -1);
 
     // Test fv().