sysdeps/ieee754/dbl-64/s_tan.c

   1 /*
   2  * IBM Accurate Mathematical Library
   3  * written by International Business Machines Corp.
   4  * Copyright (C) 2001, 2009, 2011 Free Software Foundation
   5  *
   6  * This program is free software; you can redistribute it and/or modify
   7  * it under the terms of the GNU Lesser General Public License as published by
   8  * the Free Software Foundation; either version 2.1 of the License, or
   9  * (at your option) any later version.
  10  *
  11  * This program is distributed in the hope that it will be useful,
  12  * but WITHOUT ANY WARRANTY; without even the implied warranty of
  13  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
  14  * GNU Lesser General Public License for more details.
  15  *
  16  * You should have received a copy of the GNU Lesser General Public License
  17  * along with this program; if not, write to the Free Software
  18  * Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA.
  19  */
  20 /*********************************************************************/
  21 /*  MODULE_NAME: utan.c                                              */
  22 /*                                                                   */
  23 /*  FUNCTIONS: utan                                                  */
  24 /*             tanMp                                                 */
  25 /*                                                                   */
  26 /*  FILES NEEDED:dla.h endian.h mpa.h mydefs.h utan.h                */
  27 /*               branred.c sincos32.c mptan.c                        */
  28 /*               utan.tbl                                            */
  29 /*                                                                   */
  30 /* An ultimate tan routine. Given an IEEE double machine number x    */
  31 /* it computes the correctly rounded (to nearest) value of tan(x).   */
  32 /* Assumption: Machine arithmetic operations are performed in        */
  33 /* round to nearest mode of IEEE 754 standard.                       */
  34 /*                                                                   */
  35 /*********************************************************************/
  36
  37 #include <errno.h>
  38 #include "endian.h"
  39 #include "dla.h"
  40 #include "mpa.h"
  41 #include "MathLib.h"
  42 #include "math.h"
  43
  44 static double tanMp(double);
  45 void __mptan(double, mp_no *, int);
  46
  47 double tan(double x) {
  48 #include "utan.h"
  49 #include "utan.tbl"
  50
  51   int ux,i,n;
  52   double a,da,a2,b,db,c,dc,c1,cc1,c2,cc2,c3,cc3,fi,ffi,gi,pz,s,sy,
  53   t,t1,t2,t3,t4,t7,t8,t9,t10,w,x2,xn,xx2,y,ya,yya,z0,z,zz,z2,zz2;
  54 #ifndef DLA_FMA
  55   double t5,t6;
  56 #endif
  57   int p;
  58   number num,v;
  59   mp_no mpa,mpt1,mpt2;
  60 #if 0
  61   mp_no mpy;
  62 #endif
  63
  64   int __branred(double, double *, double *);
  65   int __mpranred(double, mp_no *, int);
  66
  67   /* x=+-INF, x=NaN */
  68   num.d = x;  ux = num.i[HIGH_HALF];
  69   if ((ux&0x7ff00000)==0x7ff00000) {
  70     if ((ux&0x7fffffff)==0x7ff00000)
  71       __set_errno (EDOM);
  72     return x-x;
  73   }
  74
  75   w=(x<ZERO) ? -x : x;
  76
  77   /* (I) The case abs(x) <= 1.259e-8 */
  78   if (w<=g1.d)  return x;
  79
  80   /* (II) The case 1.259e-8 < abs(x) <= 0.0608 */
  81   if (w<=g2.d) {
  82
  83     /* First stage */
  84     x2 = x*x;
  85     t2 = x*x2*(d3.d+x2*(d5.d+x2*(d7.d+x2*(d9.d+x2*d11.d))));
  86     if ((y=x+(t2-u1.d*t2)) == x+(t2+u1.d*t2))  return y;
  87
  88     /* Second stage */
  89     c1 = x2*(a15.d+x2*(a17.d+x2*(a19.d+x2*(a21.d+x2*(a23.d+x2*(a25.d+
  90          x2*a27.d))))));
  91     EMULV(x,x,x2,xx2,t1,t2,t3,t4,t5)
  92     ADD2(a13.d,aa13.d,c1,zero.d,c2,cc2,t1,t2)
  93     MUL2(x2,xx2,c2,cc2,c1,cc1,t1,t2,t3,t4,t5,t6,t7,t8)
  94     ADD2(a11.d,aa11.d,c1,cc1,c2,cc2,t1,t2)
  95     MUL2(x2,xx2,c2,cc2,c1,cc1,t1,t2,t3,t4,t5,t6,t7,t8)
  96     ADD2(a9.d ,aa9.d ,c1,cc1,c2,cc2,t1,t2)
  97     MUL2(x2,xx2,c2,cc2,c1,cc1,t1,t2,t3,t4,t5,t6,t7,t8)
  98     ADD2(a7.d ,aa7.d ,c1,cc1,c2,cc2,t1,t2)
  99     MUL2(x2,xx2,c2,cc2,c1,cc1,t1,t2,t3,t4,t5,t6,t7,t8)
 100     ADD2(a5.d ,aa5.d ,c1,cc1,c2,cc2,t1,t2)
 101     MUL2(x2,xx2,c2,cc2,c1,cc1,t1,t2,t3,t4,t5,t6,t7,t8)
 102     ADD2(a3.d ,aa3.d ,c1,cc1,c2,cc2,t1,t2)
 103     MUL2(x2,xx2,c2,cc2,c1,cc1,t1,t2,t3,t4,t5,t6,t7,t8)
 104     MUL2(x ,zero.d,c1,cc1,c2,cc2,t1,t2,t3,t4,t5,t6,t7,t8)
 105     ADD2(x    ,zero.d,c2,cc2,c1,cc1,t1,t2)
 106     if ((y=c1+(cc1-u2.d*c1)) == c1+(cc1+u2.d*c1))  return y;
 107     return tanMp(x);
 108   }
 109
 110   /* (III) The case 0.0608 < abs(x) <= 0.787 */
 111   if (w<=g3.d) {
 112
 113     /* First stage */
 114     i = ((int) (mfftnhf.d+TWO8*w));
 115     z = w-xfg[i][0].d;  z2 = z*z;   s = (x<ZERO) ? MONE : ONE;
 116     pz = z+z*z2*(e0.d+z2*e1.d);
 117     fi = xfg[i][1].d;   gi = xfg[i][2].d;   t2 = pz*(gi+fi)/(gi-pz);
 118     if ((y=fi+(t2-fi*u3.d))==fi+(t2+fi*u3.d))  return (s*y);
 119     t3 = (t2<ZERO) ? -t2 : t2;
 120     t4 = fi*ua3.d+t3*ub3.d;
 121     if ((y=fi+(t2-t4))==fi+(t2+t4))  return (s*y);
 122
 123     /* Second stage */
 124     ffi = xfg[i][3].d;
 125     c1 = z2*(a7.d+z2*(a9.d+z2*a11.d));
 126     EMULV(z,z,z2,zz2,t1,t2,t3,t4,t5)
 127     ADD2(a5.d,aa5.d,c1,zero.d,c2,cc2,t1,t2)
 128     MUL2(z2,zz2,c2,cc2,c1,cc1,t1,t2,t3,t4,t5,t6,t7,t8)
 129     ADD2(a3.d,aa3.d,c1,cc1,c2,cc2,t1,t2)
 130     MUL2(z2,zz2,c2,cc2,c1,cc1,t1,t2,t3,t4,t5,t6,t7,t8)
 131     MUL2(z ,zero.d,c1,cc1,c2,cc2,t1,t2,t3,t4,t5,t6,t7,t8)
 132     ADD2(z ,zero.d,c2,cc2,c1,cc1,t1,t2)
 133
 134     ADD2(fi ,ffi,c1,cc1,c2,cc2,t1,t2)
 135     MUL2(fi ,ffi,c1,cc1,c3,cc3,t1,t2,t3,t4,t5,t6,t7,t8)
 136     SUB2(one.d,zero.d,c3,cc3,c1,cc1,t1,t2)
 137     DIV2(c2,cc2,c1,cc1,c3,cc3,t1,t2,t3,t4,t5,t6,t7,t8,t9,t10)
 138
 139     if ((y=c3+(cc3-u4.d*c3))==c3+(cc3+u4.d*c3))  return (s*y);
 140     return tanMp(x);
 141   }
 142
 143   /* (---) The case 0.787 < abs(x) <= 25 */
 144   if (w<=g4.d) {
 145     /* Range reduction by algorithm i */
 146     t = (x*hpinv.d + toint.d);
 147     xn = t - toint.d;
 148     v.d = t;
 149     t1 = (x - xn*mp1.d) - xn*mp2.d;
 150     n =v.i[LOW_HALF] & 0x00000001;
 151     da = xn*mp3.d;
 152     a=t1-da;
 153     da = (t1-a)-da;
 154     if (a<ZERO)  {ya=-a;  yya=-da;  sy=MONE;}
 155     else         {ya= a;  yya= da;  sy= ONE;}
 156
 157     /* (IV),(V) The case 0.787 < abs(x) <= 25,    abs(y) <= 1e-7 */
 158     if (ya<=gy1.d)  return tanMp(x);
 159
 160     /* (VI) The case 0.787 < abs(x) <= 25,    1e-7 < abs(y) <= 0.0608 */
 161     if (ya<=gy2.d) {
 162       a2 = a*a;
 163       t2 = da+a*a2*(d3.d+a2*(d5.d+a2*(d7.d+a2*(d9.d+a2*d11.d))));
 164       if (n) {
 165         /* First stage -cot */
 166         EADD(a,t2,b,db)
 167         DIV2(one.d,zero.d,b,db,c,dc,t1,t2,t3,t4,t5,t6,t7,t8,t9,t10)
 168         if ((y=c+(dc-u6.d*c))==c+(dc+u6.d*c))  return (-y); }
 169       else {
 170         /* First stage tan */
 171         if ((y=a+(t2-u5.d*a))==a+(t2+u5.d*a))  return y; }
 172       /* Second stage */
 173       /* Range reduction by algorithm ii */
 174       t = (x*hpinv.d + toint.d);
 175       xn = t - toint.d;
 176       v.d = t;
 177       t1 = (x - xn*mp1.d) - xn*mp2.d;
 178       n =v.i[LOW_HALF] & 0x00000001;
 179       da = xn*pp3.d;
 180       t=t1-da;
 181       da = (t1-t)-da;
 182       t1 = xn*pp4.d;
 183       a = t - t1;
 184       da = ((t-a)-t1)+da;
 185
 186       /* Second stage */
 187       EADD(a,da,t1,t2)   a=t1;  da=t2;
 188       MUL2(a,da,a,da,x2,xx2,t1,t2,t3,t4,t5,t6,t7,t8)
 189       c1 = x2*(a15.d+x2*(a17.d+x2*(a19.d+x2*(a21.d+x2*(a23.d+x2*(a25.d+
 190            x2*a27.d))))));
 191       ADD2(a13.d,aa13.d,c1,zero.d,c2,cc2,t1,t2)
 192       MUL2(x2,xx2,c2,cc2,c1,cc1,t1,t2,t3,t4,t5,t6,t7,t8)
 193       ADD2(a11.d,aa11.d,c1,cc1,c2,cc2,t1,t2)
 194       MUL2(x2,xx2,c2,cc2,c1,cc1,t1,t2,t3,t4,t5,t6,t7,t8)
 195       ADD2(a9.d ,aa9.d ,c1,cc1,c2,cc2,t1,t2)
 196       MUL2(x2,xx2,c2,cc2,c1,cc1,t1,t2,t3,t4,t5,t6,t7,t8)
 197       ADD2(a7.d ,aa7.d ,c1,cc1,c2,cc2,t1,t2)
 198       MUL2(x2,xx2,c2,cc2,c1,cc1,t1,t2,t3,t4,t5,t6,t7,t8)
 199       ADD2(a5.d ,aa5.d ,c1,cc1,c2,cc2,t1,t2)
 200       MUL2(x2,xx2,c2,cc2,c1,cc1,t1,t2,t3,t4,t5,t6,t7,t8)
 201       ADD2(a3.d ,aa3.d ,c1,cc1,c2,cc2,t1,t2)
 202       MUL2(x2,xx2,c2,cc2,c1,cc1,t1,t2,t3,t4,t5,t6,t7,t8)
 203       MUL2(a ,da ,c1,cc1,c2,cc2,t1,t2,t3,t4,t5,t6,t7,t8)
 204       ADD2(a  ,da  ,c2,cc2,c1,cc1,t1,t2)
 205
 206       if (n) {
 207         /* Second stage -cot */
 208         DIV2(one.d,zero.d,c1,cc1,c2,cc2,t1,t2,t3,t4,t5,t6,t7,t8,t9,t10)
 209         if ((y=c2+(cc2-u8.d*c2)) == c2+(cc2+u8.d*c2))  return (-y); }
 210       else {
 211         /* Second stage tan */
 212         if ((y=c1+(cc1-u7.d*c1)) == c1+(cc1+u7.d*c1))  return y; }
 213       return tanMp(x);
 214     }
 215
 216     /* (VII) The case 0.787 < abs(x) <= 25,    0.0608 < abs(y) <= 0.787 */
 217
 218     /* First stage */
 219     i = ((int) (mfftnhf.d+TWO8*ya));
 220     z = (z0=(ya-xfg[i][0].d))+yya;  z2 = z*z;
 221     pz = z+z*z2*(e0.d+z2*e1.d);
 222     fi = xfg[i][1].d;   gi = xfg[i][2].d;
 223
 224     if (n) {
 225       /* -cot */
 226       t2 = pz*(fi+gi)/(fi+pz);
 227       if ((y=gi-(t2-gi*u10.d))==gi-(t2+gi*u10.d))  return (-sy*y);
 228       t3 = (t2<ZERO) ? -t2 : t2;
 229       t4 = gi*ua10.d+t3*ub10.d;
 230       if ((y=gi-(t2-t4))==gi-(t2+t4))  return (-sy*y); }
 231     else   {
 232       /* tan */
 233       t2 = pz*(gi+fi)/(gi-pz);
 234       if ((y=fi+(t2-fi*u9.d))==fi+(t2+fi*u9.d))  return (sy*y);
 235       t3 = (t2<ZERO) ? -t2 : t2;
 236       t4 = fi*ua9.d+t3*ub9.d;
 237       if ((y=fi+(t2-t4))==fi+(t2+t4))  return (sy*y); }
 238
 239     /* Second stage */
 240     ffi = xfg[i][3].d;
 241     EADD(z0,yya,z,zz)
 242     MUL2(z,zz,z,zz,z2,zz2,t1,t2,t3,t4,t5,t6,t7,t8)
 243     c1 = z2*(a7.d+z2*(a9.d+z2*a11.d));
 244     ADD2(a5.d,aa5.d,c1,zero.d,c2,cc2,t1,t2)
 245     MUL2(z2,zz2,c2,cc2,c1,cc1,t1,t2,t3,t4,t5,t6,t7,t8)
 246     ADD2(a3.d,aa3.d,c1,cc1,c2,cc2,t1,t2)
 247     MUL2(z2,zz2,c2,cc2,c1,cc1,t1,t2,t3,t4,t5,t6,t7,t8)
 248     MUL2(z ,zz ,c1,cc1,c2,cc2,t1,t2,t3,t4,t5,t6,t7,t8)
 249     ADD2(z ,zz ,c2,cc2,c1,cc1,t1,t2)
 250
 251     ADD2(fi ,ffi,c1,cc1,c2,cc2,t1,t2)
 252     MUL2(fi ,ffi,c1,cc1,c3,cc3,t1,t2,t3,t4,t5,t6,t7,t8)
 253     SUB2(one.d,zero.d,c3,cc3,c1,cc1,t1,t2)
 254
 255     if (n) {
 256       /* -cot */
 257       DIV2(c1,cc1,c2,cc2,c3,cc3,t1,t2,t3,t4,t5,t6,t7,t8,t9,t10)
 258       if ((y=c3+(cc3-u12.d*c3))==c3+(cc3+u12.d*c3))  return (-sy*y); }
 259     else {
 260       /* tan */
 261       DIV2(c2,cc2,c1,cc1,c3,cc3,t1,t2,t3,t4,t5,t6,t7,t8,t9,t10)
 262       if ((y=c3+(cc3-u11.d*c3))==c3+(cc3+u11.d*c3))  return (sy*y); }
 263
 264     return tanMp(x);
 265   }
 266
 267   /* (---) The case 25 < abs(x) <= 1e8 */
 268   if (w<=g5.d) {
 269     /* Range reduction by algorithm ii */
 270     t = (x*hpinv.d + toint.d);
 271     xn = t - toint.d;
 272     v.d = t;
 273     t1 = (x - xn*mp1.d) - xn*mp2.d;
 274     n =v.i[LOW_HALF] & 0x00000001;
 275     da = xn*pp3.d;
 276     t=t1-da;
 277     da = (t1-t)-da;
 278     t1 = xn*pp4.d;
 279     a = t - t1;
 280     da = ((t-a)-t1)+da;
 281     EADD(a,da,t1,t2)   a=t1;  da=t2;
 282     if (a<ZERO)  {ya=-a;  yya=-da;  sy=MONE;}
 283     else         {ya= a;  yya= da;  sy= ONE;}
 284
 285     /* (+++) The case 25 < abs(x) <= 1e8,    abs(y) <= 1e-7 */
 286     if (ya<=gy1.d)  return tanMp(x);
 287
 288     /* (VIII) The case 25 < abs(x) <= 1e8,    1e-7 < abs(y) <= 0.0608 */
 289     if (ya<=gy2.d) {
 290       a2 = a*a;
 291       t2 = da+a*a2*(d3.d+a2*(d5.d+a2*(d7.d+a2*(d9.d+a2*d11.d))));
 292       if (n) {
 293         /* First stage -cot */
 294         EADD(a,t2,b,db)
 295         DIV2(one.d,zero.d,b,db,c,dc,t1,t2,t3,t4,t5,t6,t7,t8,t9,t10)
 296         if ((y=c+(dc-u14.d*c))==c+(dc+u14.d*c))  return (-y); }
 297       else {
 298         /* First stage tan */
 299         if ((y=a+(t2-u13.d*a))==a+(t2+u13.d*a))  return y; }
 300
 301       /* Second stage */
 302       MUL2(a,da,a,da,x2,xx2,t1,t2,t3,t4,t5,t6,t7,t8)
 303       c1 = x2*(a15.d+x2*(a17.d+x2*(a19.d+x2*(a21.d+x2*(a23.d+x2*(a25.d+
 304            x2*a27.d))))));
 305       ADD2(a13.d,aa13.d,c1,zero.d,c2,cc2,t1,t2)
 306       MUL2(x2,xx2,c2,cc2,c1,cc1,t1,t2,t3,t4,t5,t6,t7,t8)
 307       ADD2(a11.d,aa11.d,c1,cc1,c2,cc2,t1,t2)
 308       MUL2(x2,xx2,c2,cc2,c1,cc1,t1,t2,t3,t4,t5,t6,t7,t8)
 309       ADD2(a9.d ,aa9.d ,c1,cc1,c2,cc2,t1,t2)
 310       MUL2(x2,xx2,c2,cc2,c1,cc1,t1,t2,t3,t4,t5,t6,t7,t8)
 311       ADD2(a7.d ,aa7.d ,c1,cc1,c2,cc2,t1,t2)
 312       MUL2(x2,xx2,c2,cc2,c1,cc1,t1,t2,t3,t4,t5,t6,t7,t8)
 313       ADD2(a5.d ,aa5.d ,c1,cc1,c2,cc2,t1,t2)
 314       MUL2(x2,xx2,c2,cc2,c1,cc1,t1,t2,t3,t4,t5,t6,t7,t8)
 315       ADD2(a3.d ,aa3.d ,c1,cc1,c2,cc2,t1,t2)
 316       MUL2(x2,xx2,c2,cc2,c1,cc1,t1,t2,t3,t4,t5,t6,t7,t8)
 317       MUL2(a ,da ,c1,cc1,c2,cc2,t1,t2,t3,t4,t5,t6,t7,t8)
 318       ADD2(a  ,da  ,c2,cc2,c1,cc1,t1,t2)
 319
 320       if (n) {
 321         /* Second stage -cot */
 322         DIV2(one.d,zero.d,c1,cc1,c2,cc2,t1,t2,t3,t4,t5,t6,t7,t8,t9,t10)
 323         if ((y=c2+(cc2-u16.d*c2)) == c2+(cc2+u16.d*c2))  return (-y); }
 324       else {
 325         /* Second stage tan */
 326         if ((y=c1+(cc1-u15.d*c1)) == c1+(cc1+u15.d*c1))  return (y); }
 327       return tanMp(x);
 328     }
 329
 330     /* (IX) The case 25 < abs(x) <= 1e8,    0.0608 < abs(y) <= 0.787 */
 331     /* First stage */
 332     i = ((int) (mfftnhf.d+TWO8*ya));
 333     z = (z0=(ya-xfg[i][0].d))+yya;  z2 = z*z;
 334     pz = z+z*z2*(e0.d+z2*e1.d);
 335     fi = xfg[i][1].d;   gi = xfg[i][2].d;
 336
 337     if (n) {
 338       /* -cot */
 339       t2 = pz*(fi+gi)/(fi+pz);
 340       if ((y=gi-(t2-gi*u18.d))==gi-(t2+gi*u18.d))  return (-sy*y);
 341       t3 = (t2<ZERO) ? -t2 : t2;
 342       t4 = gi*ua18.d+t3*ub18.d;
 343       if ((y=gi-(t2-t4))==gi-(t2+t4))  return (-sy*y); }
 344     else   {
 345       /* tan */
 346       t2 = pz*(gi+fi)/(gi-pz);
 347       if ((y=fi+(t2-fi*u17.d))==fi+(t2+fi*u17.d))  return (sy*y);
 348       t3 = (t2<ZERO) ? -t2 : t2;
 349       t4 = fi*ua17.d+t3*ub17.d;
 350       if ((y=fi+(t2-t4))==fi+(t2+t4))  return (sy*y); }
 351
 352     /* Second stage */
 353     ffi = xfg[i][3].d;
 354     EADD(z0,yya,z,zz)
 355     MUL2(z,zz,z,zz,z2,zz2,t1,t2,t3,t4,t5,t6,t7,t8)
 356     c1 = z2*(a7.d+z2*(a9.d+z2*a11.d));
 357     ADD2(a5.d,aa5.d,c1,zero.d,c2,cc2,t1,t2)
 358     MUL2(z2,zz2,c2,cc2,c1,cc1,t1,t2,t3,t4,t5,t6,t7,t8)
 359     ADD2(a3.d,aa3.d,c1,cc1,c2,cc2,t1,t2)
 360     MUL2(z2,zz2,c2,cc2,c1,cc1,t1,t2,t3,t4,t5,t6,t7,t8)
 361     MUL2(z ,zz ,c1,cc1,c2,cc2,t1,t2,t3,t4,t5,t6,t7,t8)
 362     ADD2(z ,zz ,c2,cc2,c1,cc1,t1,t2)
 363
 364     ADD2(fi ,ffi,c1,cc1,c2,cc2,t1,t2)
 365     MUL2(fi ,ffi,c1,cc1,c3,cc3,t1,t2,t3,t4,t5,t6,t7,t8)
 366     SUB2(one.d,zero.d,c3,cc3,c1,cc1,t1,t2)
 367
 368     if (n) {
 369       /* -cot */
 370       DIV2(c1,cc1,c2,cc2,c3,cc3,t1,t2,t3,t4,t5,t6,t7,t8,t9,t10)
 371       if ((y=c3+(cc3-u20.d*c3))==c3+(cc3+u20.d*c3))  return (-sy*y); }
 372     else {
 373       /* tan */
 374       DIV2(c2,cc2,c1,cc1,c3,cc3,t1,t2,t3,t4,t5,t6,t7,t8,t9,t10)
 375       if ((y=c3+(cc3-u19.d*c3))==c3+(cc3+u19.d*c3))  return (sy*y); }
 376     return tanMp(x);
 377   }
 378
 379   /* (---) The case 1e8 < abs(x) < 2**1024 */
 380   /* Range reduction by algorithm iii */
 381   n = (__branred(x,&a,&da)) & 0x00000001;
 382   EADD(a,da,t1,t2)   a=t1;  da=t2;
 383   if (a<ZERO)  {ya=-a;  yya=-da;  sy=MONE;}
 384   else         {ya= a;  yya= da;  sy= ONE;}
 385
 386   /* (+++) The case 1e8 < abs(x) < 2**1024,    abs(y) <= 1e-7 */
 387   if (ya<=gy1.d)  return tanMp(x);
 388
 389   /* (X) The case 1e8 < abs(x) < 2**1024,    1e-7 < abs(y) <= 0.0608 */
 390   if (ya<=gy2.d) {
 391     a2 = a*a;
 392     t2 = da+a*a2*(d3.d+a2*(d5.d+a2*(d7.d+a2*(d9.d+a2*d11.d))));
 393     if (n) {
 394       /* First stage -cot */
 395       EADD(a,t2,b,db)
 396       DIV2(one.d,zero.d,b,db,c,dc,t1,t2,t3,t4,t5,t6,t7,t8,t9,t10)
 397       if ((y=c+(dc-u22.d*c))==c+(dc+u22.d*c))  return (-y); }
 398     else {
 399       /* First stage tan */
 400       if ((y=a+(t2-u21.d*a))==a+(t2+u21.d*a))  return y; }
 401
 402     /* Second stage */
 403     /* Reduction by algorithm iv */
 404     p=10;    n = (__mpranred(x,&mpa,p)) & 0x00000001;
 405     __mp_dbl(&mpa,&a,p);        __dbl_mp(a,&mpt1,p);
 406     __sub(&mpa,&mpt1,&mpt2,p);  __mp_dbl(&mpt2,&da,p);
 407
 408     MUL2(a,da,a,da,x2,xx2,t1,t2,t3,t4,t5,t6,t7,t8)
 409     c1 = x2*(a15.d+x2*(a17.d+x2*(a19.d+x2*(a21.d+x2*(a23.d+x2*(a25.d+
 410          x2*a27.d))))));
 411     ADD2(a13.d,aa13.d,c1,zero.d,c2,cc2,t1,t2)
 412     MUL2(x2,xx2,c2,cc2,c1,cc1,t1,t2,t3,t4,t5,t6,t7,t8)
 413     ADD2(a11.d,aa11.d,c1,cc1,c2,cc2,t1,t2)
 414     MUL2(x2,xx2,c2,cc2,c1,cc1,t1,t2,t3,t4,t5,t6,t7,t8)
 415     ADD2(a9.d ,aa9.d ,c1,cc1,c2,cc2,t1,t2)
 416     MUL2(x2,xx2,c2,cc2,c1,cc1,t1,t2,t3,t4,t5,t6,t7,t8)
 417     ADD2(a7.d ,aa7.d ,c1,cc1,c2,cc2,t1,t2)
 418     MUL2(x2,xx2,c2,cc2,c1,cc1,t1,t2,t3,t4,t5,t6,t7,t8)
 419     ADD2(a5.d ,aa5.d ,c1,cc1,c2,cc2,t1,t2)
 420     MUL2(x2,xx2,c2,cc2,c1,cc1,t1,t2,t3,t4,t5,t6,t7,t8)
 421     ADD2(a3.d ,aa3.d ,c1,cc1,c2,cc2,t1,t2)
 422     MUL2(x2,xx2,c2,cc2,c1,cc1,t1,t2,t3,t4,t5,t6,t7,t8)
 423     MUL2(a ,da ,c1,cc1,c2,cc2,t1,t2,t3,t4,t5,t6,t7,t8)
 424     ADD2(a    ,da    ,c2,cc2,c1,cc1,t1,t2)
 425
 426     if (n) {
 427       /* Second stage -cot */
 428       DIV2(one.d,zero.d,c1,cc1,c2,cc2,t1,t2,t3,t4,t5,t6,t7,t8,t9,t10)
 429       if ((y=c2+(cc2-u24.d*c2)) == c2+(cc2+u24.d*c2))  return (-y); }
 430     else {
 431       /* Second stage tan */
 432       if ((y=c1+(cc1-u23.d*c1)) == c1+(cc1+u23.d*c1))  return y; }
 433     return tanMp(x);
 434   }
 435
 436   /* (XI) The case 1e8 < abs(x) < 2**1024,    0.0608 < abs(y) <= 0.787 */
 437   /* First stage */
 438   i = ((int) (mfftnhf.d+TWO8*ya));
 439   z = (z0=(ya-xfg[i][0].d))+yya;  z2 = z*z;
 440   pz = z+z*z2*(e0.d+z2*e1.d);
 441   fi = xfg[i][1].d;   gi = xfg[i][2].d;
 442
 443   if (n) {
 444     /* -cot */
 445     t2 = pz*(fi+gi)/(fi+pz);
 446     if ((y=gi-(t2-gi*u26.d))==gi-(t2+gi*u26.d))  return (-sy*y);
 447     t3 = (t2<ZERO) ? -t2 : t2;
 448     t4 = gi*ua26.d+t3*ub26.d;
 449     if ((y=gi-(t2-t4))==gi-(t2+t4))  return (-sy*y); }
 450   else   {
 451     /* tan */
 452     t2 = pz*(gi+fi)/(gi-pz);
 453     if ((y=fi+(t2-fi*u25.d))==fi+(t2+fi*u25.d))  return (sy*y);
 454     t3 = (t2<ZERO) ? -t2 : t2;
 455     t4 = fi*ua25.d+t3*ub25.d;
 456     if ((y=fi+(t2-t4))==fi+(t2+t4))  return (sy*y); }
 457
 458   /* Second stage */
 459   ffi = xfg[i][3].d;
 460   EADD(z0,yya,z,zz)
 461   MUL2(z,zz,z,zz,z2,zz2,t1,t2,t3,t4,t5,t6,t7,t8)
 462   c1 = z2*(a7.d+z2*(a9.d+z2*a11.d));
 463   ADD2(a5.d,aa5.d,c1,zero.d,c2,cc2,t1,t2)
 464   MUL2(z2,zz2,c2,cc2,c1,cc1,t1,t2,t3,t4,t5,t6,t7,t8)
 465   ADD2(a3.d,aa3.d,c1,cc1,c2,cc2,t1,t2)
 466   MUL2(z2,zz2,c2,cc2,c1,cc1,t1,t2,t3,t4,t5,t6,t7,t8)
 467   MUL2(z ,zz ,c1,cc1,c2,cc2,t1,t2,t3,t4,t5,t6,t7,t8)
 468   ADD2(z ,zz ,c2,cc2,c1,cc1,t1,t2)
 469
 470   ADD2(fi ,ffi,c1,cc1,c2,cc2,t1,t2)
 471   MUL2(fi ,ffi,c1,cc1,c3,cc3,t1,t2,t3,t4,t5,t6,t7,t8)
 472   SUB2(one.d,zero.d,c3,cc3,c1,cc1,t1,t2)
 473
 474   if (n) {
 475     /* -cot */
 476     DIV2(c1,cc1,c2,cc2,c3,cc3,t1,t2,t3,t4,t5,t6,t7,t8,t9,t10)
 477     if ((y=c3+(cc3-u28.d*c3))==c3+(cc3+u28.d*c3))  return (-sy*y); }
 478   else {
 479     /* tan */
 480     DIV2(c2,cc2,c1,cc1,c3,cc3,t1,t2,t3,t4,t5,t6,t7,t8,t9,t10)
 481     if ((y=c3+(cc3-u27.d*c3))==c3+(cc3+u27.d*c3))  return (sy*y); }
 482   return tanMp(x);
 483 }
 484
 485
 486 /* multiple precision stage                                              */
 487 /* Convert x to multi precision number,compute tan(x) by mptan() routine */
 488 /* and converts result back to double                                    */
 489 static double tanMp(double x)
 490 {
 491   int p;
 492   double y;
 493   mp_no mpy;
 494   p=32;
 495   __mptan(x, &mpy, p);
 496   __mp_dbl(&mpy,&y,p);
 497   return y;
 498 }
 499
 500 #ifdef NO_LONG_DOUBLE
 501 weak_alias (tan, tanl)
 502 #endif