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_acosl.c | 37 +++++++++++++++-------------------- sysdeps/ieee754/ldbl-128ibm/e_asinl.c | 31 ++++++++++++----------------- 2 files changed, 28 insertions(+), 40 deletions(-) (limited to 'sysdeps/ieee754') diff --git a/sysdeps/ieee754/ldbl-128ibm/e_acosl.c b/sysdeps/ieee754/ldbl-128ibm/e_acosl.c index 533b59786d..5d2af30346 100644 --- a/sysdeps/ieee754/ldbl-128ibm/e_acosl.c +++ b/sysdeps/ieee754/ldbl-128ibm/e_acosl.c @@ -152,30 +152,25 @@ long double __ieee754_acosl (long double x) { long double z, r, w, p, q, s, t, f2; - int32_t ix, sign; ieee854_long_double_shape_type u; - 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) + { + if (x > 0.0L) + return 0.0; /* acos(1) = 0 */ + else + return (2.0 * pio2_hi) + (2.0 * pio2_lo); /* acos(-1)= pi */ + } + else if (u.value > 1.0L) { - if (ix == 0x3ff00000 - && (u.parts32.w1 | (u.parts32.w2&0x7fffffff) | u.parts32.w3) == 0) - { /* |x| == 1 */ - if ((sign & 0x80000000) == 0) - return 0.0; /* acos(1) = 0 */ - else - return (2.0 * pio2_hi) + (2.0 * pio2_lo); /* acos(-1)= pi */ - } return (x - x) / (x - x); /* acos(|x| > 1) is NaN */ } - else if (ix < 0x3fe00000) /* |x| < 0.5 */ + if (u.value < 0.5L) { - if (ix < 0x3c600000) /* |x| < 2**-57 */ + if (u.value < 6.938893903907228e-18L) /* |x| < 2**-57 */ return pio2_hi + pio2_lo; - if (ix < 0x3fde0000) /* |x| < .4375 */ + if (u.value < 0.4375L) { /* Arcsine of x. */ z = x * x; @@ -229,13 +224,13 @@ __ieee754_acosl (long double x) + Q1) * t + Q0; r = p / q; - if (sign & 0x80000000) + if (x < 0.0L) r = pimacosr4375 - r; else r = acosr4375 + r; return r; } - else if (ix < 0x3fe40000) /* |x| < 0.625 */ + else if (u.value < 0.625L) { t = u.value - 0.5625L; p = ((((((((((rS10 * t @@ -261,7 +256,7 @@ __ieee754_acosl (long double x) + sS2) * t + sS1) * t + sS0; - if (sign & 0x80000000) + if (x < 0.0L) r = pimacosr5625 - p / q; else r = acosr5625 + p / q; @@ -309,7 +304,7 @@ __ieee754_acosl (long double x) + qS0; r = s + (w + s * p / q); - if (sign & 0x80000000) + if (x < 0.0L) w = pio2_hi + (pio2_lo - r); else w = r; diff --git a/sysdeps/ieee754/ldbl-128ibm/e_asinl.c b/sysdeps/ieee754/ldbl-128ibm/e_asinl.c index fb6f5720d0..b395439495 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 1.4.1