diff mbox series

[1/5] Remove slow paths from asin

Message ID VE1PR08MB5599F6B916F2EE343AA29F5683D20@VE1PR08MB5599.eurprd08.prod.outlook.com
State Superseded
Headers
Series Remove remaining slow paths from math functions |

Commit Message

Wilco Dijkstra Jan. 4, 2021, 12:15 p.m. UTC 
  Remove slow paths from asin - there is quite a lot of redundant slow code that can
be removed while keeping ULP below 1. Add ULP annotations. Update AArch64
libm-test-ulps for asin.

Passes GLIBC testsuite.

---

Comments

Paul Zimmermann Jan. 4, 2021, 4:25 p.m. UTC | #1 
Dear Wilco,

I have a few comments on this first patch:

- the title only mentions asin, but the patch also deals with acos

- it would be nice to explain where the "Max ULP" come from, to allow
  future developers to debug the code if needed.

- it would be nice to have a summary with the maximal ulp error at the
  top of each function (__ieee754_asin and __ieee754_acos)

- for __ieee754_acos the largest claimed error is 0.52 ULP for
  0.5 <= |x| < 0.96875, but I find an error of more than 0.522 ulps
  for x = 0x1.dffff58d80953p-1 (on x86_64).

Best regards,
Paul
Szabolcs Nagy Jan. 4, 2021, 5:06 p.m. UTC | #2 
The 01/04/2021 17:25, Paul Zimmermann wrote:
>        Dear Wilco,
> 
> I have a few comments on this first patch:
> 
> - the title only mentions asin, but the patch also deals with acos
> 
> - it would be nice to explain where the "Max ULP" come from, to allow
>   future developers to debug the code if needed.
> 
> - it would be nice to have a summary with the maximal ulp error at the
>   top of each function (__ieee754_asin and __ieee754_acos)
> 
> - for __ieee754_acos the largest claimed error is 0.52 ULP for
>   0.5 <= |x| < 0.96875, but I find an error of more than 0.522 ulps
>   for x = 0x1.dffff58d80953p-1 (on x86_64).

the Max ULP annotations are based on the assumption that
the bounds checks in the removed code paths caught all the
incorrectly rounded half-way cases.

previously it was noted that there are incorrectly rounded
cases in asin that were not caught by those checks:
https://sourceware.org/legacy-ml/libc-alpha/2020-01/msg00574.html

so either the Max ULP annotation was computed incorrectly
or the same errors should happen in the original acos too.

in my experience the original code is fairly conservative
so in most cases i'd expect the Max ULP to be an upper
bound with significant margin, not the actual worst case
error.

i think it's worth keeping the annotation even if the
original code was not entirely correct.
Paul Zimmermann Jan. 4, 2021, 5:41 p.m. UTC | #3 
Dear Szabolcs,

> previously it was noted that there are incorrectly rounded
> cases in asin that were not caught by those checks:
> https://sourceware.org/legacy-ml/libc-alpha/2020-01/msg00574.html

the patch you mention (commit f67f9c9) was already entitled
"ieee754: Remove slow paths from asin and acos". This was confusing.
I believe it did only remove the "very slow" paths.

> in my experience the original code is fairly conservative
> so in most cases i'd expect the Max ULP to be an upper
> bound with significant margin, not the actual worst case
> error.

here are the claimed bounds and the maximal errors I found so far after
applying Wilco's patches (rounded up to 3 digits):

      claimed found so far
acos  0.52    0.523
asin  0.523   0.516
atan  0.56    0.523
tan   0.62    0.612
atan2 0.56    0.524

Paul
Wilco Dijkstra Jan. 7, 2021, 7:36 p.m. UTC | #4 
Hi Paul, 

Thanks for the reviews, I have updated the patches accordingly.

> the patch you mention (commit f67f9c9) was already entitled
> "ieee754: Remove slow paths from asin and acos". This was confusing.
> I believe it did only remove the "very slow" paths.

I've updated the title in v2.

> here are the claimed bounds and the maximal errors I found so far after
> applying Wilco's patches (rounded up to 3 digits):
>
>      claimed found so far
> acos  0.52    0.523
> asin  0.523   0.516
> atan  0.56    0.523
> tan   0.62    0.612
> atan2 0.56    0.524

Thanks for that, I have added these bounds at the start of each function.

