Re: [lmi] How to deal with this instance of UB?

[Greg Chicares]

On 6/16/22 21:39, Vadim Zeitlin wrote:
> On Thu, 16 Jun 2022 20:30:09 +0000 Greg Chicares <gchicares@sbcglobal.net> 
> wrote:
> 
> GC> Vadim--The patch below shows two ways of dealing with UB at a
> GC> particular point in fdlibm; of course want to choose one way.
> GC> 
> GC> The first way is to use an __attribute__ that blinds UBSan to
> GC> the problem. I dislike doing that, except as a last resort
[...]
> GC> The second way is to use casts to replace UB with defined
> GC> behavior. That seems like the righteous path, and I think
> GC> I've followed it correctly below, but I'd rather not commit
> GC> this without your prior review.
> 
>  This is definitely not going to result in undefined behaviour, but I'm not
> really sure what is the intention of the current code. I.e. is it really
> expected that shifting a negative number left can result in a positive
> number?

AFAICT, the argument is transformed so that it lies within a
range where the polynomial approximation works well, and the
transformation multiplies it by the kth power of two, so the
answer is some inverse transformation applied to the value
obtained for the transformed argument.

Thus, at this stage, they just want to multiply by something
like 2^k, and they do that by bit-shifting an appropriate
subfield of the binary64 representation. I believe they're
thinking of 'k' as unsigned by that point, so that we're
just rewriting what they meant in a non-UB way.

I can't imagine I'm the only person to whom it ever occurred
to test log1p(-0.999999), or to call expm1() with an argument
about an decimal order of magnitude smaller than that. But
this test case is interesting because...

    //   https://www.wolframalpha.com/input?i=log1p(-0.999999)/12
    //   -1.151292546497022842008995727342182103800550744314386488016

    // In this ill-conditioned case, we get something like eleven
    // digits of precision--an error of about one million ulp.
    double const i0 = std::log1p(-0.999999) / 12.0;

    // i0 = -1.15129254649462642312585 i1 = -0.68377223398316200331237 MinGW-w32
    // i0 = -1.15129254649462642312585 i1 = -0.68377223398316200331237 glibc

The fdlibm authors claim to be within less than one ulp,
for each of exmp1() and log1p(). And glibc's results exactly
match MinGW-w32's, but IIRC MinGW-w32 uses the x87 for these
functions, while glibc uses fdlibm--though maybe MinGW-w32
detects extreme cases and falls back on some other technique
(which just might happen to be fdlibm).

                    1111  tens
          1 234567890123  digits
          - ------------
    //   -1.151292546497  wolfram
    //   -1.151292546494  glibc

That looks like a discrepancy in the thirteenth digit.
That seems to clash with the comment above, where I seemed
to have only had eleven matching digits--but maybe that was
after composing these two transcendental functions and
losing a thousand ulp in each, so a million in the
composed function?

> GC>  - Do you agree that the second way is better in principle?
> 
>  Yes, of course.
[...]
>  Sorry if you hoped for something more definitive

No, I was hoping only for a sanity check. Thanks.