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) _at__at_ -136,7 +136,7 _at__at_ __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 */ -- SteveReceived on Sun Feb 07 2021 - 21:53:12 UTC
This archive was generated by hypermail 2.4.0 : Wed May 19 2021 - 11:41:27 UTC