Fix sin, cos, tan in non-default rounding modes (bug 3976).

1 files changed, 68 insertions, 49 deletions

diff --git a/sysdeps/ieee754/dbl-64/s_tan.c b/sysdeps/ieee754/dbl-64/s_tan.c
index 962a4eba6b..2c26756ff2 100644
--- a/sysdeps/ieee754/dbl-64/s_tan.c
+++ b/sysdeps/ieee754/dbl-64/s_tan.c
@@ -39,6 +39,8 @@
 #include "mpa.h"
 #include "MathLib.h"
 #include "math.h"
+#include "math_private.h"
+#include <fenv.h>
 
 #ifndef SECTION
 # define SECTION
@@ -66,21 +68,27 @@ tan(double x) {
   mp_no mpy;
 #endif
 
+  fenv_t env;
+  double retval;
+
   int __branred(double, double *, double *);
   int __mpranred(double, mp_no *, int);
 
+  libc_feholdexcept_setround (&env, FE_TONEAREST);
+
   /* x=+-INF, x=NaN */
   num.d = x;  ux = num.i[HIGH_HALF];
   if ((ux&0x7ff00000)==0x7ff00000) {
     if ((ux&0x7fffffff)==0x7ff00000)
       __set_errno (EDOM);
-    return x-x;
+    retval = x-x;
+    goto ret;
   }
 
   w=(x<ZERO) ? -x : x;
 
   /* (I) The case abs(x) <= 1.259e-8 */
-  if (w<=g1.d)  return x;
+  if (w<=g1.d) { retval = x; goto ret; }
 
   /* (II) The case 1.259e-8 < abs(x) <= 0.0608 */
   if (w<=g2.d) {
@@ -88,7 +96,7 @@ tan(double x) {
     /* First stage */
     x2 = x*x;
     t2 = x*x2*(d3.d+x2*(d5.d+x2*(d7.d+x2*(d9.d+x2*d11.d))));
-    if ((y=x+(t2-u1.d*t2)) == x+(t2+u1.d*t2))  return y;
+    if ((y=x+(t2-u1.d*t2)) == x+(t2+u1.d*t2)) { retval = y; goto ret; }
 
     /* Second stage */
     c1 = x2*(a15.d+x2*(a17.d+x2*(a19.d+x2*(a21.d+x2*(a23.d+x2*(a25.d+
@@ -108,8 +116,9 @@ tan(double x) {
     MUL2(x2,xx2,c2,cc2,c1,cc1,t1,t2,t3,t4,t5,t6,t7,t8)
     MUL2(x ,zero.d,c1,cc1,c2,cc2,t1,t2,t3,t4,t5,t6,t7,t8)
     ADD2(x    ,zero.d,c2,cc2,c1,cc1,t1,t2)
-    if ((y=c1+(cc1-u2.d*c1)) == c1+(cc1+u2.d*c1))  return y;
-    return tanMp(x);
+    if ((y=c1+(cc1-u2.d*c1)) == c1+(cc1+u2.d*c1)) { retval = y; goto ret; }
+    retval = tanMp(x);
+    goto ret;
   }
 
   /* (III) The case 0.0608 < abs(x) <= 0.787 */
@@ -120,10 +129,10 @@ tan(double x) {
     z = w-xfg[i][0].d;  z2 = z*z;   s = (x<ZERO) ? MONE : ONE;
     pz = z+z*z2*(e0.d+z2*e1.d);
     fi = xfg[i][1].d;   gi = xfg[i][2].d;   t2 = pz*(gi+fi)/(gi-pz);
-    if ((y=fi+(t2-fi*u3.d))==fi+(t2+fi*u3.d))  return (s*y);
+    if ((y=fi+(t2-fi*u3.d))==fi+(t2+fi*u3.d)) { retval = (s*y); goto ret; }
     t3 = (t2<ZERO) ? -t2 : t2;
     t4 = fi*ua3.d+t3*ub3.d;
-    if ((y=fi+(t2-t4))==fi+(t2+t4))  return (s*y);
+    if ((y=fi+(t2-t4))==fi+(t2+t4)) { retval = (s*y); goto ret; }
 
     /* Second stage */
     ffi = xfg[i][3].d;
@@ -141,8 +150,9 @@ tan(double x) {
     SUB2(one.d,zero.d,c3,cc3,c1,cc1,t1,t2)
     DIV2(c2,cc2,c1,cc1,c3,cc3,t1,t2,t3,t4,t5,t6,t7,t8,t9,t10)
 
-    if ((y=c3+(cc3-u4.d*c3))==c3+(cc3+u4.d*c3))  return (s*y);
-    return tanMp(x);
+      if ((y=c3+(cc3-u4.d*c3))==c3+(cc3+u4.d*c3)) { retval = (s*y); goto ret; }
+    retval = tanMp(x);
+    goto ret;
   }
 
   /* (---) The case 0.787 < abs(x) <= 25 */
@@ -160,7 +170,7 @@ tan(double x) {
     else         {ya= a;  yya= da;  sy= ONE;}
 
     /* (IV),(V) The case 0.787 < abs(x) <= 25,    abs(y) <= 1e-7 */
-    if (ya<=gy1.d)  return tanMp(x);
+    if (ya<=gy1.d) { retval = tanMp(x); goto ret; }
 
     /* (VI) The case 0.787 < abs(x) <= 25,    1e-7 < abs(y) <= 0.0608 */
     if (ya<=gy2.d) {
@@ -170,10 +180,10 @@ tan(double x) {
 	/* First stage -cot */
 	EADD(a,t2,b,db)
 	DIV2(one.d,zero.d,b,db,c,dc,t1,t2,t3,t4,t5,t6,t7,t8,t9,t10)
-	if ((y=c+(dc-u6.d*c))==c+(dc+u6.d*c))  return (-y); }
+	if ((y=c+(dc-u6.d*c))==c+(dc+u6.d*c)) { retval = (-y); goto ret; } }
       else {
 	/* First stage tan */
-	if ((y=a+(t2-u5.d*a))==a+(t2+u5.d*a))  return y; }
+	if ((y=a+(t2-u5.d*a))==a+(t2+u5.d*a)) { retval = y; goto ret; } }
       /* Second stage */
       /* Range reduction by algorithm ii */
       t = (x*hpinv.d + toint.d);
@@ -211,11 +221,12 @@ tan(double x) {
       if (n) {
 	/* Second stage -cot */
 	DIV2(one.d,zero.d,c1,cc1,c2,cc2,t1,t2,t3,t4,t5,t6,t7,t8,t9,t10)
-	if ((y=c2+(cc2-u8.d*c2)) == c2+(cc2+u8.d*c2))  return (-y); }
+	if ((y=c2+(cc2-u8.d*c2)) == c2+(cc2+u8.d*c2)) { retval = (-y); goto ret; } }
       else {
 	/* Second stage tan */
-	if ((y=c1+(cc1-u7.d*c1)) == c1+(cc1+u7.d*c1))  return y; }
-      return tanMp(x);
+	if ((y=c1+(cc1-u7.d*c1)) == c1+(cc1+u7.d*c1)) { retval = y; goto ret; } }
+      retval = tanMp(x);
+      goto ret;
     }
 
     /* (VII) The case 0.787 < abs(x) <= 25,    0.0608 < abs(y) <= 0.787 */
@@ -229,17 +240,17 @@ tan(double x) {
     if (n) {
       /* -cot */
       t2 = pz*(fi+gi)/(fi+pz);
-      if ((y=gi-(t2-gi*u10.d))==gi-(t2+gi*u10.d))  return (-sy*y);
+      if ((y=gi-(t2-gi*u10.d))==gi-(t2+gi*u10.d)) { retval = (-sy*y); goto ret; }
       t3 = (t2<ZERO) ? -t2 : t2;
       t4 = gi*ua10.d+t3*ub10.d;
-      if ((y=gi-(t2-t4))==gi-(t2+t4))  return (-sy*y); }
+      if ((y=gi-(t2-t4))==gi-(t2+t4)) { retval = (-sy*y); goto ret; } }
     else   {
       /* tan */
       t2 = pz*(gi+fi)/(gi-pz);
-      if ((y=fi+(t2-fi*u9.d))==fi+(t2+fi*u9.d))  return (sy*y);
+      if ((y=fi+(t2-fi*u9.d))==fi+(t2+fi*u9.d)) { retval = (sy*y); goto ret; }
       t3 = (t2<ZERO) ? -t2 : t2;
       t4 = fi*ua9.d+t3*ub9.d;
-      if ((y=fi+(t2-t4))==fi+(t2+t4))  return (sy*y); }
+      if ((y=fi+(t2-t4))==fi+(t2+t4)) { retval = (sy*y); goto ret; } }
 
     /* Second stage */
     ffi = xfg[i][3].d;
@@ -260,13 +271,14 @@ tan(double x) {
     if (n) {
       /* -cot */
       DIV2(c1,cc1,c2,cc2,c3,cc3,t1,t2,t3,t4,t5,t6,t7,t8,t9,t10)
-      if ((y=c3+(cc3-u12.d*c3))==c3+(cc3+u12.d*c3))  return (-sy*y); }
+      if ((y=c3+(cc3-u12.d*c3))==c3+(cc3+u12.d*c3)) { retval = (-sy*y); goto ret; } }
     else {
       /* tan */
       DIV2(c2,cc2,c1,cc1,c3,cc3,t1,t2,t3,t4,t5,t6,t7,t8,t9,t10)
-      if ((y=c3+(cc3-u11.d*c3))==c3+(cc3+u11.d*c3))  return (sy*y); }
+      if ((y=c3+(cc3-u11.d*c3))==c3+(cc3+u11.d*c3)) { retval = (sy*y); goto ret; } }
 
-    return tanMp(x);
+    retval = tanMp(x);
+    goto ret;
   }
 
   /* (---) The case 25 < abs(x) <= 1e8 */
@@ -288,7 +300,7 @@ tan(double x) {
     else         {ya= a;  yya= da;  sy= ONE;}
 
     /* (+++) The case 25 < abs(x) <= 1e8,    abs(y) <= 1e-7 */
-    if (ya<=gy1.d)  return tanMp(x);
+    if (ya<=gy1.d) { retval = tanMp(x); goto ret; }
 
     /* (VIII) The case 25 < abs(x) <= 1e8,    1e-7 < abs(y) <= 0.0608 */
     if (ya<=gy2.d) {
@@ -298,10 +310,10 @@ tan(double x) {
 	/* First stage -cot */
 	EADD(a,t2,b,db)
 	DIV2(one.d,zero.d,b,db,c,dc,t1,t2,t3,t4,t5,t6,t7,t8,t9,t10)
-	if ((y=c+(dc-u14.d*c))==c+(dc+u14.d*c))  return (-y); }
+	if ((y=c+(dc-u14.d*c))==c+(dc+u14.d*c)) { retval = (-y); goto ret; } }
       else {
 	/* First stage tan */
-	if ((y=a+(t2-u13.d*a))==a+(t2+u13.d*a))  return y; }
+	if ((y=a+(t2-u13.d*a))==a+(t2+u13.d*a)) { retval = y; goto ret; } }
 
       /* Second stage */
       MUL2(a,da,a,da,x2,xx2,t1,t2,t3,t4,t5,t6,t7,t8)
@@ -325,11 +337,12 @@ tan(double x) {
       if (n) {
 	/* Second stage -cot */
 	DIV2(one.d,zero.d,c1,cc1,c2,cc2,t1,t2,t3,t4,t5,t6,t7,t8,t9,t10)
-	if ((y=c2+(cc2-u16.d*c2)) == c2+(cc2+u16.d*c2))  return (-y); }
+	if ((y=c2+(cc2-u16.d*c2)) == c2+(cc2+u16.d*c2)) { retval = (-y); goto ret; } }
       else {
 	/* Second stage tan */
-	if ((y=c1+(cc1-u15.d*c1)) == c1+(cc1+u15.d*c1))  return (y); }
-      return tanMp(x);
+	if ((y=c1+(cc1-u15.d*c1)) == c1+(cc1+u15.d*c1)) { retval = (y); goto ret; } }
+      retval = tanMp(x);
+      goto ret;
     }
 
     /* (IX) The case 25 < abs(x) <= 1e8,    0.0608 < abs(y) <= 0.787 */
@@ -342,17 +355,17 @@ tan(double x) {
     if (n) {
       /* -cot */
       t2 = pz*(fi+gi)/(fi+pz);
-      if ((y=gi-(t2-gi*u18.d))==gi-(t2+gi*u18.d))  return (-sy*y);
+      if ((y=gi-(t2-gi*u18.d))==gi-(t2+gi*u18.d)) { retval = (-sy*y); goto ret; }
       t3 = (t2<ZERO) ? -t2 : t2;
       t4 = gi*ua18.d+t3*ub18.d;
-      if ((y=gi-(t2-t4))==gi-(t2+t4))  return (-sy*y); }
+      if ((y=gi-(t2-t4))==gi-(t2+t4)) { retval = (-sy*y); goto ret; } }
     else   {
       /* tan */
       t2 = pz*(gi+fi)/(gi-pz);
-      if ((y=fi+(t2-fi*u17.d))==fi+(t2+fi*u17.d))  return (sy*y);
+      if ((y=fi+(t2-fi*u17.d))==fi+(t2+fi*u17.d)) { retval = (sy*y); goto ret; }
       t3 = (t2<ZERO) ? -t2 : t2;
       t4 = fi*ua17.d+t3*ub17.d;
-      if ((y=fi+(t2-t4))==fi+(t2+t4))  return (sy*y); }
+      if ((y=fi+(t2-t4))==fi+(t2+t4)) { retval = (sy*y); goto ret; } }
 
     /* Second stage */
     ffi = xfg[i][3].d;
@@ -373,12 +386,13 @@ tan(double x) {
     if (n) {
       /* -cot */
       DIV2(c1,cc1,c2,cc2,c3,cc3,t1,t2,t3,t4,t5,t6,t7,t8,t9,t10)
-      if ((y=c3+(cc3-u20.d*c3))==c3+(cc3+u20.d*c3))  return (-sy*y); }
+      if ((y=c3+(cc3-u20.d*c3))==c3+(cc3+u20.d*c3)) { retval = (-sy*y); goto ret; } }
     else {
       /* tan */
       DIV2(c2,cc2,c1,cc1,c3,cc3,t1,t2,t3,t4,t5,t6,t7,t8,t9,t10)
-      if ((y=c3+(cc3-u19.d*c3))==c3+(cc3+u19.d*c3))  return (sy*y); }
-    return tanMp(x);
+      if ((y=c3+(cc3-u19.d*c3))==c3+(cc3+u19.d*c3)) { retval = (sy*y); goto ret; } }
+    retval = tanMp(x);
+    goto ret;
   }
 
   /* (---) The case 1e8 < abs(x) < 2**1024 */
