Back to index

glibc  2.9
s_tan.c
Go to the documentation of this file.
00001 /*
00002  * IBM Accurate Mathematical Library
00003  * written by International Business Machines Corp.
00004  * Copyright (C) 2001 Free Software Foundation
00005  *
00006  * This program is free software; you can redistribute it and/or modify
00007  * it under the terms of the GNU Lesser General Public License as published by
00008  * the Free Software Foundation; either version 2.1 of the License, or
00009  * (at your option) any later version.
00010  *
00011  * This program is distributed in the hope that it will be useful,
00012  * but WITHOUT ANY WARRANTY; without even the implied warranty of
00013  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
00014  * GNU Lesser General Public License for more details.
00015  *
00016  * You should have received a copy of the GNU Lesser General Public License
00017  * along with this program; if not, write to the Free Software
00018  * Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA.
00019  */
00020 /*********************************************************************/
00021 /*  MODULE_NAME: utan.c                                              */
00022 /*                                                                   */
00023 /*  FUNCTIONS: utan                                                  */
00024 /*             tanMp                                                 */
00025 /*                                                                   */
00026 /*  FILES NEEDED:dla.h endian.h mpa.h mydefs.h utan.h                */
00027 /*               branred.c sincos32.c mptan.c                        */
00028 /*               utan.tbl                                            */
00029 /*                                                                   */
00030 /* An ultimate tan routine. Given an IEEE double machine number x    */
00031 /* it computes the correctly rounded (to nearest) value of tan(x).   */
00032 /* Assumption: Machine arithmetic operations are performed in        */
00033 /* round to nearest mode of IEEE 754 standard.                       */
00034 /*                                                                   */
00035 /*********************************************************************/
00036 #include "endian.h"
00037 #include "dla.h"
00038 #include "mpa.h"
00039 #include "MathLib.h"
00040 #include "math.h"
00041 
00042 static double tanMp(double);
00043 void __mptan(double, mp_no *, int);
00044 
00045 double tan(double x) {
00046 #include "utan.h"
00047 #include "utan.tbl"
00048 
00049   int ux,i,n;
00050   double a,da,a2,b,db,c,dc,c1,cc1,c2,cc2,c3,cc3,fi,ffi,gi,pz,s,sy,
00051   t,t1,t2,t3,t4,t5,t6,t7,t8,t9,t10,w,x2,xn,xx2,y,ya,yya,z0,z,zz,z2,zz2;
00052   int p;
00053   number num,v;
00054   mp_no mpa,mpt1,mpt2;
00055 #if 0
00056   mp_no mpy;
00057 #endif
00058 
00059   int __branred(double, double *, double *);
00060   int __mpranred(double, mp_no *, int);
00061 
00062   /* x=+-INF, x=NaN */
00063   num.d = x;  ux = num.i[HIGH_HALF];
00064   if ((ux&0x7ff00000)==0x7ff00000) return x-x;
00065 
00066   w=(x<ZERO) ? -x : x;
00067 
00068   /* (I) The case abs(x) <= 1.259e-8 */
00069   if (w<=g1.d)  return x;
00070 
00071   /* (II) The case 1.259e-8 < abs(x) <= 0.0608 */
00072   if (w<=g2.d) {
00073 
00074     /* First stage */
00075     x2 = x*x;
00076     t2 = x*x2*(d3.d+x2*(d5.d+x2*(d7.d+x2*(d9.d+x2*d11.d))));
00077     if ((y=x+(t2-u1.d*t2)) == x+(t2+u1.d*t2))  return y;
00078 
00079     /* Second stage */
00080     c1 = x2*(a15.d+x2*(a17.d+x2*(a19.d+x2*(a21.d+x2*(a23.d+x2*(a25.d+
00081          x2*a27.d))))));
00082     EMULV(x,x,x2,xx2,t1,t2,t3,t4,t5)
00083     ADD2(a13.d,aa13.d,c1,zero.d,c2,cc2,t1,t2)
00084     MUL2(x2,xx2,c2,cc2,c1,cc1,t1,t2,t3,t4,t5,t6,t7,t8)
00085     ADD2(a11.d,aa11.d,c1,cc1,c2,cc2,t1,t2)
00086     MUL2(x2,xx2,c2,cc2,c1,cc1,t1,t2,t3,t4,t5,t6,t7,t8)
00087     ADD2(a9.d ,aa9.d ,c1,cc1,c2,cc2,t1,t2)
00088     MUL2(x2,xx2,c2,cc2,c1,cc1,t1,t2,t3,t4,t5,t6,t7,t8)
00089     ADD2(a7.d ,aa7.d ,c1,cc1,c2,cc2,t1,t2)
00090     MUL2(x2,xx2,c2,cc2,c1,cc1,t1,t2,t3,t4,t5,t6,t7,t8)
00091     ADD2(a5.d ,aa5.d ,c1,cc1,c2,cc2,t1,t2)
00092     MUL2(x2,xx2,c2,cc2,c1,cc1,t1,t2,t3,t4,t5,t6,t7,t8)
00093     ADD2(a3.d ,aa3.d ,c1,cc1,c2,cc2,t1,t2)
00094     MUL2(x2,xx2,c2,cc2,c1,cc1,t1,t2,t3,t4,t5,t6,t7,t8)
00095     MUL2(x ,zero.d,c1,cc1,c2,cc2,t1,t2,t3,t4,t5,t6,t7,t8)
00096     ADD2(x    ,zero.d,c2,cc2,c1,cc1,t1,t2)
00097     if ((y=c1+(cc1-u2.d*c1)) == c1+(cc1+u2.d*c1))  return y;
00098     return tanMp(x);
00099   }
00100 
00101   /* (III) The case 0.0608 < abs(x) <= 0.787 */
00102   if (w<=g3.d) {
00103 
00104     /* First stage */
00105     i = ((int) (mfftnhf.d+TWO8*w));
00106     z = w-xfg[i][0].d;  z2 = z*z;   s = (x<ZERO) ? MONE : ONE;
00107     pz = z+z*z2*(e0.d+z2*e1.d);
00108     fi = xfg[i][1].d;   gi = xfg[i][2].d;   t2 = pz*(gi+fi)/(gi-pz);
00109     if ((y=fi+(t2-fi*u3.d))==fi+(t2+fi*u3.d))  return (s*y);
00110     t3 = (t2<ZERO) ? -t2 : t2;
00111     if ((y=fi+(t2-(t4=fi*ua3.d+t3*ub3.d)))==fi+(t2+t4))  return (s*y);
00112 
00113     /* Second stage */
00114     ffi = xfg[i][3].d;
00115     c1 = z2*(a7.d+z2*(a9.d+z2*a11.d));
00116     EMULV(z,z,z2,zz2,t1,t2,t3,t4,t5)
00117     ADD2(a5.d,aa5.d,c1,zero.d,c2,cc2,t1,t2)
00118     MUL2(z2,zz2,c2,cc2,c1,cc1,t1,t2,t3,t4,t5,t6,t7,t8)
00119     ADD2(a3.d,aa3.d,c1,cc1,c2,cc2,t1,t2)
00120     MUL2(z2,zz2,c2,cc2,c1,cc1,t1,t2,t3,t4,t5,t6,t7,t8)
00121     MUL2(z ,zero.d,c1,cc1,c2,cc2,t1,t2,t3,t4,t5,t6,t7,t8)
00122     ADD2(z ,zero.d,c2,cc2,c1,cc1,t1,t2)
00123 
00124     ADD2(fi ,ffi,c1,cc1,c2,cc2,t1,t2)
00125     MUL2(fi ,ffi,c1,cc1,c3,cc3,t1,t2,t3,t4,t5,t6,t7,t8)
00126     SUB2(one.d,zero.d,c3,cc3,c1,cc1,t1,t2)
00127     DIV2(c2,cc2,c1,cc1,c3,cc3,t1,t2,t3,t4,t5,t6,t7,t8,t9,t10)
00128 
00129     if ((y=c3+(cc3-u4.d*c3))==c3+(cc3+u4.d*c3))  return (s*y);
00130     return tanMp(x);
00131   }
00132 
00133   /* (---) The case 0.787 < abs(x) <= 25 */
00134   if (w<=g4.d) {
00135     /* Range reduction by algorithm i */
00136     t = (x*hpinv.d + toint.d);
00137     xn = t - toint.d;
00138     v.d = t;
00139     t1 = (x - xn*mp1.d) - xn*mp2.d;
00140     n =v.i[LOW_HALF] & 0x00000001;
00141     da = xn*mp3.d;
00142     a=t1-da;
00143     da = (t1-a)-da;
00144     if (a<ZERO)  {ya=-a;  yya=-da;  sy=MONE;}
00145     else         {ya= a;  yya= da;  sy= ONE;}
00146 
00147     /* (IV),(V) The case 0.787 < abs(x) <= 25,    abs(y) <= 1e-7 */
00148     if (ya<=gy1.d)  return tanMp(x);
00149 
00150     /* (VI) The case 0.787 < abs(x) <= 25,    1e-7 < abs(y) <= 0.0608 */
00151     if (ya<=gy2.d) {
00152       a2 = a*a;
00153       t2 = da+a*a2*(d3.d+a2*(d5.d+a2*(d7.d+a2*(d9.d+a2*d11.d))));
00154       if (n) {
00155         /* First stage -cot */
00156         EADD(a,t2,b,db)
00157         DIV2(one.d,zero.d,b,db,c,dc,t1,t2,t3,t4,t5,t6,t7,t8,t9,t10)
00158         if ((y=c+(dc-u6.d*c))==c+(dc+u6.d*c))  return (-y); }
00159       else {
00160         /* First stage tan */
00161         if ((y=a+(t2-u5.d*a))==a+(t2+u5.d*a))  return y; }
00162       /* Second stage */
00163       /* Range reduction by algorithm ii */
00164       t = (x*hpinv.d + toint.d);
00165       xn = t - toint.d;
00166       v.d = t;
00167       t1 = (x - xn*mp1.d) - xn*mp2.d;
00168       n =v.i[LOW_HALF] & 0x00000001;
00169       da = xn*pp3.d;
00170       t=t1-da;
00171       da = (t1-t)-da;
00172       t1 = xn*pp4.d;
00173       a = t - t1;
00174       da = ((t-a)-t1)+da;
00175 
00176       /* Second stage */
00177       EADD(a,da,t1,t2)   a=t1;  da=t2;
00178       MUL2(a,da,a,da,x2,xx2,t1,t2,t3,t4,t5,t6,t7,t8)
00179       c1 = x2*(a15.d+x2*(a17.d+x2*(a19.d+x2*(a21.d+x2*(a23.d+x2*(a25.d+
00180            x2*a27.d))))));
00181       ADD2(a13.d,aa13.d,c1,zero.d,c2,cc2,t1,t2)
00182       MUL2(x2,xx2,c2,cc2,c1,cc1,t1,t2,t3,t4,t5,t6,t7,t8)
00183       ADD2(a11.d,aa11.d,c1,cc1,c2,cc2,t1,t2)
00184       MUL2(x2,xx2,c2,cc2,c1,cc1,t1,t2,t3,t4,t5,t6,t7,t8)
00185       ADD2(a9.d ,aa9.d ,c1,cc1,c2,cc2,t1,t2)
00186       MUL2(x2,xx2,c2,cc2,c1,cc1,t1,t2,t3,t4,t5,t6,t7,t8)
00187       ADD2(a7.d ,aa7.d ,c1,cc1,c2,cc2,t1,t2)
00188       MUL2(x2,xx2,c2,cc2,c1,cc1,t1,t2,t3,t4,t5,t6,t7,t8)
00189       ADD2(a5.d ,aa5.d ,c1,cc1,c2,cc2,t1,t2)
00190       MUL2(x2,xx2,c2,cc2,c1,cc1,t1,t2,t3,t4,t5,t6,t7,t8)
00191       ADD2(a3.d ,aa3.d ,c1,cc1,c2,cc2,t1,t2)
00192       MUL2(x2,xx2,c2,cc2,c1,cc1,t1,t2,t3,t4,t5,t6,t7,t8)
00193       MUL2(a ,da ,c1,cc1,c2,cc2,t1,t2,t3,t4,t5,t6,t7,t8)
00194       ADD2(a  ,da  ,c2,cc2,c1,cc1,t1,t2)
00195 
00196       if (n) {
00197         /* Second stage -cot */
00198         DIV2(one.d,zero.d,c1,cc1,c2,cc2,t1,t2,t3,t4,t5,t6,t7,t8,t9,t10)
00199         if ((y=c2+(cc2-u8.d*c2)) == c2+(cc2+u8.d*c2))  return (-y); }
00200       else {
00201         /* Second stage tan */
00202         if ((y=c1+(cc1-u7.d*c1)) == c1+(cc1+u7.d*c1))  return y; }
00203       return tanMp(x);
00204     }
00205 
00206     /* (VII) The case 0.787 < abs(x) <= 25,    0.0608 < abs(y) <= 0.787 */
00207 
00208     /* First stage */
00209     i = ((int) (mfftnhf.d+TWO8*ya));
00210     z = (z0=(ya-xfg[i][0].d))+yya;  z2 = z*z;
00211     pz = z+z*z2*(e0.d+z2*e1.d);
00212     fi = xfg[i][1].d;   gi = xfg[i][2].d;
00213 
00214     if (n) {
00215       /* -cot */
00216       t2 = pz*(fi+gi)/(fi+pz);
00217       if ((y=gi-(t2-gi*u10.d))==gi-(t2+gi*u10.d))  return (-sy*y);
00218       t3 = (t2<ZERO) ? -t2 : t2;
00219       if ((y=gi-(t2-(t4=gi*ua10.d+t3*ub10.d)))==gi-(t2+t4))  return (-sy*y); }
00220     else   {
00221       /* tan */
00222       t2 = pz*(gi+fi)/(gi-pz);
00223       if ((y=fi+(t2-fi*u9.d))==fi+(t2+fi*u9.d))  return (sy*y);
00224       t3 = (t2<ZERO) ? -t2 : t2;
00225       if ((y=fi+(t2-(t4=fi*ua9.d+t3*ub9.d)))==fi+(t2+t4))  return (sy*y); }
00226 
00227     /* Second stage */
00228     ffi = xfg[i][3].d;
00229     EADD(z0,yya,z,zz)
00230     MUL2(z,zz,z,zz,z2,zz2,t1,t2,t3,t4,t5,t6,t7,t8)
00231     c1 = z2*(a7.d+z2*(a9.d+z2*a11.d));
00232     ADD2(a5.d,aa5.d,c1,zero.d,c2,cc2,t1,t2)
00233     MUL2(z2,zz2,c2,cc2,c1,cc1,t1,t2,t3,t4,t5,t6,t7,t8)
00234     ADD2(a3.d,aa3.d,c1,cc1,c2,cc2,t1,t2)
00235     MUL2(z2,zz2,c2,cc2,c1,cc1,t1,t2,t3,t4,t5,t6,t7,t8)
00236     MUL2(z ,zz ,c1,cc1,c2,cc2,t1,t2,t3,t4,t5,t6,t7,t8)
00237     ADD2(z ,zz ,c2,cc2,c1,cc1,t1,t2)
00238 
00239     ADD2(fi ,ffi,c1,cc1,c2,cc2,t1,t2)
00240     MUL2(fi ,ffi,c1,cc1,c3,cc3,t1,t2,t3,t4,t5,t6,t7,t8)
00241     SUB2(one.d,zero.d,c3,cc3,c1,cc1,t1,t2)
00242 
00243     if (n) {
00244       /* -cot */
00245       DIV2(c1,cc1,c2,cc2,c3,cc3,t1,t2,t3,t4,t5,t6,t7,t8,t9,t10)
00246       if ((y=c3+(cc3-u12.d*c3))==c3+(cc3+u12.d*c3))  return (-sy*y); }
00247     else {
00248       /* tan */
00249       DIV2(c2,cc2,c1,cc1,c3,cc3,t1,t2,t3,t4,t5,t6,t7,t8,t9,t10)
00250       if ((y=c3+(cc3-u11.d*c3))==c3+(cc3+u11.d*c3))  return (sy*y); }
00251 
00252     return tanMp(x);
00253   }
00254 
00255   /* (---) The case 25 < abs(x) <= 1e8 */
00256   if (w<=g5.d) {
00257     /* Range reduction by algorithm ii */
00258     t = (x*hpinv.d + toint.d);
00259     xn = t - toint.d;
00260     v.d = t;
00261     t1 = (x - xn*mp1.d) - xn*mp2.d;
00262     n =v.i[LOW_HALF] & 0x00000001;
00263     da = xn*pp3.d;
00264     t=t1-da;
00265     da = (t1-t)-da;
00266     t1 = xn*pp4.d;
00267     a = t - t1;
00268     da = ((t-a)-t1)+da;
00269     EADD(a,da,t1,t2)   a=t1;  da=t2;
00270     if (a<ZERO)  {ya=-a;  yya=-da;  sy=MONE;}
00271     else         {ya= a;  yya= da;  sy= ONE;}
00272 
00273     /* (+++) The case 25 < abs(x) <= 1e8,    abs(y) <= 1e-7 */
00274     if (ya<=gy1.d)  return tanMp(x);
00275 
00276     /* (VIII) The case 25 < abs(x) <= 1e8,    1e-7 < abs(y) <= 0.0608 */
00277     if (ya<=gy2.d) {
00278       a2 = a*a;
00279       t2 = da+a*a2*(d3.d+a2*(d5.d+a2*(d7.d+a2*(d9.d+a2*d11.d))));
00280       if (n) {
00281         /* First stage -cot */
00282         EADD(a,t2,b,db)
00283         DIV2(one.d,zero.d,b,db,c,dc,t1,t2,t3,t4,t5,t6,t7,t8,t9,t10)
00284         if ((y=c+(dc-u14.d*c))==c+(dc+u14.d*c))  return (-y); }
00285       else {
00286         /* First stage tan */
00287         if ((y=a+(t2-u13.d*a))==a+(t2+u13.d*a))  return y; }
00288 
00289       /* Second stage */
00290       MUL2(a,da,a,da,x2,xx2,t1,t2,t3,t4,t5,t6,t7,t8)
00291       c1 = x2*(a15.d+x2*(a17.d+x2*(a19.d+x2*(a21.d+x2*(a23.d+x2*(a25.d+
00292            x2*a27.d))))));
00293       ADD2(a13.d,aa13.d,c1,zero.d,c2,cc2,t1,t2)
00294       MUL2(x2,xx2,c2,cc2,c1,cc1,t1,t2,t3,t4,t5,t6,t7,t8)
00295       ADD2(a11.d,aa11.d,c1,cc1,c2,cc2,t1,t2)
00296       MUL2(x2,xx2,c2,cc2,c1,cc1,t1,t2,t3,t4,t5,t6,t7,t8)
00297       ADD2(a9.d ,aa9.d ,c1,cc1,c2,cc2,t1,t2)
00298       MUL2(x2,xx2,c2,cc2,c1,cc1,t1,t2,t3,t4,t5,t6,t7,t8)
00299       ADD2(a7.d ,aa7.d ,c1,cc1,c2,cc2,t1,t2)
00300       MUL2(x2,xx2,c2,cc2,c1,cc1,t1,t2,t3,t4,t5,t6,t7,t8)
00301       ADD2(a5.d ,aa5.d ,c1,cc1,c2,cc2,t1,t2)
00302       MUL2(x2,xx2,c2,cc2,c1,cc1,t1,t2,t3,t4,t5,t6,t7,t8)
00303       ADD2(a3.d ,aa3.d ,c1,cc1,c2,cc2,t1,t2)
00304       MUL2(x2,xx2,c2,cc2,c1,cc1,t1,t2,t3,t4,t5,t6,t7,t8)
00305       MUL2(a ,da ,c1,cc1,c2,cc2,t1,t2,t3,t4,t5,t6,t7,t8)
00306       ADD2(a  ,da  ,c2,cc2,c1,cc1,t1,t2)
00307 
00308       if (n) {
00309         /* Second stage -cot */
00310         DIV2(one.d,zero.d,c1,cc1,c2,cc2,t1,t2,t3,t4,t5,t6,t7,t8,t9,t10)
00311         if ((y=c2+(cc2-u16.d*c2)) == c2+(cc2+u16.d*c2))  return (-y); }
00312       else {
00313         /* Second stage tan */
00314         if ((y=c1+(cc1-u15.d*c1)) == c1+(cc1+u15.d*c1))  return (y); }
00315       return tanMp(x);
00316     }
00317 
00318     /* (IX) The case 25 < abs(x) <= 1e8,    0.0608 < abs(y) <= 0.787 */
00319     /* First stage */
00320     i = ((int) (mfftnhf.d+TWO8*ya));
00321     z = (z0=(ya-xfg[i][0].d))+yya;  z2 = z*z;
00322     pz = z+z*z2*(e0.d+z2*e1.d);
00323     fi = xfg[i][1].d;   gi = xfg[i][2].d;
00324 
00325     if (n) {
00326       /* -cot */
00327       t2 = pz*(fi+gi)/(fi+pz);
00328       if ((y=gi-(t2-gi*u18.d))==gi-(t2+gi*u18.d))  return (-sy*y);
00329       t3 = (t2<ZERO) ? -t2 : t2;
00330       if ((y=gi-(t2-(t4=gi*ua18.d+t3*ub18.d)))==gi-(t2+t4))  return (-sy*y); }
00331     else   {
00332       /* tan */
00333       t2 = pz*(gi+fi)/(gi-pz);
00334       if ((y=fi+(t2-fi*u17.d))==fi+(t2+fi*u17.d))  return (sy*y);
00335       t3 = (t2<ZERO) ? -t2 : t2;
00336       if ((y=fi+(t2-(t4=fi*ua17.d+t3*ub17.d)))==fi+(t2+t4))  return (sy*y); }
00337 
00338     /* Second stage */
00339     ffi = xfg[i][3].d;
00340     EADD(z0,yya,z,zz)
00341     MUL2(z,zz,z,zz,z2,zz2,t1,t2,t3,t4,t5,t6,t7,t8)
00342     c1 = z2*(a7.d+z2*(a9.d+z2*a11.d));
00343     ADD2(a5.d,aa5.d,c1,zero.d,c2,cc2,t1,t2)
00344     MUL2(z2,zz2,c2,cc2,c1,cc1,t1,t2,t3,t4,t5,t6,t7,t8)
00345     ADD2(a3.d,aa3.d,c1,cc1,c2,cc2,t1,t2)
00346     MUL2(z2,zz2,c2,cc2,c1,cc1,t1,t2,t3,t4,t5,t6,t7,t8)
00347     MUL2(z ,zz ,c1,cc1,c2,cc2,t1,t2,t3,t4,t5,t6,t7,t8)
00348     ADD2(z ,zz ,c2,cc2,c1,cc1,t1,t2)
00349 
00350     ADD2(fi ,ffi,c1,cc1,c2,cc2,t1,t2)
00351     MUL2(fi ,ffi,c1,cc1,c3,cc3,t1,t2,t3,t4,t5,t6,t7,t8)
00352     SUB2(one.d,zero.d,c3,cc3,c1,cc1,t1,t2)
00353 
00354     if (n) {
00355       /* -cot */
00356       DIV2(c1,cc1,c2,cc2,c3,cc3,t1,t2,t3,t4,t5,t6,t7,t8,t9,t10)
00357       if ((y=c3+(cc3-u20.d*c3))==c3+(cc3+u20.d*c3))  return (-sy*y); }
00358     else {
00359       /* tan */
00360       DIV2(c2,cc2,c1,cc1,c3,cc3,t1,t2,t3,t4,t5,t6,t7,t8,t9,t10)
00361       if ((y=c3+(cc3-u19.d*c3))==c3+(cc3+u19.d*c3))  return (sy*y); }
00362     return tanMp(x);
00363   }
00364 
00365   /* (---) The case 1e8 < abs(x) < 2**1024 */
00366   /* Range reduction by algorithm iii */
00367   n = (__branred(x,&a,&da)) & 0x00000001;
00368   EADD(a,da,t1,t2)   a=t1;  da=t2;
00369   if (a<ZERO)  {ya=-a;  yya=-da;  sy=MONE;}
00370   else         {ya= a;  yya= da;  sy= ONE;}
00371 
00372   /* (+++) The case 1e8 < abs(x) < 2**1024,    abs(y) <= 1e-7 */
00373   if (ya<=gy1.d)  return tanMp(x);
00374 
00375   /* (X) The case 1e8 < abs(x) < 2**1024,    1e-7 < abs(y) <= 0.0608 */
00376   if (ya<=gy2.d) {
00377     a2 = a*a;
00378     t2 = da+a*a2*(d3.d+a2*(d5.d+a2*(d7.d+a2*(d9.d+a2*d11.d))));
00379     if (n) {
00380       /* First stage -cot */
00381       EADD(a,t2,b,db)
00382       DIV2(one.d,zero.d,b,db,c,dc,t1,t2,t3,t4,t5,t6,t7,t8,t9,t10)
00383       if ((y=c+(dc-u22.d*c))==c+(dc+u22.d*c))  return (-y); }
00384     else {
00385       /* First stage tan */
00386       if ((y=a+(t2-u21.d*a))==a+(t2+u21.d*a))  return y; }
00387 
00388     /* Second stage */
00389     /* Reduction by algorithm iv */
00390     p=10;    n = (__mpranred(x,&mpa,p)) & 0x00000001;
00391     __mp_dbl(&mpa,&a,p);        __dbl_mp(a,&mpt1,p);
00392     __sub(&mpa,&mpt1,&mpt2,p);  __mp_dbl(&mpt2,&da,p);
00393 
00394     MUL2(a,da,a,da,x2,xx2,t1,t2,t3,t4,t5,t6,t7,t8)
00395     c1 = x2*(a15.d+x2*(a17.d+x2*(a19.d+x2*(a21.d+x2*(a23.d+x2*(a25.d+
00396          x2*a27.d))))));
00397     ADD2(a13.d,aa13.d,c1,zero.d,c2,cc2,t1,t2)
00398     MUL2(x2,xx2,c2,cc2,c1,cc1,t1,t2,t3,t4,t5,t6,t7,t8)
00399     ADD2(a11.d,aa11.d,c1,cc1,c2,cc2,t1,t2)
00400     MUL2(x2,xx2,c2,cc2,c1,cc1,t1,t2,t3,t4,t5,t6,t7,t8)
00401     ADD2(a9.d ,aa9.d ,c1,cc1,c2,cc2,t1,t2)
00402     MUL2(x2,xx2,c2,cc2,c1,cc1,t1,t2,t3,t4,t5,t6,t7,t8)
00403     ADD2(a7.d ,aa7.d ,c1,cc1,c2,cc2,t1,t2)
00404     MUL2(x2,xx2,c2,cc2,c1,cc1,t1,t2,t3,t4,t5,t6,t7,t8)
00405     ADD2(a5.d ,aa5.d ,c1,cc1,c2,cc2,t1,t2)
00406     MUL2(x2,xx2,c2,cc2,c1,cc1,t1,t2,t3,t4,t5,t6,t7,t8)
00407     ADD2(a3.d ,aa3.d ,c1,cc1,c2,cc2,t1,t2)
00408     MUL2(x2,xx2,c2,cc2,c1,cc1,t1,t2,t3,t4,t5,t6,t7,t8)
00409     MUL2(a ,da ,c1,cc1,c2,cc2,t1,t2,t3,t4,t5,t6,t7,t8)
00410     ADD2(a    ,da    ,c2,cc2,c1,cc1,t1,t2)
00411 
00412     if (n) {
00413       /* Second stage -cot */
00414       DIV2(one.d,zero.d,c1,cc1,c2,cc2,t1,t2,t3,t4,t5,t6,t7,t8,t9,t10)
00415       if ((y=c2+(cc2-u24.d*c2)) == c2+(cc2+u24.d*c2))  return (-y); }
00416     else {
00417       /* Second stage tan */
00418       if ((y=c1+(cc1-u23.d*c1)) == c1+(cc1+u23.d*c1))  return y; }
00419     return tanMp(x);
00420   }
00421 
00422   /* (XI) The case 1e8 < abs(x) < 2**1024,    0.0608 < abs(y) <= 0.787 */
00423   /* First stage */
00424   i = ((int) (mfftnhf.d+TWO8*ya));
00425   z = (z0=(ya-xfg[i][0].d))+yya;  z2 = z*z;
00426   pz = z+z*z2*(e0.d+z2*e1.d);
00427   fi = xfg[i][1].d;   gi = xfg[i][2].d;
00428 
00429   if (n) {
00430     /* -cot */
00431     t2 = pz*(fi+gi)/(fi+pz);
00432     if ((y=gi-(t2-gi*u26.d))==gi-(t2+gi*u26.d))  return (-sy*y);
00433     t3 = (t2<ZERO) ? -t2 : t2;
00434     if ((y=gi-(t2-(t4=gi*ua26.d+t3*ub26.d)))==gi-(t2+t4))  return (-sy*y); }
00435   else   {
00436     /* tan */
00437     t2 = pz*(gi+fi)/(gi-pz);
00438     if ((y=fi+(t2-fi*u25.d))==fi+(t2+fi*u25.d))  return (sy*y);
00439     t3 = (t2<ZERO) ? -t2 : t2;
00440     if ((y=fi+(t2-(t4=fi*ua25.d+t3*ub25.d)))==fi+(t2+t4))  return (sy*y); }
00441 
00442   /* Second stage */
00443   ffi = xfg[i][3].d;
00444   EADD(z0,yya,z,zz)
00445   MUL2(z,zz,z,zz,z2,zz2,t1,t2,t3,t4,t5,t6,t7,t8)
00446   c1 = z2*(a7.d+z2*(a9.d+z2*a11.d));
00447   ADD2(a5.d,aa5.d,c1,zero.d,c2,cc2,t1,t2)
00448   MUL2(z2,zz2,c2,cc2,c1,cc1,t1,t2,t3,t4,t5,t6,t7,t8)
00449   ADD2(a3.d,aa3.d,c1,cc1,c2,cc2,t1,t2)
00450   MUL2(z2,zz2,c2,cc2,c1,cc1,t1,t2,t3,t4,t5,t6,t7,t8)
00451   MUL2(z ,zz ,c1,cc1,c2,cc2,t1,t2,t3,t4,t5,t6,t7,t8)
00452   ADD2(z ,zz ,c2,cc2,c1,cc1,t1,t2)
00453 
00454   ADD2(fi ,ffi,c1,cc1,c2,cc2,t1,t2)
00455   MUL2(fi ,ffi,c1,cc1,c3,cc3,t1,t2,t3,t4,t5,t6,t7,t8)
00456   SUB2(one.d,zero.d,c3,cc3,c1,cc1,t1,t2)
00457 
00458   if (n) {
00459     /* -cot */
00460     DIV2(c1,cc1,c2,cc2,c3,cc3,t1,t2,t3,t4,t5,t6,t7,t8,t9,t10)
00461     if ((y=c3+(cc3-u28.d*c3))==c3+(cc3+u28.d*c3))  return (-sy*y); }
00462   else {
00463     /* tan */
00464     DIV2(c2,cc2,c1,cc1,c3,cc3,t1,t2,t3,t4,t5,t6,t7,t8,t9,t10)
00465     if ((y=c3+(cc3-u27.d*c3))==c3+(cc3+u27.d*c3))  return (sy*y); }
00466   return tanMp(x);
00467 }
00468 
00469 
00470 /* multiple precision stage                                              */
00471 /* Convert x to multi precision number,compute tan(x) by mptan() routine */
00472 /* and converts result back to double                                    */
00473 static double tanMp(double x)
00474 {
00475   int p;
00476   double y;
00477   mp_no mpy;
00478   p=32;
00479   __mptan(x, &mpy, p);
00480   __mp_dbl(&mpy,&y,p);
00481   return y;
00482 }
00483 
00484 #ifdef NO_LONG_DOUBLE
00485 weak_alias (tan, tanl)
00486 #endif