> I have only one further remark:
>
> -typedef union { int4 i[2]; double x; } mynumber;
> +typedef union { int4 i[2]; double x; double d; } mynumber;
>
> Why did you add another 'double' field d? Only the 'd' field is used in
> e_atan2.c, why not replace it by 'x'? Maybe to avoid a too large diff?

Yes, the purpose is just to remove the slow paths and minimise the
changes, so this was the quickest and simplest way to remove the uses
of 'number'.

Cheers,
Wilco
Paul Zimmermann Jan. 8, 2021, 10:29 a.m. UTC | #5 
Dear Wilco,

the v2 of this series is ok for me, except in the meantime I have found a
larger error for tan:

tan 0 -1 0x1.c673a473b3503p+3 [0.619] 0.618365 0.6183652837454632

thus you could replace ~0.612 by ~0.619.

Let me summarize the largest errors we get with 2.32 and after this patch
(on random sampling):

            2.32      after patch (with corresponding input)
asin       0.500      0.516 0x1.07fffeadab7ecp-3
acos       0.501      0.523 0x1.dffff776c7505p-1
tan        0.500      0.619 0x1.c673a473b3503p+3
atan       0.500      0.523 -0x1.f9004176e07bep-4
atan2      0.500      0.524 0x1.4aa3aada7c44fp-996,0x1.4f32ff56386e8p-993

Best regards,
Paul
diff mbox series

Patch

diff --git a/sysdeps/aarch64/libm-test-ulps b/sysdeps/aarch64/libm-test-ulps
index 22fcf8db73dc444c25e0c356b1e0036571edd112..bbadf667ee4b7a0cf80506d321553f064049c516 100644
--- a/sysdeps/aarch64/libm-test-ulps
+++ b/sysdeps/aarch64/libm-test-ulps
@@ -41,6 +41,7 @@  float: 2
 ldouble: 2
 
 Function: "asin":
+double: 1
 float: 1
 ldouble: 1
 
@@ -55,7 +56,7 @@  float: 1
 ldouble: 1
 
 Function: "asin_upward":
-double: 1
+double: 2
 float: 1
 ldouble: 2
 
diff --git a/sysdeps/ieee754/dbl-64/e_asin.c b/sysdeps/ieee754/dbl-64/e_asin.c
index 8a3b26f6645b66818ad0a57fe833355a3e9961e6..c01e8a34517e4a33f126bbab5c132379b4b58a4d 100644
--- a/sysdeps/ieee754/dbl-64/e_asin.c
+++ b/sysdeps/ieee754/dbl-64/e_asin.c
@@ -21,8 +21,7 @@ 
 /*                                                                */
 /*     FUNCTIONS: uasin                                           */
 /*                uacos                                           */
-/* FILES NEEDED: dla.h endian.h mpa.h mydefs.h  usncs.h           */
-/*               doasin.c sincos32.c dosincos.c mpa.c             */
+/* FILES NEEDED: dla.h endian.h mydefs.h  usncs.h                 */
 /*               sincos.tbl  asincos.tbl  powtwo.tbl root.tbl     */
 /*                                                                */
 /******************************************************************/
@@ -31,7 +30,6 @@ 
 #include "asincos.tbl"
 #include "root.tbl"
 #include "powtwo.tbl"
-#include "MathLib.h"
 #include "uasncs.h"
 #include <float.h>
 #include <math.h>
@@ -43,15 +41,10 @@ 
 # define SECTION
 #endif
 
