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,

It probably flowed from this line of thought:

    /// Like C99 nextafter(), but prevents range error, and returns
    /// exact integral values unchanged.

    template<typename T>
    T adjust_bound(T t, T direction)

where the implementation does something kind of like this:

    return static_cast<T> (t * (1 + std::numeric_limits<T>::epsilon()));

though obviously I didn't write the test as carefully as the implementation.

(Here's the motivation. If a textcontrol is to accept an interest rate in
a range like [0.03, 0.07], and treat input outside that range as invalid,
then it would be naive to compare the input directly to either bound: if
"0.07" is correctly entered and the machine translates that text string to
0.07000000000000001 (e.g.), and compares it to a bound represented as, say,
0.06999999999999999, then we shouldn't reject it, because it's the user's
best attempt at specifying the upper-limit interest rate. Therefore, the
bounds are "adjusted" outward by one ulp.)

But (if we take adequate precautions, like casting to the original type as
is done in the implementation (though formerly not in the test)), then I do
think that multiplying by (1 + ε) gives the "next" value of the given
floating-point type (except for special cases like +/-INFINITY or zero,
which are handled in the implementation though elided in the part loosely
quoted above).

Clearly C99 nextafter() would be preferable, but I think it was unavailable
with MinGW gcc-3.x at the time I wrote that. Yet is the "× (1 + ε)" trick
actually incorrect (except for special values as noted above)? Isn't it
generally true (modulo those exceptions) that nextafter(X) / X is 1 + ε?

>  But IMO, to find the next double, it would be better to rely on its
> representation mandated by IEEE-754 and do the usual trick with casting to
> int64_t and incrementing it and then casting back to double. This is
> guaranteed to work, unlike multiplying by epsilon. Of course, hard coding
> the values, as r6523 does, works too, but I'd get rid of the remaining uses
> of epsilon too.

I imagine that I preferred multiplying by one plus epsilon because it's
easier to read and write, and especially because it's harder to get wrong
than twiddling the bits. Of course, nextafter() has those advantages with
none of the disadvantages, and we do seem to have it in MinGW-w64 gcc now:

  i686-w64-mingw32/include/c++/tr1/math.h:using std::tr1::nextafter;

Getting rid of the remaining uses of epsilon would be good, including these:
  $grep 'epsilon_plus_one\|ldbl_eps_plus_one' *.?pp
which have the same rationale as described above.

More idealistically, we should probably use integral cents as our basic
currency unit instead of floating-point dollars rounded to the closest
approximation to integral hundredths, because in the real world we can
have exactly seven cents, and (double)(0.07) is not exactly the same.