Date: Tue, 9 Feb 2021 17:42:14 -0800 From: Steve Kargl <sgk@troutmask.apl.washington.edu> To: Dimitry Andric <dim@freebsd.org> Cc: FreeBSD Current <freebsd-current@freebsd.org> Subject: Re: [PACTH,libm] hypothl(x) mishandles subnormal numbers. Message-ID: <20210210014214.GA60428@troutmask.apl.washington.edu> In-Reply-To: <051767C6-85FB-48E8-AFE1-546DC41E8D17@FreeBSD.org> References: <20210206203929.GA45801@troutmask.apl.washington.edu> <20210206210448.GC45801@troutmask.apl.washington.edu> <C5652F7E-A943-42AA-A70B-1D4C1C763E06@FreeBSD.org> <20210209231529.GA59891@troutmask.apl.washington.edu> <051767C6-85FB-48E8-AFE1-546DC41E8D17@FreeBSD.org>
next in thread | previous in thread | raw e-mail | index | archive | help
On Wed, Feb 10, 2021 at 12:26:29AM +0100, Dimitry Andric wrote: > > > On 10 Feb 2021, at 00:15, Steve Kargl <sgk@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
Want to link to this message? Use this URL: <https://mail-archive.FreeBSD.org/cgi/mid.cgi?20210210014214.GA60428>