-void __doasin(double x, double dx, double w[]);
-void __dubsin(double x, double dx, double v[]);
-void __dubcos(double x, double dx, double v[]);
-void __docos(double x, double dx, double v[]);
-
 double
 SECTION
 __ieee754_asin(double x){
-  double x1,x2,xx,s1,s2,res1,p,t,res,r,cor,cc,y,c,z,w[2];
+  double x2,xx,res1,p,t,res,r,cor,cc,y,c,z;
   mynumber u,v;
   int4 k,m,n;
 
@@ -70,27 +63,8 @@  __ieee754_asin(double x){
     x2 = x*x;
     t = (((((f6*x2 + f5)*x2 + f4)*x2 + f3)*x2 + f2)*x2 + f1)*(x2*x);
     res = x+t;         /*  res=arcsin(x) according to Taylor series  */
-    cor = (x-res)+t;
-    if (res == res+1.025*cor) return res;
-    else {
-      x1 = x+big;
-      xx = x*x;
-      x1 -= big;
-      x2 = x - x1;
-      p = x1*x1*x1;
-      s1 = a1.x*p;
-      s2 = ((((((c7*xx + c6)*xx + c5)*xx + c4)*xx + c3)*xx + c2)*xx*xx*x +
-	     ((a1.x+a2.x)*x2*x2+ 0.5*x1*x)*x2) + a2.x*p;
-      res1 = x+s1;
-      s2 = ((x-res1)+s1)+s2;
-      res = res1+s2;
-      cor = (res1-res)+s2;
-      if (res == res+1.00014*cor) return res;
-      else {
-	__doasin(x,0,w);
-	return w[0];
-      }
-    }
+    /* Max ULP is 0.512.  */
+    return res;
   }
   /*---------------------0.125 <= |x| < 0.5 -----------------------------*/
   else if (k < 0x3fe00000) {
@@ -103,26 +77,8 @@  __ieee754_asin(double x){
      +xx*asncs.x[n+6]))))+asncs.x[n+7];
     t+=p;
     res =asncs.x[n+8] +t;
-    cor = (asncs.x[n+8]-res)+t;
-    if (res == res+1.05*cor) return (m>0)?res:-res;
-    else {
-      r=asncs.x[n+8]+xx*asncs.x[n+9];
-      t=((asncs.x[n+8]-r)+xx*asncs.x[n+9])+(p+xx*asncs.x[n+10]);
-      res = r+t;
-      cor = (r-res)+t;
-      if (res == res+1.0005*cor) return (m>0)?res:-res;
-      else {
-	res1=res+1.1*cor;
-	z=0.5*(res1-res);
-	__dubsin(res,z,w);
-	z=(w[0]-fabs(x))+w[1];
-	if (z>1.0e-27) return (m>0)?min(res,res1):-min(res,res1);
-	else if (z<-1.0e-27) return (m>0)?max(res,res1):-max(res,res1);
-	else {
-	  return (m>0)?res:-res;
-	}
-      }
-    }
+    /* Max ULP is 0.523.  */
+    return (m>0)?res:-res;
   }    /*   else  if (k < 0x3fe00000)    */
   /*-------------------- 0.5 <= |x| < 0.75 -----------------------------*/
   else
@@ -135,26 +91,8 @@  __ieee754_asin(double x){
 	   +xx*(asncs.x[n+6]+xx*asncs.x[n+7])))))+asncs.x[n+8];
     t+=p;
     res =asncs.x[n+9] +t;
-    cor = (asncs.x[n+9]-res)+t;
-    if (res == res+1.01*cor) return (m>0)?res:-res;
-    else {
-      r=asncs.x[n+9]+xx*asncs.x[n+10];
-      t=((asncs.x[n+9]-r)+xx*asncs.x[n+10])+(p+xx*asncs.x[n+11]);
-      res = r+t;
-      cor = (r-res)+t;
-      if (res == res+1.0005*cor) return (m>0)?res:-res;
-      else {
-	res1=res+1.1*cor;
-	z=0.5*(res1-res);
-	__dubsin(res,z,w);
-	z=(w[0]-fabs(x))+w[1];
-	if (z>1.0e-27) return (m>0)?min(res,res1):-min(res,res1);
-	else if (z<-1.0e-27) return (m>0)?max(res,res1):-max(res,res1);
-	else {
-	  return (m>0)?res:-res;
-	}
-      }
-    }
+    /* Max ULP is 0.505.  */
+    return (m>0)?res:-res;
   }    /*   else  if (k < 0x3fe80000)    */
   /*--------------------- 0.75 <= |x|< 0.921875 ----------------------*/
   else
@@ -167,28 +105,8 @@  __ieee754_asin(double x){
      +xx*(asncs.x[n+6]+xx*(asncs.x[n+7]+xx*asncs.x[n+8]))))))+asncs.x[n+9];
     t+=p;
     res =asncs.x[n+10] +t;