@@ -389,7 +403,7 @@ tan(double x) {
   else         {ya= a;  yya= da;  sy= ONE;}
 
   /* (+++) The case 1e8 < abs(x) < 2**1024,    abs(y) <= 1e-7 */
-  if (ya<=gy1.d)  return tanMp(x);
+  if (ya<=gy1.d) { retval = tanMp(x); goto ret; }
 
   /* (X) The case 1e8 < abs(x) < 2**1024,    1e-7 < abs(y) <= 0.0608 */
   if (ya<=gy2.d) {
@@ -399,10 +413,10 @@ tan(double x) {
       /* First stage -cot */
       EADD(a,t2,b,db)
       DIV2(one.d,zero.d,b,db,c,dc,t1,t2,t3,t4,t5,t6,t7,t8,t9,t10)
-      if ((y=c+(dc-u22.d*c))==c+(dc+u22.d*c))  return (-y); }
+      if ((y=c+(dc-u22.d*c))==c+(dc+u22.d*c)) { retval = (-y); goto ret; } }
     else {
       /* First stage tan */
-      if ((y=a+(t2-u21.d*a))==a+(t2+u21.d*a))  return y; }
+      if ((y=a+(t2-u21.d*a))==a+(t2+u21.d*a)) { retval = y; goto ret; } }
 
     /* Second stage */
     /* Reduction by algorithm iv */
