[PACTH,libm] hypothl(x) mishandles subnormal numbers.
Steve Kargl
sgk at troutmask.apl.washington.edu
Wed Feb 10 01:42:17 UTC 2021
On Wed, Feb 10, 2021 at 12:26:29AM +0100, Dimitry Andric wrote:
>
> > On 10 Feb 2021, at 00:15, Steve Kargl <sgk at troutmask.apl.washington.edu> wrote:
> >
> > This patch fixes the issue. t1 is used to scale the subnormal
> > numbers. It is generated by scaling the exponent, but that
> > only works if t1 is 1 not 0.
> >
> > Index: src/e_hypotl.c
> > ===================================================================
> > --- src/e_hypotl.c (revision 2342)
> > +++ src/e_hypotl.c (working copy)
> > @@ -82,7 +82,7 @@ hypotl(long double x, long double y)
> > man_t manh, manl;
> > GET_LDBL_MAN(manh,manl,b);
> > if((manh|manl)==0) return a;
> > - t1=0;
> > + t1=1;
> > SET_HIGH_WORD(t1,ESW(MAX_EXP-2)); /* t1=2^(MAX_EXP-2) */
> > b *= t1;
> > a *= t1;
> >
>
> Ah, having looked at the glibc code, I had concluded something similar
> in https://bugs.freebsd.org/bugzilla/show_bug.cgi?id=253313#c2, but
> using INSERT_LDBL80_WORDS(t1,ESW(MAX_EXP-2),0x8000000000000000).
>
> Your solution is a lot simpler though. :) Note that to make it work, I
> also needed to insert a volatile into the INSERT_LDBL80_WORDS() macro.
I don't look at glibc's libm code due to possible Copyright
issues (long/short story that I rather not get into here).
I do, however, have a math_sgk.h that appears at the
end of math_privated.h with a bunch of instrumentation code
hidden behind -DEBUG. For the above, it becomes
t1=0;
SET_HIGH_WORD(t1,ESW(MAX_EXP-2)); /* t1=2^(MAX_EXP-2) */
PRNL(t1);
which yields inf as output. Clearly, not a scaling of a subnormal
to a normal number.
> There are more places where this is apparently needed, due to the way
> recent clang versions tend to over-optimize floating point code at
> compile time. And specifically for the case where one union field is
> written, and then another field read, sometimes leading to unexpected
> results...
Hmmm. This is a bummer. I suppose someone will say the code
in msun is using undefined behavior or some such standardese.
--
Steve
More information about the freebsd-current
mailing list