Re: [lmi] Numerics

[Greg Chicares]

On 2016-12-14 00:52, Vadim Zeitlin wrote:
> On Tue, 13 Dec 2016 15:57:34 +0000 Greg Chicares <address@hidden> wrote:
> 
> GC> On 2016-03-27 22:41, Vadim Zeitlin wrote:
> GC> > On Thu, 24 Mar 2016 21:20:54 +0000 Greg Chicares <address@hidden> wrote:
> GC> > 
> GC> > GC> On 2016-03-18 00:56, Vadim Zeitlin wrote:
[...an apparently feckless attempt to avoid some rounding problems...]
> GC> 
> GC> Here's the motivation. Some rate table contains 0.01130, and we multiply
> GC> that by some number like 1375000, which would yield 15537.50 in exact
> GC> decimal math. But in binary math, the rate has no exact representation.
> 
>  Yes, I understand this.

I wrote that only in case Kim wanted a mailing-list post to refer to in
documentation that's potentially read by auditors who don't understand it.
I didn't mean to imply that you were unaware of this.

> GC> One of the two closest representations must be used. One is slightly
> GC> higher, and the other slightly lower. The result is to be rounded down
> GC> to cents. The hope was that
> GC>   0.01130 * 1375000 * (1 + ε)
> GC> would give a value no less than the exact decimal number 15537.50, and
> GC> no more than that number times (1 + 2ε),
> 
>  But I don't understand what is this hope based on. I may be missing
> something, but I just don't see any reason for this to be true... To go to
> the next representable floating point number we need to increase the last
> bit of the significand, but why would multiplying by 1+ε do that?

If we want the next higher representable floating point number after unity,
that's 1+ε by definition; so multiplying 1.0 by 1+ε give the correct answer.

If we want the successor of 2.0, that should be 2+2ε, right?

And the successor of any general floating-point number d would be d+dε, no?
But that's d*(1+ε), so why isn't this algorithm valid for every number with
the exception of zeros, NaNs, infinities, and maybe denormals?

It seems so obvious to me that I figure someone else must have thought of
it, but...
  https://www.google.com/search?q=%22one+plus+epsilon%22+multiply
...it looks like this might be my very own original idea--which worries me,
but I don't see what's essentially wrong with it.

> GC> So this trick to obtain the "next" floating-point number doesn't work
> GC> in this actual case that we encountered in regression testing today.

This observed failure doesn't disprove the idea behind the 1+ε algorithm.
It might be that I should have used DBL_EPSILON instead of LDBL_EPSILON.
(I just figured that, in every case where this technique is used, all
operands would be on the x87 stack, so LDBL_EPSILON would suffice. But
I guess popping off the stack to a double could return the next-lower
double.)

> GC> The ideal solution is to use integer math: 1130 * 1375000 produces an
> GC> exact integer result.
> 
>  Yes, absolutely. But if we just want to get the next representable
> floating point number we could just use the standard (in C++11) nextafter()
> function or, in case we have any QoI problems with it, do it ourselves
> easily enough by bit fiddling.
> 
>  But, ceterum autem censeo, ldbl_eps_plus_one() esse delendam.

The only reason I didn't use nextafter() is that it didn't exist when
this code was originally written. But now we're using no old compiler
whose RTL would lack it, so there's no reason not to use it unless we
find a QoI issue.

I'd like to try nextafter() with the particular regression error we
found today and see what happens. I'll do that tomorrow morning.