[lmi-commits] [lmi] master 205497ae 8/9: Reimplement max_modal_premium()

[Greg Chicares]

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

    Reimplement max_modal_premium()
    
    The number of places where ldbl_eps_plus_one_times() is called has
    attained its ideal value.
---
 ieee754.hpp           | 18 -----------
 ul_utilities.cpp      | 75 +++++++++++++++++++++++++++++++++++++++++--
 ul_utilities_test.cpp | 88 +++++++++++++++++++++++++++++++++++++++++++++++++++
 3 files changed, 161 insertions(+), 20 deletions(-)

diff --git a/ieee754.hpp b/ieee754.hpp
index 59d2d371..922f6cef 100644
--- a/ieee754.hpp
+++ b/ieee754.hpp
@@ -110,22 +110,4 @@ inline bool is_infinite(T t)
     return has_inf && (pos_inf == t || neg_inf == t);
 }
 
-/// Floating-point numbers that represent integers scaled by negative
-/// powers of ten are inexact. For example, a premium rate of $2.40
-/// per $1000 is notionally 0.0024, but to the hardware may look like:
-///   0.0023999999999999998 [0x3ff69d495182a9930800]
-/// Multiplying that number by a million dollars and rounding down to
-/// cents yields 2399.99, where 2400.00 is wanted.
-///
-/// SOMEDAY !! The best way to handle this is to store integers. For
-/// the time being, multiplying by 1 + LDBL_EPSILON in problematic
-/// circumstances avoids this embarrassment while introducing an error
-/// that shouldn't matter.
-
-inline double ldbl_eps_plus_one_times(double z)
-{
-    static long double const y = 1.0L + std::numeric_limits<long 
double>::epsilon();
-    return static_cast<double>(y * z);
-}
-
 #endif // ieee754_hpp
diff --git a/ul_utilities.cpp b/ul_utilities.cpp
index 2c521232..625eca6a 100644
--- a/ul_utilities.cpp
+++ b/ul_utilities.cpp
@@ -23,11 +23,14 @@
 
 #include "ul_utilities.hpp"
 
+#include "assert_lmi.hpp"
+#include "bourn_cast.hpp"
 #include "calendar_date.hpp"
-#include "ieee754.hpp"                  // ldbl_eps_plus_one_times()
+#include "materially_equal.hpp"
 #include "round_to.hpp"
 
 #include <algorithm>                    // generate(), min()
+#include <cfenv>                        // fesetround()
 #include <cmath>                        // pow()
 #include <numeric>                      // inner_product()
 #include <vector>
@@ -98,5 +101,73 @@ currency max_modal_premium
     ,round_to<double> const& rounder
     )
 {
-    return rounder.c(ldbl_eps_plus_one_times(specamt * rate / mode));
+    // Assume premium-rate argument is precise to at most eight decimals,
+    // any further digits being representation error.
+    constexpr int radix {100'000'000};
+    // Premium rate and specified amount are nonnegative by their nature.
+    LMI_ASSERT(0.0 <= rate);
+    LMI_ASSERT(C0  <= specamt);
+    // Do not save and restore prior rounding direction, because lmi
+    // generally expects rounding to nearest everywhere.
+    std::fesetround(FE_TONEAREST);
+    // Make 'rate' a shifted integer.
+    // Shift the decimal point eight places, discarding anything further.
+    // Store the result as a wide integer, to be used in integer math.
+    // Use bourn_cast<>() for conversions here and elsewhere: it
+    // implicitly asserts that values are preserved.
+    std::int64_t irate = bourn_cast<std::int64_t>(std::nearbyint(rate * 
radix));
+    // If the rate really has more than eight significant (non-erroneous)
+    // digits, then treat them all as significant. In that case, there
+    // is no representation error to be removed. Here, 'tol' is just a
+    // guess; it may need refinement.
+    constexpr double tol = 1.0e-12;
+    if(!materially_equal(bourn_cast<double>(irate), rate * radix, tol))
+        {
+#if 0
+        // Enable this (including <iostream>) for research.
+        std::cout.precision(21);
+        std::cout
+            << "Excessive precision in rate; check the table\n"
+            << bourn_cast<double>(irate) << " bourn_cast<double>(irate)\n"
+            << rate * radix << " rate * radix\n"
+            << std::flush
+            ;
+#endif // 0
+        return rounder.c(specamt * rate / mode);
+        }
+#if 0
+    // Enable this assertion, adjusting the tolerance (last) argument
+    // p.r.n., if no table is allowed to have more than eight decimals.
+    LMI_ASSERT(materially_equal(bourn_cast<double>(irate), rate * radix, tol));
+#endif // 0
+    // Multiply integer rate by integral-cents specamt.
+    // Use a large integer type to avoid overflow.
+    std::int64_t iprod = irate * bourn_cast<std::int64_t>(specamt.cents());
+    // Result is an integer--safe to represent as double now.
+    // Function from_cents() has its own value-preservation test.
+    currency cprod = from_cents(bourn_cast<double>(iprod));
+    // Unshift the result, and round it in the specified direction.
+    // Dividing two integers generally yields a nonzero remainder,
+    // in which case do the division in floating point and round its
+    // result. However, if the remainder of integer division is zero,
+    // then the result is exact, in which case the corresponding
+    // rounded floating-point division may give the wrong answer.
+    std::int64_t quotient  = iprod / radix;
+    std::int64_t remainder = iprod % radix;
+    currency const annual_premium =
+        ((0 == remainder)
+        ? from_cents(bourn_cast<double>(quotient))
+        : rounder.c(cprod / radix)
+        );
+    // Calculate modal premium from annual as a separate step,
+    // using integer division to discard any fractional part.
+    // In a sense, this is double rounding, which is often a
+    // mistake, but here it's correct: the invariant
+    //   mode * max_modal_premium <= max_annual premium
+    // is explicitly desired. For example, if the maximum annual
+    // premium is 12.30, then the monthly maximum is 1.02,
+    // which is the highest level premium that can be paid twelve
+    // times without exceeding the annual maximum: 12.24 <= 12.30 .
+    std::int64_t annual_int = 
static_cast<std::int64_t>(annual_premium.cents());
+    return from_cents(bourn_cast<double>(annual_int / mode));
 }