-    cor = (asncs.x[n+10]-res)+t;
-    if (res == res+1.01*cor) return (m>0)?res:-res;
-    else {
-      r=asncs.x[n+10]+xx*asncs.x[n+11];
-      t=((asncs.x[n+10]-r)+xx*asncs.x[n+11])+(p+xx*asncs.x[n+12]);
-      res = r+t;
-      cor = (r-res)+t;
-      if (res == res+1.0008*cor) return (m>0)?res:-res;
-      else {
-	res1=res+1.1*cor;
-	z=0.5*(res1-res);
-	y=hp0.x-res;
-	z=((hp0.x-y)-res)+(hp1.x-z);
-	__dubcos(y,z,w);
-	z=(w[0]-fabs(x))+w[1];
-	if (z>1.0e-27) return (m>0)?min(res,res1):-min(res,res1);
-	else if (z<-1.0e-27) return (m>0)?max(res,res1):-max(res,res1);
-	else {
-	  return (m>0)?res:-res;
-	}
-      }
-    }
+    /* Max ULP is 0.505.  */
+    return (m>0)?res:-res;
   }    /*   else  if (k < 0x3fed8000)    */
   /*-------------------0.921875 <= |x| < 0.953125 ------------------------*/
   else
@@ -203,29 +121,8 @@  __ieee754_asin(double x){
 		      xx*asncs.x[n+9])))))))+asncs.x[n+10];
     t+=p;
     res =asncs.x[n+11] +t;
-    cor = (asncs.x[n+11]-res)+t;
-    if (res == res+1.01*cor) return (m>0)?res:-res;
-    else {
-      r=asncs.x[n+11]+xx*asncs.x[n+12];
-      t=((asncs.x[n+11]-r)+xx*asncs.x[n+12])+(p+xx*asncs.x[n+13]);
-      res = r+t;
-      cor = (r-res)+t;
-      if (res == res+1.0007*cor) return (m>0)?res:-res;
-      else {
-	res1=res+1.1*cor;
-	z=0.5*(res1-res);
-	y=(hp0.x-res)-z;
-	z=y+hp1.x;
-	y=(y-z)+hp1.x;
-	__dubcos(z,y,w);
-	z=(w[0]-fabs(x))+w[1];
-	if (z>1.0e-27) return (m>0)?min(res,res1):-min(res,res1);
-	else if (z<-1.0e-27) return (m>0)?max(res,res1):-max(res,res1);
-	else {
-	  return (m>0)?res:-res;
-	}
-      }
-    }
+    /* Max ULP is 0.505.  */
+    return (m>0)?res:-res;
   }    /*   else  if (k < 0x3fee8000)    */
 
   /*--------------------0.953125 <= |x| < 0.96875 ------------------------*/
@@ -241,29 +138,8 @@  __ieee754_asin(double x){
 		    xx*(asncs.x[n+9]+xx*asncs.x[n+10]))))))))+asncs.x[n+11];
     t+=p;
     res =asncs.x[n+12] +t;
-    cor = (asncs.x[n+12]-res)+t;
-    if (res == res+1.01*cor) return (m>0)?res:-res;
-    else {
-      r=asncs.x[n+12]+xx*asncs.x[n+13];
-      t=((asncs.x[n+12]-r)+xx*asncs.x[n+13])+(p+xx*asncs.x[n+14]);
-      res = r+t;
-      cor = (r-res)+t;
-      if (res == res+1.0007*cor) return (m>0)?res:-res;
-      else {
-	res1=res+1.1*cor;
-	z=0.5*(res1-res);
-	y=(hp0.x-res)-z;
-	z=y+hp1.x;
-	y=(y-z)+hp1.x;
-	__dubcos(z,y,w);
-	z=(w[0]-fabs(x))+w[1];
-	if (z>1.0e-27) return (m>0)?min(res,res1):-min(res,res1);
-	else if (z<-1.0e-27) return (m>0)?max(res,res1):-max(res,res1);
-	else {
-	  return (m>0)?res:-res;
-	}
-      }
-    }
+    /* Max ULP is 0.505.  */
+    return (m>0)?res:-res;
   }    /*   else  if (k < 0x3fef0000)    */
   /*--------------------0.96875 <= |x| < 1 --------------------------------*/
   else
@@ -282,16 +158,8 @@  __ieee754_asin(double x){
     cor = (hp1.x - 2.0*cc)-2.0*(y+cc)*p;
     res1 = hp0.x - 2.0*y;
     res =res1 + cor;
-    if (res == res+1.003*((res1-res)+cor)) return (m>0)?res:-res;
-    else {
-      c=y+cc;
-      cc=(y-c)+cc;
-      __doasin(c,cc,w);
-      res1=hp0.x-2.0*w[0];
-      cor=((hp0.x-res1)-2.0*w[0])+(hp1.x-2.0*w[1]);
-      res = res1+cor;
-      return (m>0)?res:-res;
-    }
+    /* Max ULP is 0.5015.  */
+    return (m>0)?res:-res;
   }    /*   else  if (k < 0x3ff00000)    */
   /*---------------------------- |x|>=1 -------------------------------*/
   else if (k==0x3ff00000 && u.i[LOW_HALF]==0) return (m>0)?hp0.x:-hp0.x;
