This is the mail archive of the libc-alpha@sourceware.org mailing list for the glibc project.

Index Nav: [Date Index] [Subject Index] [Author Index] [Thread Index]
Message Nav: [Date Prev] [Date Next] [Thread Prev] [Thread Next]
Other format: [Raw text]

Re: [PATCH] Fix for ldbl-128ibm acosl/asinl inaccuracies

On 05/02/2012 04:32 PM, Adhemerval Zanella wrote: 
Commit 'adfbc8ac9e192b6e3007f7a47852df937afa2145' (Fix x86 acos near 1) and
5ba3cc691c856e5c67a7d4cd4713f20a79f7ba81 (Fix acos (-1) in round-downwards mode on x86)
added some additional asinl/acosl tests that triggered some issues on ldbl-128ibm
implementation.

The issues are due current implementation does a wrong fabsl and long double comparisons
and values like 0x0.ffffffffffffffffp0L are treated it were higher than |x|. The patch
corrects it by using __fabsl correctly and by doing long double comparisons, which
uses both high and low double instead of just the high part in current implementation.

I also added the ULP update for above commits.

---

2012-05-02 Adhemerval Zanella<azanella@linux.vnet.ibm.com>

	* sysdeps/ieee754/ldbl-128ibm/e_acosl.c: Fix long double comparison
	inaccuracies.


Please always add the function name, e.g.:
 	* sysdeps/ieee754/ldbl-128ibm/e_acosl.c (__ieee754_acosl):

I'll make those changes and commit now,

thanks,
Andreas

	* sysdeps/ieee754/ldbl-128ibm/e_asinl.c: Likewise.
	* sysdeps/powerpc/fpu/libm-test-ulps: Update.

diff --git a/sysdeps/ieee754/ldbl-128ibm/e_acosl.c b/sysdeps/ieee754/ldbl-128ibm/e_acosl.c
index 533b597..5d2af30 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 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;
diff --git a/sysdeps/powerpc/fpu/libm-test-ulps b/sysdeps/powerpc/fpu/libm-test-ulps
index 43fa642..5abff41 100644
--- a/sysdeps/powerpc/fpu/libm-test-ulps
+++ b/sysdeps/powerpc/fpu/libm-test-ulps
@@ -4,12 +4,113 @@
  Test "acos (2e-17) == 1.57079632679489659923132169163975144":
  ildouble: 1
  ldouble: 1
+Test "acos (-0x0.ffffffff8p0) == 3.1415773948007305904329067627145550395696":
+ldouble: 1
+ildouble: 1
+Test "acos (-0x0.ffffffp0) == 3.1412473866050770348750401337968641476999":
+ldouble: 1
+ildouble: 1
+
+# acos_downward
+Test "acos_downward (-0.5) == M_PI_6l*4.0":
+double: 1
+idouble: 1
+ldouble: 1
+ildouble: 1
+Test "acos_downward (0.5) == M_PI_6l*2.0":
+float: 1
+ifloat: 1
+double: 1
+idouble: 1
+ldouble: 1
+ildouble: 1
+
+# acos_towardzero
+Test "acos_towardzero (-0.5) == M_PI_6l*4.0":
+double: 1
+idouble: 1
+ldouble: 1
+ildouble: 1
+Test "acos_towardzero (0.5) == M_PI_6l*2.0":
+float: 1
+ifloat: 1
+double: 1
+idouble: 1
+ldouble: 1
+ildouble: 1
+
+# acos_upward
+Test "acos_upward (-0) == pi/2":
+ldouble: 2
+ildouble: 2
+Test "acos_upward (-1) == pi":
+ldouble: 2
+ildouble: 2
+Test "acos_upward (0) == pi/2":
+ldouble: 2
+ildouble: 2

  # asin
