Re: [PACTH,libm] hypothl(x) mishandles subnormal numbers.

2021-02-09 Thread Steve Kargl 
On Wed, Feb 10, 2021 at 12:26:29AM +0100, Dimitry Andric wrote:
> 
> > On 10 Feb 2021, at 00:15, Steve Kargl  
> > 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),0x8000).
> 
> 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
___
freebsd-current@freebsd.org mailing list
https://lists.freebsd.org/mailman/listinfo/freebsd-current
To unsubscribe, send any mail to "freebsd-current-unsubscr...@freebsd.org"

Re: [PACTH,libm] hypothl(x) mishandles subnormal numbers.

2021-02-09 Thread Dimitry Andric 


> On 10 Feb 2021, at 00:15, Steve Kargl  
> wrote:
> 
> On Sat, Feb 06, 2021 at 10:32:33PM +0100, Dimitry Andric wrote:
>> On 6 Feb 2021, at 22:04, Steve Kargl  
>> wrote:
>>> 
>>> On Sat, Feb 06, 2021 at 12:39:29PM -0800, Steve Kargl wrote:
 I've long forgotten by freebsd bugzilla password.
 So, if someone would like to submit a bug report,
 here's a test program.
 
>>> Forgot to include that issue was identified from
>>> a bug report in the OpenLibm bug mailing list.
>>> 
>>> https://github.com/JuliaMath/openlibm/issues/224
>> 
>> I put this in . Now the trick is to
>> figure out what is going on in e_hypotl.c... :)
>> 
> 
> 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),0x8000).

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.

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...

In any case, I will apply this soon, and also add a test case. Thanks!

-Dimitry



signature.asc
Description: Message signed with OpenPGP

[PACTH,libm] hypothl(x) mishandles subnormal numbers.

2021-02-09 Thread Steve Kargl 
On Sat, Feb 06, 2021 at 10:32:33PM +0100, Dimitry Andric wrote:
> On 6 Feb 2021, at 22:04, Steve Kargl  
> wrote:
> > 
> > On Sat, Feb 06, 2021 at 12:39:29PM -0800, Steve Kargl wrote:
> >> I've long forgotten by freebsd bugzilla password.
> >> So, if someone would like to submit a bug report,
> >> here's a test program.
> >> 
> > Forgot to include that issue was identified from
> > a bug report in the OpenLibm bug mailing list.
> > 
> > https://github.com/JuliaMath/openlibm/issues/224
> 
> I put this in . Now the trick is to
> figure out what is going on in e_hypotl.c... :)
> 

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;


-- 
Steve
___
freebsd-current@freebsd.org mailing list
https://lists.freebsd.org/mailman/listinfo/freebsd-current
To unsubscribe, send any mail to "freebsd-current-unsubscr...@freebsd.org"

Re: hypothl(x) mishandles subnormal numbers.

2021-02-06 Thread Steve Kargl 
On Sat, Feb 06, 2021 at 10:32:33PM +0100, Dimitry Andric wrote:
> On 6 Feb 2021, at 22:04, Steve Kargl  
> wrote:
> > 
> > On Sat, Feb 06, 2021 at 12:39:29PM -0800, Steve Kargl wrote:
> >> I've long forgotten by freebsd bugzilla password.
> >> So, if someone would like to submit a bug report,
> >> here's a test program.
> >> 
> > Forgot to include that issue was identified from
> > a bug report in the OpenLibm bug mailing list.
> > 
> > https://github.com/JuliaMath/openlibm/issues/224
> 
> I put this in . Now the trick is to
> figure out what is going on in e_hypotl.c... :)
> 

Thanks.  I took a quick look, and there seems to
be some magic scaling happen.  hypotl(x,y) should
compute sqrtl(x*x+y*y) while avoiding overflows
and underflows.  I suspect that the scaling needs
to be tweaked.  It may take me a bit to work out
how to fix it.

-- 
Steve
___
freebsd-current@freebsd.org mailing list
https://lists.freebsd.org/mailman/listinfo/freebsd-current
To unsubscribe, send any mail to "freebsd-current-unsubscr...@freebsd.org"

Re: hypothl(x) mishandles subnormal numbers.

2021-02-06 Thread Dimitry Andric 
On 6 Feb 2021, at 22:04, Steve Kargl  wrote:
> 
> On Sat, Feb 06, 2021 at 12:39:29PM -0800, Steve Kargl wrote:
>> I've long forgotten by freebsd bugzilla password.
>> So, if someone would like to submit a bug report,
>> here's a test program.
>> 
> Forgot to include that issue was identified from
> a bug report in the OpenLibm bug mailing list.
> 
> https://github.com/JuliaMath/openlibm/issues/224

I put this in . Now the trick is to
figure out what is going on in e_hypotl.c... :)

-Dimitry



signature.asc
Description: Message signed with OpenPGP

Re: hypothl(x) mishandles subnormal numbers.

2021-02-06 Thread Steve Kargl 
On Sat, Feb 06, 2021 at 12:39:29PM -0800, Steve Kargl wrote:
> I've long forgotten by freebsd bugzilla password.  
> So, if someone would like to submit a bug report,
> here's a test program.
> 
Forgot to include that issue was identified from
a bug report in the OpenLibm bug mailing list.

https://github.com/JuliaMath/openlibm/issues/224

-- 
Steve
___
freebsd-current@freebsd.org mailing list
https://lists.freebsd.org/mailman/listinfo/freebsd-current
To unsubscribe, send any mail to "freebsd-current-unsubscr...@freebsd.org"

hypothl(x) mishandles subnormal numbers.

2021-02-06 Thread Steve Kargl 
I've long forgotten by freebsd bugzilla password.  
So, if someone would like to submit a bug report,
here's a test program.

% cc -o z -O2 -fno-builtin z.c -lm
% ./z
Via scaling and sqrtl
x=2.853684e-4932 y=1.444012e-4922 a=1.444012e-4922
x=0x1.b2933cafa0bb8p-16383 y=0x1.fp-16351 
a=0x1.6ca4c8p-16350

Via hypotl
x=2.853684e-4932 y=1.444012e-4922 z=nan
x=0x1.b2933cafa0bb8p-16383 y=0x1.fp-16351 z=nan

% cat z.c
#include 
#include 

#if defined(__i386__)
#include 
#endif

int main()
{
  long double x, y, z, a;

#if defined(__i386__)
fpsetprec(FP_PE);
#endif

   x = 0x3.6526795f4176ep-16384L;
   y = 0x3.ep-16352L;
   z = hypotl(x, y);

   if (x > y) {
  a = y / x;
  a = x * sqrtl(1 + a);
   } else {
  a = x / y;
  a = y * sqrtl(a + 1);
   }

   printf("Via scaling and sqrtl\n");
   printf("x=%Le y=%Le a=%Le\n", x, y, a);
   printf("x=%La y=%La a=%La\n", x, y, a);

   printf("\nVia hypotl\n");
   printf("x=%Le y=%Le z=%Le\n", x, y, z);
   printf("x=%La y=%La z=%La\n", x, y, z);

   return 0;
}

-- 
Steve
___
freebsd-current@freebsd.org mailing list
https://lists.freebsd.org/mailman/listinfo/freebsd-current
To unsubscribe, send any mail to "freebsd-current-unsubscr...@freebsd.org"