@@ -319,7 +187,7 @@  double
 SECTION
 __ieee754_acos(double x)
 {
-  double x1,x2,xx,s1,s2,res1,p,t,res,r,cor,cc,y,c,z,w[2],eps;
+  double x2,xx,res1,p,t,res,r,cor,cc,y,c,z;
   mynumber u,v;
   int4 k,m,n;
   u.x = x;
@@ -336,32 +204,8 @@  __ieee754_acos(double x)
     r=hp0.x-x;
     cor=(((hp0.x-r)-x)+hp1.x)-t;
     res = r+cor;
-    cor = (r-res)+cor;
-    if (res == res+1.004*cor) return res;
-    else {
-      x1 = x+big;
-      xx = x*x;
-      x1 -= big;
-      x2 = x - x1;
-      p = x1*x1*x1;
-      s1 = a1.x*p;
-      s2 = ((((((c7*xx + c6)*xx + c5)*xx + c4)*xx + c3)*xx + c2)*xx*xx*x +
-	    ((a1.x+a2.x)*x2*x2+ 0.5*x1*x)*x2) + a2.x*p;
-      res1 = x+s1;
-      s2 = ((x-res1)+s1)+s2;
-      r=hp0.x-res1;
-      cor=(((hp0.x-r)-res1)+hp1.x)-s2;
-      res = r+cor;
-      cor = (r-res)+cor;
-      if (res == res+1.00004*cor) return res;
-      else {
-	__doasin(x,0,w);
-	r=hp0.x-w[0];
-	cor=((hp0.x-r)-w[0])+(hp1.x-w[1]);
-	res=r+cor;
-	return res;
-      }
-    }
+    /* Max ULP is 0.502.  */
+    return res;
   }    /*   else  if (k < 0x3fc00000)    */
   /*----------------------  0.125 <= |x| < 0.5 --------------------*/
   else
@@ -377,35 +221,16 @@  __ieee754_acos(double x)
     y = (m>0)?(hp0.x-asncs.x[n+8]):(hp0.x+asncs.x[n+8]);
     t = (m>0)?(hp1.x-t):(hp1.x+t);
     res = y+t;