+Test "asin (-0x0.ffffffff8p0) == -1.5707810680058339712015850710748035974710":
+ldouble: 1
+ildouble: 1
+Test "asin (-0x0.ffffffp0) == -1.5704510598101804156437184421571127056013":
+ldouble: 1
+ildouble: 1
+Test "asin (0x0.ffffffff8p0) == 1.5707810680058339712015850710748035974710":
+ldouble: 1
+ildouble: 1
+Test "asin (0x0.ffffffp0) == 1.5704510598101804156437184421571127056013":
+ldouble: 1
+ildouble: 1
  Test "asin (0.75) == 0.848062078981481008052944338998418080":
  ildouble: 2
  ldouble: 2

+# asin_downward
+Test "asin_downward (-0.5) == -pi/6":
+double: 1
+idouble: 1
+ldouble: 1
+ildouble: 1
+Test "asin_downward (0.5) == pi/6":
+double: 1
+idouble: 1
+ldouble: 1
+ildouble: 1
+Test "asin_downward (-1.0) == -pi/2":
+ldouble: 1
+ildouble: 1
+Test "asin_downward (1.0) == pi/2":
+float: 1
+ifloat: 1
+
+# asin_towardzero
+Test "asin_towardzero (-0.5) == -pi/6":
+double: 1
+idouble: 1
+ldouble: 1
+ildouble: 1
+Test "asin_towardzero (0.5) == pi/6":
+double: 1
+idouble: 1
+ldouble: 1
+ildouble: 1
+Test "asin_towardzero (-1.0) == -pi/2":
+float: 1
+ifloat:1
+Test "asin_towardzero (1.0) == pi/2":
+float: 1
+ifloat: 1
+
+# asin_upward
+Test "asin_upward (-1.0) == -pi/2":
+float: 1
+ifloat: 1
+Test "asin_upward (1.0) == pi/2":
+ldouble: 1
+ildouble: 1
+
  # atan2
  Test "atan2 (-0.00756827042671106339, -.001792735857538728036) == -1.80338464113663849327153994379639112":
  ildouble: 1
@@ -2061,6 +2162,30 @@ Function: "acos":
  ildouble: 1
  ldouble: 1

+Function: "acos_downward":
+float: 1
+ifloat: 1
+double: 1
+idouble: 1
+ldouble: 1
+ildouble: 1
+
+Function: "acos_tonearest":
+ldouble: 1
+ildouble: 1
+
+Function: "acos_towardzero":
+float: 1
+ifloat: 1
+double: 1
+idouble: 1
+ldouble: 1
+ildouble: 1
+
+Function: "acos_upward":
+ldouble: 2
+ildouble: 2
+
  Function: "acosh":
  ildouble: 1
  ldouble: 1
@@ -2069,6 +2194,32 @@ Function: "asin":
  ildouble: 2
  ldouble: 2

+Function: "asin_downward":
+float: 1
+ifloat: 1
+double: 1
+idouble: 1
+ldouble: 1
+ildouble: 1
+
+Function: "asin_tonearest":
+ldouble: 1
+ildouble: 1
+
+Function: "asin_towardzero":
+float: 1
+ifloat: 1
+double: 1
+idouble: 1
+ldouble: 1
+ildouble: 1
+
+Function: "asin_upward":
+float: 1
+ifloat: 1
+ldouble: 1
+ildouble: 1
+
  Function: "asinh":
  ildouble: 1
  ldouble: 1



--
 Andreas Jaeger aj@{suse.com,opensuse.org} Twitter/Identica: jaegerandi
  SUSE LINUX Products GmbH, Maxfeldstr. 5, 90409 Nürnberg, Germany
   GF: Jeff Hawn,Jennifer Guild,Felix Imendörffer,HRB16746 (AG Nürnberg)
    GPG fingerprint = 93A3 365E CE47 B889 DF7F  FED1 389A 563C C272 A126

Index Nav: [Date Index] [Subject Index] [Author Index] [Thread Index]
Message Nav: [Date Prev] [Date Next] [Thread Prev] [Thread Next]