Re: [lmi] Numerics

[Greg Chicares]

On 2016-03-27 22:41, Vadim Zeitlin wrote:
> On Thu, 24 Mar 2016 21:20:54 +0000 Greg Chicares <address@hidden> wrote:
> 
> GC> On 2016-03-18 00:56, Vadim Zeitlin wrote:
> GC> [...]
> GC> > [*] See 
> https://randomascii.wordpress.com/2013/02/07/float-precision-revisited-nine-digit-float-portability/
> GC> 
> GC> And now the expression
> GC>   0.07 * (1.0 + 1.0 * epsilon)
> GC> doesn't seem to equal
> GC>   0.07 × (1 + ε)
> GC> with gcc-4.9.1 ; see:
> GC>   http://lists.nongnu.org/archive/html/lmi-commits/2016-03/msg00008.html
> GC> for a unit test that always worked with gcc-3.4.5, but failed with 
> gcc-4.9.1
> GC> until I substituted a literal for that expression. I doubt gcc is wrong.
> 
>  Yes, I think that gcc (4.9) is indeed correct and it's the test itself
> that was wrong. I don't really understand where did the idea that
> multiplying by "1.0 + epsilon" (I discard the "1.0*" before epsilon as it's
> completely optimized away by the compiler anyhow) is supposed to give the
> "next" floating point number come from, but it's clearly not true in
> practice, at least when using x87.

Here's the motivation. Some rate table contains 0.01130, and we multiply
that by some number like 1375000, which would yield 15537.50 in exact
decimal math. But in binary math, the rate has no exact representation.
One of the two closest representations must be used. One is slightly
higher, and the other slightly lower. The result is to be rounded down
to cents. The hope was that
  0.01130 * 1375000 * (1 + ε)
would give a value no less than the exact decimal number 15537.50, and
no more than that number times (1 + 2ε), so that rounding it down to
cents would yield 15537.50 and not 15537.49 . In this particular
example, BasicValues::GetModalPremTgtFromTable() does this:

    return round_max_premium()
        (
            (   TgtPremMonthlyPolFee * 12
            +       a_specamt
                *   ldbl_eps_plus_one()
                *   MortalityRates_->TargetPremiumRates()[0]
            )
        /   a_mode
        );

yet it yields the unwanted 15537.49 .

So this trick to obtain the "next" floating-point number doesn't work
in this actual case that we encountered in regression testing today.
The ideal solution is to use integer math: 1130 * 1375000 produces an
exact integer result. We just haven't had time to do rewrite these
calculations in terms of the currency class discussed earlier.