lib/libm/tan.c

   1 /*                                                      tan.c
   2  *
   3  *      Circular tangent
   4  *
   5  *
   6  *
   7  * SYNOPSIS:
   8  *
   9  * double x, y, tan();
  10  *
  11  * y = tan( x );
  12  *
  13  *
  14  *
  15  * DESCRIPTION:
  16  *
  17  * Returns the circular tangent of the radian argument x.
  18  *
  19  * Range reduction is modulo pi/4.  A rational function
  20  *       x + x**3 P(x**2)/Q(x**2)
  21  * is employed in the basic interval [0, pi/4].
  22  *
  23  *
  24  *
  25  * ACCURACY:
  26  *
  27  *                      Relative error:
  28  * arithmetic   domain     # trials      peak         rms
  29  *    DEC      +-1.07e9      44000      4.1e-17     1.0e-17
  30  *    IEEE     +-1.07e9      30000      2.9e-16     8.1e-17
  31  *
  32  * ERROR MESSAGES:
  33  *
  34  *   message         condition          value returned
  35  * tan total loss   x > 1.073741824e9     0.0
  36  *
  37  */
  38 \f/*                                                     cot.c
  39  *
  40  *      Circular cotangent
  41  *
  42  *
  43  *
  44  * SYNOPSIS:
  45  *
  46  * double x, y, cot();
  47  *
  48  * y = cot( x );
  49  *
  50  *
  51  *
  52  * DESCRIPTION:
  53  *
  54  * Returns the circular cotangent of the radian argument x.
  55  *
  56  * Range reduction is modulo pi/4.  A rational function
  57  *       x + x**3 P(x**2)/Q(x**2)
  58  * is employed in the basic interval [0, pi/4].
  59  *
  60  *
  61  *
  62  * ACCURACY:
  63  *
  64  *                      Relative error:
  65  * arithmetic   domain     # trials      peak         rms
  66  *    IEEE     +-1.07e9      30000      2.9e-16     8.2e-17
  67  *
  68  *
  69  * ERROR MESSAGES:
  70  *
  71  *   message         condition          value returned
  72  * cot total loss   x > 1.073741824e9       0.0
  73  * cot singularity  x = 0                  INFINITY
  74  *
  75  */
  76 \f
  77 /*
  78 Cephes Math Library Release 2.8:  June, 2000
  79 yright 1984, 1995, 2000 by Stephen L. Moshier
  80 */
  81
  82 #include "mconf.h"
  83
  84 #ifdef UNK
  85 static double P[] = {
  86 -1.30936939181383777646E4,
  87  1.15351664838587416140E6,
  88 -1.79565251976484877988E7
  89 };
  90 static double Q[] = {
  91 /* 1.00000000000000000000E0,*/
  92  1.36812963470692954678E4,
  93 -1.32089234440210967447E6,
  94  2.50083801823357915839E7,
  95 -5.38695755929454629881E7
  96 };
  97 static double DP1 = 7.853981554508209228515625E-1;
  98 static double DP2 = 7.94662735614792836714E-9;
  99 static double DP3 = 3.06161699786838294307E-17;
 100 static double lossth = 1.073741824e9;
 101 #endif
 102
 103 #ifdef DEC
 104 static unsigned short P[] = {
 105 0143514,0113306,0111171,0174674,
 106 0045214,0147545,0027744,0167346,
 107 0146210,0177526,0114514,0105660
 108 };
 109 static unsigned short Q[] = {
 110 /*0040200,0000000,0000000,0000000,*/
 111 0043525,0142457,0072633,0025617,
 112 0145241,0036742,0140525,0162256,
 113 0046276,0146176,0013526,0143573,
 114 0146515,0077401,0162762,0150607
 115 };
 116 /*  7.853981629014015197753906250000E-1 */
 117 static unsigned short P1[] = {0040111,0007732,0120000,0000000,};
 118 /*  4.960467869796758577649598009884E-10 */
 119 static unsigned short P2[] = {0030410,0055060,0100000,0000000,};
 120 /*  2.860594363054915898381331279295E-18 */
 121 static unsigned short P3[] = {0021523,0011431,0105056,0001560,};
 122 #define DP1 *(double *)P1
 123 #define DP2 *(double *)P2
 124 #define DP3 *(double *)P3
 125 static double lossth = 1.073741824e9;
 126 #endif
 127
 128 #ifdef IBMPC
 129 static unsigned short P[] = {
 130 0x3f38,0xd24f,0x92d8,0xc0c9,
 131 0x9ddd,0xa5fc,0x99ec,0x4131,
 132 0x9176,0xd329,0x1fea,0xc171
 133 };
 134 static unsigned short Q[] = {
 135 /*0x0000,0x0000,0x0000,0x3ff0,*/
 136 0x6572,0xeeb3,0xb8a5,0x40ca,
 137 0xbc96,0x582a,0x27bc,0xc134,
 138 0xd8ef,0xc2ea,0xd98f,0x4177,
 139 0x5a31,0x3cbe,0xafe0,0xc189
 140 };
 141 /*
 142   7.85398125648498535156E-1,
 143   3.77489470793079817668E-8,
 144   2.69515142907905952645E-15,
 145 */
 146 static unsigned short P1[] = {0x0000,0x4000,0x21fb,0x3fe9};
 147 static unsigned short P2[] = {0x0000,0x0000,0x442d,0x3e64};
 148 static unsigned short P3[] = {0x5170,0x98cc,0x4698,0x3ce8};
 149 #define DP1 *(double *)P1
 150 #define DP2 *(double *)P2
 151 #define DP3 *(double *)P3
 152 static double lossth = 1.073741824e9;
 153 #endif
 154
 155 #ifdef MIEEE
 156 static unsigned short P[] = {
 157 0xc0c9,0x92d8,0xd24f,0x3f38,
 158 0x4131,0x99ec,0xa5fc,0x9ddd,
 159 0xc171,0x1fea,0xd329,0x9176
 160 };
 161 static unsigned short Q[] = {
 162 0x40ca,0xb8a5,0xeeb3,0x6572,
 163 0xc134,0x27bc,0x582a,0xbc96,
 164 0x4177,0xd98f,0xc2ea,0xd8ef,
 165 0xc189,0xafe0,0x3cbe,0x5a31
 166 };
 167 static unsigned short P1[] = {
 168 0x3fe9,0x21fb,0x4000,0x0000
 169 };
 170 static unsigned short P2[] = {
 171 0x3e64,0x442d,0x0000,0x0000
 172 };
 173 static unsigned short P3[] = {
 174 0x3ce8,0x4698,0x98cc,0x5170,
 175 };
 176 #define DP1 *(double *)P1
 177 #define DP2 *(double *)P2
 178 #define DP3 *(double *)P3
 179 static double lossth = 1.073741824e9;
 180 #endif
 181
 182 #ifdef ANSIPROT
 183 extern double polevl ( double, void *, int );
 184 extern double p1evl ( double, void *, int );
 185 extern double floor ( double );
 186 extern double ldexp ( double, int );
 187 extern int isnan ( double );
 188 extern int isfinite ( double );
 189 static double tancot(double, int);
 190 #else
 191 double polevl(), p1evl(), floor(), ldexp();
 192 static double tancot();
 193 int isnan(), isfinite();
 194 #endif
 195 extern double PIO4;
 196 extern double INFINITY;
 197 extern double NAN;
 198
 199 double tan(x)
 200 double x;
 201 {
 202 #ifdef MINUSZERO
 203 if( x == 0.0 )
 204         return(x);
 205 #endif
 206 #ifdef NANS
 207 if( isnan(x) )
 208         return(x);
 209 if( !isfinite(x) )
 210         {
 211         mtherr( "tan", DOMAIN );
 212         return(NAN);
 213         }
 214 #endif
 215 return( tancot(x,0) );
 216 }
 217
 218
 219 double cot(x)
 220 double x;
 221 {
 222
 223 if( x == 0.0 )
 224         {
 225         mtherr( "cot", SING );
 226         return( INFINITY );
 227         }
 228 return( tancot(x,1) );
 229 }
 230
 231
 232 static double tancot( xx, cotflg )
 233 double xx;
 234 int cotflg;
 235 {
 236 double x, y, z, zz;
 237 int j, sign;
 238
 239 /* make argument positive but save the sign */
 240 if( xx < 0 )
 241         {
 242         x = -xx;
 243         sign = -1;
 244         }
 245 else
 246         {
 247         x = xx;
 248         sign = 1;
 249         }
 250
 251 if( x > lossth )
 252         {
 253         if( cotflg )
 254                 mtherr( "cot", TLOSS );
 255         else
 256                 mtherr( "tan", TLOSS );
 257         return(0.0);
 258         }
 259
 260 /* compute x mod PIO4 */
 261 y = floor( x/PIO4 );
 262
 263 /* strip high bits of integer part */
 264 z = ldexp( y, -3 );
 265 z = floor(z);           /* integer part of y/8 */
 266 z = y - ldexp( z, 3 );  /* y - 16 * (y/16) */
 267
 268 /* integer and fractional part modulo one octant */
 269 j = z;
 270
 271 /* map zeros and singularities to origin */
 272 if( j & 1 )
 273         {
 274         j += 1;
 275         y += 1.0;
 276         }
 277
 278 z = ((x - y * DP1) - y * DP2) - y * DP3;
 279
 280 zz = z * z;
 281
 282 if( zz > 1.0e-14 )
 283         y = z  +  z * (zz * polevl( zz, P, 2 )/p1evl(zz, Q, 4));
 284 else
 285         y = z;
 286
 287 if( j & 2 )
 288         {
 289         if( cotflg )
 290                 y = -y;
 291         else
 292                 y = -1.0/y;
 293         }
 294 else
 295         {
 296         if( cotflg )
 297                 y = 1.0/y;
 298         }
 299
 300 if( sign < 0 )
 301         y = -y;
 302
 303 return( y );
 304 }