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