From 31dc8730af585d8e13021484752fb20decae0661 Mon Sep 17 00:00:00 2001 From: Adhemerval Zanella Date: Fri, 4 May 2012 13:05:57 +0200 Subject: Fix for ldbl-128ibm acosl/asinl inaccuracies 2012-05-02 Adhemerval Zanella * sysdeps/ieee754/ldbl-128ibm/e_acosl.c (__ieee754_acosl): Fix long double comparison inaccuracies. * sysdeps/ieee754/ldbl-128ibm/e_asinl.c (__ieee754_asinl): * Likewise. * sysdeps/powerpc/fpu/libm-test-ulps: Update. --- sysdeps/ieee754/ldbl-128ibm/e_asinl.c | 31 ++++++++++++------------------- 1 file changed, 12 insertions(+), 19 deletions(-) (limited to 'sysdeps/ieee754/ldbl-128ibm/e_asinl.c') diff --git a/sysdeps/ieee754/ldbl-128ibm/e_asinl.c b/sysdeps/ieee754/ldbl-128ibm/e_asinl.c index fb6f572..b395439 100644 --- a/sysdeps/ieee754/ldbl-128ibm/e_asinl.c +++ b/sysdeps/ieee754/ldbl-128ibm/e_asinl.c @@ -132,25 +132,18 @@ long double __ieee754_asinl (long double x) { long double t, w, p, q, c, r, s; - int32_t ix, sign, flag; + int flag; ieee854_long_double_shape_type u; flag = 0; - u.value = x; - sign = u.parts32.w0; - ix = sign & 0x7fffffff; - u.parts32.w0 = ix; /* |x| */ - if (ix >= 0x3ff00000) /* |x|>= 1 */ + u.value = __builtin_fabsl (x); + if (u.value == 1.0L) /* |x|>= 1 */ + return x * pio2_hi + x * pio2_lo; /* asin(1)=+-pi/2 with inexact */ + else if (u.value >= 1.0L) + return (x - x) / (x - x); /* asin(|x|>1) is NaN */ + else if (u.value < 0.5L) { - if (ix == 0x3ff00000 - && (u.parts32.w1 | (u.parts32.w2 & 0x7fffffff) | u.parts32.w3) == 0) - /* asin(1)=+-pi/2 with inexact */ - return x * pio2_hi + x * pio2_lo; - return (x - x) / (x - x); /* asin(|x|>1) is NaN */ - } - else if (ix < 0x3fe00000) /* |x| < 0.5 */ - { - if (ix < 0x3c600000) /* |x| < 2**-57 */ + if (u.value < 6.938893903907228e-18L) /* |x| < 2**-57 */ { if (huge + x > one) return x; /* return x with inexact if x!=0 */ @@ -162,7 +155,7 @@ __ieee754_asinl (long double x) flag = 1; } } - else if (ix < 0x3fe40000) /* 0.625 */ + else if (u.value < 0.625L) { t = u.value - 0.5625; p = ((((((((((rS10 * t @@ -189,7 +182,7 @@ __ieee754_asinl (long double x) + sS1) * t + sS0; t = asinr5625 + p / q; - if ((sign & 0x80000000) == 0) + if (x > 0.0L) return t; else return -t; @@ -230,7 +223,7 @@ __ieee754_asinl (long double x) } s = __ieee754_sqrtl (t); - if (ix >= 0x3fef3333) /* |x| > 0.975 */ + if (u.value > 0.975L) { w = p / q; t = pio2_hi - (2.0 * (s + s * w) - pio2_lo); @@ -248,7 +241,7 @@ __ieee754_asinl (long double x) t = pio4_hi - (p - q); } - if ((sign & 0x80000000) == 0) + if (x > 0.0L) return t; else return -t; -- cgit v1.1