[Date Prev][Date Next][Thread Prev][Thread Next][Date Index][Thread Index]
[lmi-commits] [lmi] master 205497ae 8/9: Reimplement max_modal_premium()
From: |
Greg Chicares |
Subject: |
[lmi-commits] [lmi] master 205497ae 8/9: Reimplement max_modal_premium() |
Date: |
Fri, 6 May 2022 19:37:36 -0400 (EDT) |
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;
}
- [lmi-commits] [lmi] master updated (13799952 -> 9e5a09bf), Greg Chicares, 2022/05/06
- [lmi-commits] [lmi] master be841c05 3/9: Refactor for clarity, Greg Chicares, 2022/05/06
- [lmi-commits] [lmi] master a24d665c 4/9: Include only appropriate headers, and say why they're included, Greg Chicares, 2022/05/06
- [lmi-commits] [lmi] master 46534540 2/9: Refactor for uniformity, Greg Chicares, 2022/05/06
- [lmi-commits] [lmi] master 9e5a09bf 9/9: Record speed measurements, Greg Chicares, 2022/05/06
- [lmi-commits] [lmi] master 78657d01 6/9: Improve physical structure, Greg Chicares, 2022/05/06
- [lmi-commits] [lmi] master 4b8bf312 1/9: Expunge workarounds for nonexistent problems, Greg Chicares, 2022/05/06
- [lmi-commits] [lmi] master 67af2e9b 5/9: Improve a debatable workaround, Greg Chicares, 2022/05/06
- [lmi-commits] [lmi] master 205497ae 8/9: Reimplement max_modal_premium(),
Greg Chicares <=
- [lmi-commits] [lmi] master 1de83b9d 7/9: Refactor, Greg Chicares, 2022/05/06