diff --git a/ul_utilities_test.cpp b/ul_utilities_test.cpp
index 2a094983..51a064eb 100644
--- a/ul_utilities_test.cpp
+++ b/ul_utilities_test.cpp
@@ -23,9 +23,97 @@
 
 #include "ul_utilities.hpp"
 
+#include "materially_equal.hpp"
+#include "round_to.hpp"
 #include "test_tools.hpp"
 
+void test_max_modal_premium()
+{
+    round_to<double> const round_down(2, r_downward);
+    round_to<double> const round_near(2, r_to_nearest);
+    round_to<double> const round_not (2, r_not_at_all);
+    round_to<double> const round_up  (2, r_upward);
+
+    double   const rate    {0.0123456700000001};
+    currency const specamt {9'876'543'21_cents};
+
+    // This affects diagnostics shown when LMI_TEST_EQUAL() fails.
+    std::cout.precision(21);
+
+    LMI_TEST(materially_equal(12193254.3211401, rate * specamt.cents()));
+
+    currency p01 = max_modal_premium(rate, specamt, mce_annual , round_down);
+    LMI_TEST_EQUAL(121'932'54_cents, p01);
+    currency p02 = max_modal_premium(rate, specamt, mce_annual , round_near);
+    LMI_TEST_EQUAL(121'932'54_cents, p02);
+    currency p03 = max_modal_premium(rate, specamt, mce_annual , round_up  );
+    LMI_TEST_EQUAL(121'932'55_cents, p03);
+
+    LMI_TEST(materially_equal(10161.045267617, rate * specamt.cents() / 1200));
+    // Annual premium 'p01' is already rounded down to cents.
+    // Monthly premium is derived from annual.
+    LMI_TEST(materially_equal(10161.0450, p01 / 12));
+
+    currency p04 = max_modal_premium(rate, specamt, mce_monthly, round_down);
+    LMI_TEST_EQUAL( 10'161'04_cents, p04);
+    currency p05 = max_modal_premium(rate, specamt, mce_monthly, round_near);
+    LMI_TEST_EQUAL( 10'161'04_cents, p05);
+    // Rounding direction pertains to annual, not monthly.
+    // Monthly is always rounded down, to preserve the
+    //   12 * monthly <= annual
+    // invariant. Therefore, instead of
+    //   X/12, rounded up,
+    // this is
+    //   (X, rounded down) / 12, discarding the remainder.
+    currency p06 = max_modal_premium(rate, specamt, mce_monthly, round_up  );
+    LMI_TEST_EQUAL( 10'161'04_cents, p06);
+
+    // Real-world examples from system test.
+
+    currency q00 = max_modal_premium
+        (0.0195527999999999986536, 1'000'000'00_cents, mce_annual , 
round_down);
+    LMI_TEST_EQUAL(19'552'80_cents, q00);
+
+    currency q01 = max_modal_premium
+        (0.0195527999999999986536, 2'000'000'00_cents, mce_annual , 
round_down);
+    LMI_TEST_EQUAL(39'105'60_cents, q01);
+
+    currency q02 = max_modal_premium
+        (0.0193523999999999987698, 1'000'000'00_cents, mce_monthly, 
round_down);
+    LMI_TEST_EQUAL( 1'612'70_cents, q02);
+
+    currency q03 = max_modal_premium
+        (0.0128891999999999999627,   500'000'00_cents, mce_monthly, 
round_down);
+    LMI_TEST_EQUAL(   537'05_cents, q03);
+
+    currency q04 = max_modal_premium
+        (0.0128891999999999999627, 1'000'000'00_cents, mce_monthly, 
round_down);
+    LMI_TEST_EQUAL( 1'074'10_cents, q04);
+
+    currency q05 = max_modal_premium
+        (0.0105983999999999991409,    50'000'00_cents, mce_annual , 
round_down);
+    LMI_TEST_EQUAL(   529'92_cents, q05);
+
+    currency q06 = max_modal_premium
+        (0.0169656000000000008188,   250'000'00_cents, mce_annual , 
round_down);
+    LMI_TEST_EQUAL( 4'241'40_cents, q06);
+
+    currency q07 = max_modal_premium
+        (0.0169656000000000008188,   250'000'00_cents, mce_monthly, 
round_down);
+    LMI_TEST_EQUAL(   353'45_cents, q07);
+
+    currency q08 = max_modal_premium
+        (0.0169656000000000008188, 1'000'000'00_cents, mce_annual , 
round_down);
+    LMI_TEST_EQUAL(16'965'60_cents, q08);
+
+    currency q09 = max_modal_premium
+        (0.0382740000000000024638, 2'100'000'00_cents, mce_annual , 
round_down);
+    LMI_TEST_EQUAL(80'375'40_cents, q09);
+}
+
 int test_main(int, char*[])
 {
+    test_max_modal_premium();
+
     return EXIT_SUCCESS;
 }