Date: Sun, 7 Feb 2021 14:53:04 -0800 From: Steve Kargl <sgk@troutmask.apl.washington.edu> To: freebsd-current@freebsd.org Subject: [PATCH, LIBM] powf is wrong with x near 1 and |y| >> 1 Message-ID: <20210207225304.GC50688@troutmask.apl.washington.edu>
next in thread | raw e-mail | index | archive | help
See the last few comments here: https://github.com/JuliaMath/openlibm/issues/211 Test program: #include <stdio.h> #include <math.h> #if defined(__i386__) #include <ieeefp.h> #endif int main() { float x, y, z, a; #if defined(__i386__) fpsetprec(FP_PE); #endif x = 0x1.ffffeep-1f; y = -0x1.000002p+27f; z = powf (x, y); printf ("x=%e y=%e z=%e\n", x, y, z); printf ("x=%a y=%a z=%a\n", x, y, z); a = expf(y * logf(x)); // 0x1.d53532p+103f; printf ("floats: a=%a a=%e\n", a, a); a = (float)exp(y * log(x)); // 0x1.d53532p+103f; printf (" dbles: a=%a a=%e\n", a, a); return 0; } Patch without too much analysis for x near 1 or |y| > 2**27. Index: src/e_powf.c =================================================================== --- src/e_powf.c (revision 2342) +++ src/e_powf.c (working copy) @@ -136,7 +136,7 @@ __ieee754_powf(float x, float y) /* |y| is huge */ if(iy>0x4d000000) { /* if |y| > 2**27 */ /* over/underflow if x is not close to one */ - if(ix<0x3f7ffff8) return (hy<0)? sn*huge*huge:sn*tiny*tiny; + if(ix<0x3f7ffff7) return (hy<0)? sn*huge*huge:sn*tiny*tiny; if(ix>0x3f800007) return (hy>0)? sn*huge*huge:sn*tiny*tiny; /* now |1-x| is tiny <= 2**-20, suffice to compute log(x) by x-x^2/2+x^3/3-x^4/4 */ -- Steve
Want to link to this message? Use this URL: <https://mail-archive.FreeBSD.org/cgi/mid.cgi?20210207225304.GC50688>