Re: [lmi] Numerics

[Vadim Zeitlin]

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

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?

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.
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.
VZ