Re: [PATCH, LIBM] powf is wrong with x near 1 and |y| >> 1

[Dimitry Andric]

On 7 Feb 2021, at 23:53, Steve Kargl  wrote:
> 
> See the last few comments here:
> 
> https://github.com/JuliaMath/openlibm/issues/211

Thanks Steve. Committed here:

https://cgit.freebsd.org/src/commit/?id=93fc67896550548f91b307dbe3053f11db5d4a8a

-Dimitry



signature.asc
Description: Message signed with OpenPGP

[PATCH, LIBM] powf is wrong with x near 1 and |y| >> 1

[Steve Kargl]

See the last few comments here:

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

Test program:

#include 
#include 
  
#if defined(__i386__)
#include 
#endif

int main()
{
  float x, y, z, a;
#if defined(__i386__)
fpsetprec(FP_PE);
#endif
  x = 0x1.eep-1f;
  y = -0x1.02p+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>0x4d00) { /* if |y| > 2**27 */
/* over/underflow if x is not close to one */
-   if(ix<0x3f78) return (hy<0)? sn*huge*huge:sn*tiny*tiny;
+   if(ix<0x3f77) return (hy<0)? sn*huge*huge:sn*tiny*tiny;
if(ix>0x3f87) 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
___
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"