Open64: osprey/libm/mips/tan.c Source File

00001 /*
00002  * Copyright 2003, 2004, 2005, 2006 PathScale, Inc.  All Rights Reserved.
00003  */
00004 
00005 
00006 /*
00007 
00008   Copyright (C) 2000, 2001 Silicon Graphics, Inc.  All Rights Reserved.
00009 
00010   This program is free software; you can redistribute it and/or modify it
00011   under the terms of version 2.1 of the GNU Lesser General Public License 
00012   as published by the Free Software Foundation.
00013 
00014   This program is distributed in the hope that it would be useful, but
00015   WITHOUT ANY WARRANTY; without even the implied warranty of
00016   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  
00017 
00018   Further, this software is distributed without any warranty that it is
00019   free of the rightful claim of any third person regarding infringement 
00020   or the like.  Any license provided herein, whether implied or 
00021   otherwise, applies only to this software file.  Patent licenses, if
00022   any, provided herein do not apply to combinations of this program with 
00023   other software, or any other product whatsoever.  
00024 
00025   You should have received a copy of the GNU Lesser General Public 
00026   License along with this program; if not, write the Free Software 
00027   Foundation, Inc., 59 Temple Place - Suite 330, Boston MA 02111-1307, 
00028   USA.
00029 
00030   Contact information:  Silicon Graphics, Inc., 1600 Amphitheatre Pky,
00031   Mountain View, CA 94043, or:
00032 
00033   http://www.sgi.com
00034 
00035   For further information regarding this notice, see:
00036 
00037   http://oss.sgi.com/projects/GenInfo/NoticeExplan
00038 
00039 */
00040 
00041 /* ====================================================================
00042  * ====================================================================
00043  *
00044  * Module: tan.c
00045  * $Revision: 1.5 $
00046  * $Date: 04/12/21 14:58:23-08:00 $
00047  * $Author: bos@eng-25.internal.keyresearch.com $
00048  * $Source: /home/bos/bk/kpro64-pending/libm/mips/SCCS/s.tan.c $
00049  *
00050  * Revision history:
00051  *  09-Jun-93 - Original Version
00052  *
00053  * Description: source code for tan function
00054  *
00055  * ====================================================================
00056  * ====================================================================
00057  */
00058 
00059 static char *rcs_id = "$Source: /home/bos/bk/kpro64-pending/libm/mips/SCCS/s.tan.c $ $Revision: 1.5 $";
00060 
00061 #ifdef _CALL_MATHERR
00062 #include <stdio.h>
00063 #include <math.h>
00064 #include <errno.h>
00065 #endif
00066 
00067 #include "libm.h"
00068 
00069 #if defined(mips) && !defined(__GNUC__)
00070 extern  double  tan(double);
00071 
00072 #pragma weak tan = __tan
00073 #endif
00074 
00075 #if defined(BUILD_OS_DARWIN) /* Mach-O doesn't support aliases */
00076 extern double __tan(double);
00077 #pragma weak tan
00078 double tan( double x ) {
00079   return __tan( x );
00080 }
00081 #elif defined(__GNUC__)
00082 extern  double  __tan(double);
00083 
00084 double    tan() __attribute__ ((weak, alias ("__tan")));
00085 
00086 #endif
00087 
00088 static const du Qnan =
00089 {D(QNANHI, QNANLO)};
00090 
00091 static const du Inf =
00092 {D(0x7ff00000, 0x00000000)};
00093 
00094 static const du half =
00095 {D(0x3fe00000, 0x00000000)};
00096 
00097 static const du one =
00098 {D(0x3ff00000, 0x00000000)};
00099 
00100 static const du twop19xpi =
00101 {D(0x413921fb, 0x54442d18)};
00102 
00103 static const du twopm26 =
00104 {D(0x3e500000, 0x00000000)};
00105 
00106 static const du rpiby2 =
00107 {D(0x3fe45f30, 0x6dc9c883)};
00108 
00109 static const du piby2hi =
00110 {D(0x3ff921fb, 0x54400000)};
00111 
00112 static const du piby2lo =
00113 {D(0x3dd0b461, 0x1a600000)};
00114 
00115 static const du piby2tiny =
00116 {D(0x3ba3198a, 0x2e037073)};
00117 
00118 static const du ph =
00119 {D(0x3ff921fb, 0x50000000)};
00120 
00121 static const du pl =
00122 {D(0x3e5110b4, 0x60000000)};
00123 
00124 static const du pt =
00125 {D(0x3c91a626, 0x30000000)};
00126 
00127 static const du pe =
00128 {D(0x3ae8a2e0, 0x30000000)};
00129 
00130 static const du pe2 =
00131 {D(0x394c1cd1, 0x29024e09)};
00132 
00133 static const du Ph =
00134 {D(0x3ff921fb, 0x54000000)};
00135 
00136 static const du Pl =
00137 {D(0x3e110b46, 0x10000000)};
00138 
00139 static const du Pt =
00140 {D(0x3c5a6263, 0x3145c06e)};
00141 
00142 /* coefficients for rational approximation of tan on +/- pi/4 */
00143 
00144 static const du p[] =
00145 {
00146 {D(0x3ff00000, 0x00000000)},
00147 {D(0xbfc06b8f, 0x5f225706)},
00148 {D(0x3f66fc34, 0x2943627f)},
00149 {D(0xbedf6255, 0x88fc315e)},
00150 };
00151 
00152 static const du q[] =
00153 {
00154 {D(0x3ff00000, 0x00000000)},
00155 {D(0xbfdd8b1d, 0x04e680b5)},
00156 {D(0x3f97e798, 0xc16cfadc)},
00157 {D(0xbf2b51d1, 0xd7f65e57)},
00158 };
00159 
00160 #ifdef _TABLE_BASED_REDUCTION
00161 
00162 /* tables used to do argument reduction for args between +/- 16 radians;
00163    the sum of the high and low values of the kth entry is (k - 10)*pi/2
00164 */
00165 
00166 static const du tblh[] =
00167 {
00168 {D(0xc02f6a7a, 0x2955385e)},
00169 {D(0xc02c463a, 0xbeccb2bb)},
00170 {D(0xc02921fb, 0x54442d18)},
00171 {D(0xc025fdbb, 0xe9bba775)},
00172 {D(0xc022d97c, 0x7f3321d2)},
00173 {D(0xc01f6a7a, 0x2955385e)},
00174 {D(0xc01921fb, 0x54442d18)},
00175 {D(0xc012d97c, 0x7f3321d2)},
00176 {D(0xc00921fb, 0x54442d18)},
00177 {D(0xbff921fb, 0x54442d18)},
00178 {D(0x00000000, 0x00000000)},
00179 {D(0x3ff921fb, 0x54442d18)},
00180 {D(0x400921fb, 0x54442d18)},
00181 {D(0x4012d97c, 0x7f3321d2)},
00182 {D(0x401921fb, 0x54442d18)},
00183 {D(0x401f6a7a, 0x2955385e)},
00184 {D(0x4022d97c, 0x7f3321d2)},
00185 {D(0x4025fdbb, 0xe9bba775)},
00186 {D(0x402921fb, 0x54442d18)},
00187 {D(0x402c463a, 0xbeccb2bb)},
00188 {D(0x402f6a7a, 0x2955385e)},
00189 };
00190 
00191 static const du tbll[] =
00192 {
00193 {D(0xbcc60faf, 0xbfd97309)},
00194 {D(0xbcc3daea, 0xf976e788)},
00195 {D(0xbcc1a626, 0x33145c07)},
00196 {D(0xbcbee2c2, 0xd963a10c)},
00197 {D(0xbcba7939, 0x4c9e8a0a)},
00198 {D(0xbcb60faf, 0xbfd97309)},
00199 {D(0xbcb1a626, 0x33145c07)},
00200 {D(0xbcaa7939, 0x4c9e8a0a)},
00201 {D(0xbca1a626, 0x33145c07)},
00202 {D(0xbc91a626, 0x33145c07)},
00203 {D(0x00000000, 0x00000000)},
00204 {D(0x3c91a626, 0x33145c07)},
00205 {D(0x3ca1a626, 0x33145c07)},
00206 {D(0x3caa7939, 0x4c9e8a0a)},
00207 {D(0x3cb1a626, 0x33145c07)},
00208 {D(0x3cb60faf, 0xbfd97309)},
00209 {D(0x3cba7939, 0x4c9e8a0a)},
00210 {D(0x3cbee2c2, 0xd963a10c)},
00211 {D(0x3cc1a626, 0x33145c07)},
00212 {D(0x3cc3daea, 0xf976e788)},
00213 {D(0x3cc60faf, 0xbfd97309)},
00214 };
00215 #endif
00216 
00217 
00218 /* ====================================================================
00219  *
00220  * FunctionName   tan
00221  *
00222  * Description    computes tangent of arg
00223  *
00224  * ====================================================================
00225  */
00226 
00227 double
00228 __tan( x )
00229 double  x;
00230 {
00231 #ifdef _32BIT_MACHINE
00232 int ix, xpt, m, l;
00233 #else
00234 long long ix, xpt, m, l;
00235 #endif
00236 
00237 int n;
00238 double  xsq;
00239 double  num, denom;
00240 double  poly1, poly2;
00241 double  result;
00242 double  z, dn;
00243 double  absx, dn1, dn2;
00244 double  s, ss;
00245 double  t, w, ww;
00246 #ifdef _CALL_MATHERR
00247 struct exception  exstruct;
00248 #endif
00249 
00250   /* extract exponent of x and 1 bit of mantissa for some quick screening */
00251 
00252 #ifdef _32BIT_MACHINE
00253 
00254   DBLHI2INT(x, ix); /* copy MSW of x to ix  */
00255 #else
00256   DBL2LL(x, ix);    /* copy x to ix */
00257 #endif
00258   xpt = (ix >> (DMANTWIDTH-1));
00259   xpt &= 0xfff;
00260 
00261   if ( xpt < 0x7fd )
00262   {
00263     /*   |x| < 0.75   */
00264 
00265     if ( fabs(x) > twopm26.d )
00266     {
00267       /*   |x| > 2^(-26)   */
00268 
00269       /* just compute tan(x) */
00270 
00271       xsq = x*x;
00272 
00273       poly1 = ((p[3].d*xsq + p[2].d)*xsq + p[1].d)*(xsq*x) + x;
00274       poly2 = ((q[3].d*xsq + q[2].d)*xsq + q[1].d)*xsq + q[0].d;
00275 
00276       return ( poly1/poly2 );
00277     }
00278 
00279     return ( x );
00280   }
00281 
00282 #ifdef _TABLE_BASED_REDUCTION
00283 
00284   if ( xpt < 0x806 )
00285   {
00286     /*   |x| < 16.0   */
00287 
00288     /*  do a table based argument reduction to +/- pi/4  */
00289 
00290     dn = x*rpiby2.d;
00291     n = ROUND(dn);
00292 
00293     /*  compute z = x - n*pi/2  */
00294 
00295     x = x - tblh[n+10].d;
00296     x = x - tbll[n+10].d;
00297 
00298     xsq = x*x;
00299   
00300     poly1 = ((p[3].d*xsq + p[2].d)*xsq + p[1].d)*(xsq*x) + x;
00301     poly2 = ((q[3].d*xsq + q[2].d)*xsq + q[1].d)*xsq + q[0].d;
00302 
00303     if ( (n&1) == 0 )
00304     {
00305       /* result is poly1/poly2 */
00306 
00307       num = poly1;
00308       denom = poly2;
00309     }
00310 
00311     if ( (n&1) != 0 )
00312     {
00313       /* result is -poly2/poly1 */
00314 
00315       denom = poly1;
00316       num = -poly2;
00317     }
00318 
00319     result = num/denom;
00320 
00321     return ( result );
00322   }
00323 #endif
00324 
00325   if ( xpt < 0x827 )
00326   {
00327     /*  |x| < 1.5*2^20  */
00328 cont:
00329     dn = x*rpiby2.d;
00330     n = ROUND(dn);
00331     dn = n;
00332 
00333     x = x - dn*piby2hi.d;
00334     x = x - dn*piby2lo.d;
00335     x = x - dn*piby2tiny.d; /* x = x - n*pi/2 */
00336 
00337     xsq = x*x;
00338   
00339     poly1 = ((p[3].d*xsq + p[2].d)*xsq + p[1].d)*(xsq*x) + x;
00340     poly2 = ((q[3].d*xsq + q[2].d)*xsq + q[1].d)*xsq + q[0].d;
00341 
00342     if ( (n&1) == 0 )
00343     {
00344       /* result is poly1/poly2 */
00345 
00346       num = poly1;
00347       denom = poly2;
00348     }
00349 
00350     if ( (n&1) != 0 )
00351     {
00352       /* result is -poly2/poly1 */
00353 
00354       denom = poly1;
00355       num = -poly2;
00356     }
00357 
00358     result = num/denom;
00359 
00360     return ( result );
00361   }
00362 
00363   if ( xpt < 0x836 )
00364   {
00365     /*  |x| < 2^28  */
00366 
00367     if ( fabs(x) < twop19xpi.d )
00368       goto cont;
00369 
00370     dn = x*rpiby2.d;
00371     n = ROUND(dn);
00372     dn = n;
00373 
00374     x = x - dn*ph.d;
00375     x = x - dn*pl.d;
00376     x = x - dn*pt.d;
00377     x = x - dn*pe.d;
00378     x = x - dn*pe2.d;
00379 
00380     xsq = x*x;
00381   
00382     poly1 = ((p[3].d*xsq + p[2].d)*xsq + p[1].d)*(xsq*x) + x;
00383     poly2 = ((q[3].d*xsq + q[2].d)*xsq + q[1].d)*xsq + q[0].d;
00384 
00385     if ( (n&1) == 0 )
00386     {
00387       /* result is poly1/poly2 */
00388 
00389       num = poly1;
00390       denom = poly2;
00391     }
00392 
00393     if ( (n&1) != 0 )
00394     {
00395       /* result is -poly2/poly1 */
00396 
00397       denom = poly1;
00398       num = -poly2;
00399     }
00400 
00401     result = num/denom;
00402 
00403     return ( result );
00404   }
00405 
00406   if ( xpt < 0x862 )
00407   {
00408     /*  |x| < 2^50  */
00409 
00410     absx = fabs(x);
00411 
00412     dn = z = absx*rpiby2.d;
00413 
00414     /* round dn to the nearest integer */
00415 
00416 #ifdef _32BIT_MACHINE
00417 
00418     DBLHI2INT(dn, l);
00419     m = (l >> DMANTWIDTH);
00420     m &= 0x7ff;
00421 
00422     /* shift off fractional bits of dn */
00423 
00424     DBLLO2INT(dn, l);
00425 
00426     l >>= (0x433 - m);
00427     n = l;
00428     l <<= (0x433 - m);
00429     INT2DBLLO(l, dn);
00430 #else
00431     DBL2LL(dn, l);
00432     m = (l >> DMANTWIDTH);
00433     m &= 0x7ff;
00434 
00435     /* shift off fractional bits of dn */
00436 
00437     l >>= (0x433 - m);
00438     n = l;
00439     l <<= (0x433 - m);
00440     LL2DBL(l, dn);
00441 #endif
00442     /* adjust dn and n if the fractional part of dn 
00443        was >= 0.5
00444     */
00445 
00446     n &= 3;
00447 
00448     if ( (z - dn) >= half.d )
00449     {
00450       dn += one.d;
00451       n += 1;
00452     }
00453 
00454   /* compute x - dn*Ph - dn*Pl - dn*Pt by dividing dn into
00455      two parts and using the Kahan summation formula
00456   */
00457 
00458     /* split dn into 2 parts */
00459 
00460     dn1 = dn;
00461 
00462 #ifdef _32BIT_MACHINE
00463 
00464     DBLLO2INT(dn1, m);
00465     m >>= 28;
00466     m <<= 28;
00467     INT2DBLLO(m, dn1);
00468 #else
00469     DBL2LL(dn1, m);
00470     m >>= 28;
00471     m <<= 28;
00472     LL2DBL(m, dn1);
00473 #endif
00474     dn2 = dn - dn1;
00475 
00476     z = absx - dn1*Ph.d;  /* this operation is exact */
00477 
00478     t = dn2*Ph.d;
00479     s = z - t;
00480     ss = z - s - t;   /* correction term */
00481 
00482     t = ss - dn1*Pl.d;
00483     w = s + t;
00484     ww = s - w + t;   /* correction term */
00485 
00486     t = ww - dn2*Pl.d;
00487     s = w + t;
00488     ss = w - s + t;   /* correction term */
00489 
00490     t = ss - dn1*Pt.d;
00491     w = s + t;
00492     ww = s - w + t;   /* correction term */
00493 
00494     t = ww - dn2*Pt.d;
00495     z = w + t;    /* z = reduced arg */
00496 
00497     if ( x < 0.0 )
00498     { /* adjust for sign of arg */
00499 
00500       z = -z;
00501       n = -n;
00502     }
00503 
00504     xsq = z*z;
00505   
00506     poly1 = ((p[3].d*xsq + p[2].d)*xsq + p[1].d)*(xsq*z) + z;
00507     poly2 = ((q[3].d*xsq + q[2].d)*xsq + q[1].d)*xsq + q[0].d;
00508 
00509     if ( (n&1) == 0 )
00510     {
00511       /* result is poly1/poly2 */
00512 
00513       num = poly1;
00514       denom = poly2;
00515     }
00516 
00517     if ( (n&1) != 0 )
00518     {
00519       /* result is -poly2/poly1 */
00520 
00521       denom = poly1;
00522       num = -poly2;
00523     }
00524 
00525     result = num/denom;
00526 
00527     return ( result );
00528   }
00529 
00530   if ( (x != x) || (fabs(x) == Inf.d) )
00531   {
00532     /* x is a NaN or +/-inf; return a quiet NaN */
00533 
00534 #ifdef _CALL_MATHERR
00535 
00536     exstruct.type = DOMAIN;
00537     exstruct.name = "tan";
00538     exstruct.arg1 = x;
00539     exstruct.retval = Qnan.d;
00540 
00541     if ( matherr( &exstruct ) == 0 )
00542     {
00543       fprintf(stderr, "domain error in tan\n");
00544       SETERRNO(EDOM);
00545     }
00546 
00547     return ( exstruct.retval );
00548 #else
00549     NAN_SETERRNO(EDOM);
00550 
00551     return ( Qnan.d );
00552 #endif
00553   }
00554 
00555   /* just give up and return 0.0, setting errno = ERANGE */
00556 
00557 #ifdef _CALL_MATHERR
00558 
00559     exstruct.type = TLOSS;
00560     exstruct.name = "tan";
00561     exstruct.arg1 = x;
00562     exstruct.retval = 0.0;
00563 
00564     if ( matherr( &exstruct ) == 0 )
00565     {
00566       fprintf(stderr, "range error in tan (total loss \
00567 of significance)\n");
00568       SETERRNO(ERANGE);
00569     }
00570 
00571     return ( exstruct.retval );
00572 #else
00573     SETERRNO(ERANGE);
00574 
00575     return ( 0.0 );
00576 #endif
00577 }
00578 
00579 #ifdef NO_LONG_DOUBLE
00580 
00581 #if defined(BUILD_OS_DARWIN) /* Mach-O doesn't support aliases */
00582 extern long double __tanl(long double);
00583 long double tanl( long double x ) { 
00584   return ( (long double)__tan((double)x) );
00585 }
00586 #elif defined(__GNUC__)
00587 extern  long double  __tanl(long double);
00588 
00589 long double    tanl() __attribute__ ((weak, alias ("__tanl")));
00590 
00591 #endif
00592 
00593 long double
00594 __tanl( long double x )
00595 { 
00596   return ( (long double)__tan((double)x) );
00597 }
00598 
00599 #endif
00600