[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