@@ -431,11 +445,12 @@ tan(double x) {
     if (n) {
       /* Second stage -cot */
       DIV2(one.d,zero.d,c1,cc1,c2,cc2,t1,t2,t3,t4,t5,t6,t7,t8,t9,t10)
-      if ((y=c2+(cc2-u24.d*c2)) == c2+(cc2+u24.d*c2))  return (-y); }
+      if ((y=c2+(cc2-u24.d*c2)) == c2+(cc2+u24.d*c2)) { retval = (-y); goto ret; } }
     else {
       /* Second stage tan */
-      if ((y=c1+(cc1-u23.d*c1)) == c1+(cc1+u23.d*c1))  return y; }
-    return tanMp(x);
+      if ((y=c1+(cc1-u23.d*c1)) == c1+(cc1+u23.d*c1)) { retval = y; goto ret; } }
+    retval = tanMp(x);
+    goto ret;
   }
 
   /* (XI) The case 1e8 < abs(x) < 2**1024,    0.0608 < abs(y) <= 0.787 */
@@ -448,17 +463,17 @@ tan(double x) {
   if (n) {
     /* -cot */
     t2 = pz*(fi+gi)/(fi+pz);
-    if ((y=gi-(t2-gi*u26.d))==gi-(t2+gi*u26.d))  return (-sy*y);
+    if ((y=gi-(t2-gi*u26.d))==gi-(t2+gi*u26.d)) { retval = (-sy*y); goto ret; }
     t3 = (t2<ZERO) ? -t2 : t2;
     t4 = gi*ua26.d+t3*ub26.d;
-    if ((y=gi-(t2-t4))==gi-(t2+t4))  return (-sy*y); }
+    if ((y=gi-(t2-t4))==gi-(t2+t4)) { retval = (-sy*y); goto ret; } }
   else   {
     /* tan */
     t2 = pz*(gi+fi)/(gi-pz);
-    if ((y=fi+(t2-fi*u25.d))==fi+(t2+fi*u25.d))  return (sy*y);
+    if ((y=fi+(t2-fi*u25.d))==fi+(t2+fi*u25.d)) { retval = (sy*y); goto ret; }
     t3 = (t2<ZERO) ? -t2 : t2;
     t4 = fi*ua25.d+t3*ub25.d;
-    if ((y=fi+(t2-t4))==fi+(t2+t4))  return (sy*y); }
+    if ((y=fi+(t2-t4))==fi+(t2+t4)) { retval = (sy*y); goto ret; } }
 
   /* Second stage */
   ffi = xfg[i][3].d;
@@ -479,14 +494,18 @@ tan(double x) {
   if (n) {
     /* -cot */
     DIV2(c1,cc1,c2,cc2,c3,cc3,t1,t2,t3,t4,t5,t6,t7,t8,t9,t10)
-    if ((y=c3+(cc3-u28.d*c3))==c3+(cc3+u28.d*c3))  return (-sy*y); }
+    if ((y=c3+(cc3-u28.d*c3))==c3+(cc3+u28.d*c3)) { retval = (-sy*y); goto ret; } }
   else {
     /* tan */
     DIV2(c2,cc2,c1,cc1,c3,cc3,t1,t2,t3,t4,t5,t6,t7,t8,t9,t10)
-    if ((y=c3+(cc3-u27.d*c3))==c3+(cc3+u27.d*c3))  return (sy*y); }
-  return tanMp(x);
-}
+    if ((y=c3+(cc3-u27.d*c3))==c3+(cc3+u27.d*c3)) { retval = (sy*y); goto ret; } }
+  retval = tanMp(x);
+  goto ret;
 
+ ret:
+  libc_feupdateenv (&env);
+  return retval;
+}
 
 /* multiple precision stage                                              */
 /* Convert x to multi precision number,compute tan(x) by mptan() routine */