-    if (res == res+1.02*((y-res)+t)) return res;
-    else {
-      r=asncs.x[n+8]+xx*asncs.x[n+9];
-      t=((asncs.x[n+8]-r)+xx*asncs.x[n+9])+(p+xx*asncs.x[n+10]);
-      if (m>0)
-	{p = hp0.x-r; t = (((hp0.x-p)-r)-t)+hp1.x; }
-      else
-	{p = hp0.x+r; t = ((hp0.x-p)+r)+(hp1.x+t); }
-      res = p+t;
-      cor = (p-res)+t;
-      if (res == (res+1.0002*cor)) return res;
-      else {
-	res1=res+1.1*cor;
-	z=0.5*(res1-res);
-	__docos(res,z,w);
-	z=(w[0]-x)+w[1];
-	if (z>1.0e-27) return max(res,res1);
-	else if (z<-1.0e-27) return min(res,res1);
-	else return res;
-      }
-    }
+   /* Max ULP is 0.51.  */
+    return res;
   }    /*   else  if (k < 0x3fe00000)    */
 
   /*--------------------------- 0.5 <= |x| < 0.75 ---------------------*/
   else
   if (k < 0x3fe80000) {
     n = 1056+((k&0x000fe000)>>11)*3;
-    if (m>0) {xx = x - asncs.x[n]; eps=1.04; }
-    else {xx = -x - asncs.x[n]; eps=1.02; }
+    if (m>0) {xx = x - asncs.x[n]; }
+    else {xx = -x - asncs.x[n]; }
     t = asncs.x[n+1]*xx;
     p=xx*xx*(asncs.x[n+2]+xx*(asncs.x[n+3]+xx*(asncs.x[n+4]+
 		   xx*(asncs.x[n+5]+xx*(asncs.x[n+6]+
@@ -414,33 +239,16 @@  __ieee754_acos(double x)
    y = (m>0)?(hp0.x-asncs.x[n+9]):(hp0.x+asncs.x[n+9]);
    t = (m>0)?(hp1.x-t):(hp1.x+t);
    res = y+t;
-   if (res == res+eps*((y-res)+t)) return res;
-   else {
-     r=asncs.x[n+9]+xx*asncs.x[n+10];
-     t=((asncs.x[n+9]-r)+xx*asncs.x[n+10])+(p+xx*asncs.x[n+11]);
-     if (m>0) {p = hp0.x-r; t = (((hp0.x-p)-r)-t)+hp1.x; eps=1.0004; }
-     else   {p = hp0.x+r; t = ((hp0.x-p)+r)+(hp1.x+t); eps=1.0002; }
-     res = p+t;
-     cor = (p-res)+t;
-     if (res == (res+eps*cor)) return res;
-     else {
-       res1=res+1.1*cor;
-       z=0.5*(res1-res);
-       __docos(res,z,w);
-       z=(w[0]-x)+w[1];
-       if (z>1.0e-27) return max(res,res1);
-       else if (z<-1.0e-27) return min(res,res1);
-       else return res;
-     }
-   }
+   /* Max ULP is 0.52.  */
+   return res;
   }    /*   else  if (k < 0x3fe80000)    */
 
 /*------------------------- 0.75 <= |x| < 0.921875 -------------*/
   else
   if (k < 0x3fed8000) {
     n = 992+((k&0x000fe000)>>13)*13;
-    if (m>0) {xx = x - asncs.x[n]; eps = 1.04; }
-    else {xx = -x - asncs.x[n]; eps = 1.01; }
+    if (m>0) {xx = x - asncs.x[n]; }
+    else {xx = -x - asncs.x[n]; }
     t = asncs.x[n+1]*xx;
     p=xx*xx*(asncs.x[n+2]+xx*(asncs.x[n+3]+xx*(asncs.x[n+4]+
 		      xx*(asncs.x[n+5]+xx*(asncs.x[n+6]+xx*(asncs.x[n+7]+
@@ -449,33 +257,16 @@  __ieee754_acos(double x)
     y = (m>0)?(hp0.x-asncs.x[n+10]):(hp0.x+asncs.x[n+10]);
     t = (m>0)?(hp1.x-t):(hp1.x+t);
     res = y+t;
-    if (res == res+eps*((y-res)+t)) return res;
-    else {
-      r=asncs.x[n+10]+xx*asncs.x[n+11];
-      t=((asncs.x[n+10]-r)+xx*asncs.x[n+11])+(p+xx*asncs.x[n+12]);
-      if (m>0) {p = hp0.x-r; t = (((hp0.x-p)-r)-t)+hp1.x; eps=1.0032; }
-      else   {p = hp0.x+r; t = ((hp0.x-p)+r)+(hp1.x+t); eps=1.0008; }
-      res = p+t;
-      cor = (p-res)+t;
-      if (res == (res+eps*cor)) return res;
-      else {
-	res1=res+1.1*cor;
-	z=0.5*(res1-res);
-	__docos(res,z,w);
-	z=(w[0]-x)+w[1];
-	if (z>1.0e-27) return max(res,res1);
-	else if (z<-1.0e-27) return min(res,res1);
-	else return res;
-      }
-    }
+   /* Max ULP is 0.52.  */
+    return res;
   }    /*   else  if (k < 0x3fed8000)    */
 
 /*-------------------0.921875 <= |x| < 0.953125 ------------------*/
   else
   if (k < 0x3fee8000) {
     n = 884+((k&0x000fe000)>>13)*14;
-    if (m>0) {xx = x - asncs.x[n]; eps=1.04; }
-    else {xx = -x - asncs.x[n]; eps =1.005; }
+    if (m>0) {xx = x - asncs.x[n]; }
+    else {xx = -x - asncs.x[n]; }
     t = asncs.x[n+1]*xx;
     p=xx*xx*(asncs.x[n+2]+xx*(asncs.x[n+3]+xx*(asncs.x[n+4]+
 		   xx*(asncs.x[n+5]+xx*(asncs.x[n+6]
@@ -485,33 +276,16 @@  __ieee754_acos(double x)
     y = (m>0)?(hp0.x-asncs.x[n+11]):(hp0.x+asncs.x[n+11]);
     t = (m>0)?(hp1.x-t):(hp1.x+t);
     res = y+t;
-    if (res == res+eps*((y-res)+t)) return res;
-    else {
-      r=asncs.x[n+11]+xx*asncs.x[n+12];
-      t=((asncs.x[n+11]-r)+xx*asncs.x[n+12])+(p+xx*asncs.x[n+13]);
-      if (m>0) {p = hp0.x-r; t = (((hp0.x-p)-r)-t)+hp1.x; eps=1.0030; }
-      else   {p = hp0.x+r; t = ((hp0.x-p)+r)+(hp1.x+t); eps=1.0005; }
-      res = p+t;
-      cor = (p-res)+t;
-      if (res == (res+eps*cor)) return res;
-      else {
-	res1=res+1.1*cor;
-	z=0.5*(res1-res);
-	__docos(res,z,w);
-	z=(w[0]-x)+w[1];
-	if (z>1.0e-27) return max(res,res1);
-	else if (z<-1.0e-27) return min(res,res1);
-	else return res;
-      }
-    }
+   /* Max ULP is 0.52.  */
+    return res;
   }    /*   else  if (k < 0x3fee8000)    */
 
   /*--------------------0.953125 <= |x| < 0.96875 ----------------*/
   else
   if (k < 0x3fef0000) {
     n = 768+((k&0x000fe000)>>13)*15;
-    if (m>0) {xx = x - asncs.x[n]; eps=1.04; }
-    else {xx = -x - asncs.x[n]; eps=1.005;}
+    if (m>0) {xx = x - asncs.x[n]; }
+    else {xx = -x - asncs.x[n]; }
     t = asncs.x[n+1]*xx;
     p=xx*xx*(asncs.x[n+2]+xx*(asncs.x[n+3]+xx*(asncs.x[n+4]+
 	    xx*(asncs.x[n+5]+xx*(asncs.x[n+6]
@@ -521,25 +295,8 @@  __ieee754_acos(double x)
     y = (m>0)?(hp0.x-asncs.x[n+12]):(hp0.x+asncs.x[n+12]);
    t = (m>0)?(hp1.x-t):(hp1.x+t);
    res = y+t;
-   if (res == res+eps*((y-res)+t)) return res;
-   else {
-     r=asncs.x[n+12]+xx*asncs.x[n+13];
-     t=((asncs.x[n+12]-r)+xx*asncs.x[n+13])+(p+xx*asncs.x[n+14]);
-     if (m>0) {p = hp0.x-r; t = (((hp0.x-p)-r)-t)+hp1.x; eps=1.0030; }
-     else   {p = hp0.x+r; t = ((hp0.x-p)+r)+(hp1.x+t); eps=1.0005; }
-     res = p+t;
-     cor = (p-res)+t;
-     if (res == (res+eps*cor)) return res;
-     else {
-       res1=res+1.1*cor;
-       z=0.5*(res1-res);
-       __docos(res,z,w);
-       z=(w[0]-x)+w[1];
-       if (z>1.0e-27) return max(res,res1);
-       else if (z<-1.0e-27) return min(res,res1);
-       else return res;
-     }
-   }
+   /* Max ULP is 0.52.  */
+   return res;
   }    /*   else  if (k < 0x3fef0000)    */
   /*-----------------0.96875 <= |x| < 1 ---------------------------*/
 
@@ -560,28 +317,14 @@  __ieee754_acos(double x)
       cor = (hp1.x - cc)-(y+cc)*p;
       res1 = hp0.x - y;
       res =res1 + cor;
-      if (res == res+1.002*((res1-res)+cor)) return (res+res);
-      else {
-	c=y+cc;
-	cc=(y-c)+cc;
-	__doasin(c,cc,w);
-	res1=hp0.x-w[0];
-	cor=((hp0.x-res1)-w[0])+(hp1.x-w[1]);
-	res = res1+cor;
-	return (res+res);
-      }
+      /* Max ULP is 0.501.  */
+      return (res+res);
     }
     else {
       cor = cc+p*(y+cc);
       res = y + cor;
-      if (res == res+1.03*((y-res)+cor)) return (res+res);
-      else {
-	c=y+cc;
-	cc=(y-c)+cc;
-	__doasin(c,cc,w);
-	res = w[0];
-	return (res+res);
-      }
+      /* Max ULP is 0.515.  */
+      return (res+res);
     }
   }    /*   else  if (k < 0x3ff00000)    */