Re: [lmi] Numerics

[Vadim Zeitlin]

On Wed, 14 Dec 2016 02:30:04 +0000 Greg Chicares <address@hidden> wrote:

GC> On 2016-12-14 00:52, Vadim Zeitlin wrote:
GC> > On Tue, 13 Dec 2016 15:57:34 +0000 Greg Chicares <address@hidden> wrote:
GC> > 
GC> > GC> On 2016-03-27 22:41, Vadim Zeitlin wrote:
GC> > GC> > On Thu, 24 Mar 2016 21:20:54 +0000 Greg Chicares <address@hidden> 
wrote:
GC> > GC> > 
GC> > GC> > GC> On 2016-03-18 00:56, Vadim Zeitlin wrote:
GC> [...an apparently feckless attempt to avoid some rounding problems...]
GC> > GC> 
GC> > GC> Here's the motivation. Some rate table contains 0.01130, and we 
multiply
GC> > GC> that by some number like 1375000, which would yield 15537.50 in exact
GC> > GC> decimal math. But in binary math, the rate has no exact 
representation.
GC> > 
GC> >  Yes, I understand this.
GC> 
GC> I wrote that only in case Kim wanted a mailing-list post to refer to in
GC> documentation that's potentially read by auditors who don't understand it.
GC> I didn't mean to imply that you were unaware of this.

 I mostly wrote the above to show the point where I stopped following, i.e.
for dramatic contrast with "I don't understand this" below.

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

 Yes, again, I agree until now.

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

 Yes, but this is not "natural": you could naïvely expect the successor of
2 to be 2+ε instead. This doesn't work, however, because of precision loss
and 2+ε turns out to the same as 2. 2+2ε does work because it has the same
significand as 1+ε and just uses a greater exponent.

GC> And the successor of any general floating-point number d would be d+dε, no?

 This is where I stop understanding. What possible reason would be for this
to be true in general?

 To hopefully settle this issue once and for all, let me notice that this
is already not true for d=3: 3+3ε == 3*(1+ε) == 0x1.8000000000002p+1 in 64
bit IEEE-754 representation, but this is, obviously, not the successor of 3
which is -- as returned by nextafter() -- 0x1.8000000000001p+1.

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

 I don't know. I don't see what is right with it, i.e. I simply don't see
any reason for this to be true. And as it's demonstrably indeed not true, I
think it's not really as obvious as you think...


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

 This would definitely be the right solution. Looking at the headers, it
seems like nextafter() is implemented in terms of a compiler intrinsic
function, so it should work just fine. But even if it doesn't, once again,
it would be better to just increase the representation of the double as a
64 bit (unsigned) integer instead of trying to multiply it by something (of
course, special values would need to be carefully accounted for in this
case).

 Thanks in advance,
VZ