JuliaMath / openlibm

High quality system independent, portable, open source libm implementation
https://openlibm.org
Other
515 stars 140 forks source link

wrong hypotl result near underflow #224

Closed zimmermann6 closed 3 years ago

zimmermann6 commented 3 years ago

with the attached program I get with OpenLibm 0.7.4:

$ gcc -fno-builtin test_hypot_openlibm.c /localdisk/zimmerma/openlibm-0.7.4/libopenlibm.a $ ./a.out x=0x6.ca4cf2be82edcp-16385 y=0xf.ffffffffffff8p-16354 z=0x0p+0

whereas the GNU libc gives:

$ ./a.out x=0x6.ca4cf2be82edcp-16385 y=0xf.ffffffffffff8p-16354 z=0xf.ffffffffffff8p-16354

Obviously hypot(x,y) >= |y|, and thus cannot be zero.

test_hypot_openlibm.c.gz

kargl commented 3 years ago

This patch fixes the issue of FreeBSD, which is the source of e_hypothl.c. t1 is intended to scale subnormal numbers, but t1=1 is needed instead of t1=0 when preparing the scaling. Don't know how to generate a pull request.

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;
zimmermann6 commented 3 years ago

I confirm this fixes the issue, thanks!