libm/ldouble/atanl.c

   1 /*                                                      atanl.c
   2  *
   3  *      Inverse circular tangent, long double precision
   4  *      (arctangent)
   5  *
   6  *
   7  *
   8  * SYNOPSIS:
   9  *
  10  * long double x, y, atanl();
  11  *
  12  * y = atanl( x );
  13  *
  14  *
  15  *
  16  * DESCRIPTION:
  17  *
  18  * Returns radian angle between -pi/2 and +pi/2 whose tangent
  19  * is x.
  20  *
  21  * Range reduction is from four intervals into the interval
  22  * from zero to  tan( pi/8 ).  The approximant uses a rational
  23  * function of degree 3/4 of the form x + x**3 P(x)/Q(x).
  24  *
  25  *
  26  *
  27  * ACCURACY:
  28  *
  29  *                      Relative error:
  30  * arithmetic   domain     # trials      peak         rms
  31  *    IEEE      -10, 10    150000       1.3e-19     3.0e-20
  32  *
  33  */
  34 \f/*                                                     atan2l()
  35  *
  36  *      Quadrant correct inverse circular tangent,
  37  *      long double precision
  38  *
  39  *
  40  *
  41  * SYNOPSIS:
  42  *
  43  * long double x, y, z, atan2l();
  44  *
  45  * z = atan2l( y, x );
  46  *
  47  *
  48  *
  49  * DESCRIPTION:
  50  *
  51  * Returns radian angle whose tangent is y/x.
  52  * Define compile time symbol ANSIC = 1 for ANSI standard,
  53  * range -PI < z <= +PI, args (y,x); else ANSIC = 0 for range
  54  * 0 to 2PI, args (x,y).
  55  *
  56  *
  57  *
  58  * ACCURACY:
  59  *
  60  *                      Relative error:
  61  * arithmetic   domain     # trials      peak         rms
  62  *    IEEE      -10, 10     60000       1.7e-19     3.2e-20
  63  * See atan.c.
  64  *
  65  */
  66 \f
  67 /*                                                      atan.c */
  68
  69
  70 /*
  71 Cephes Math Library Release 2.7:  May, 1998
  72 Copyright 1984, 1990, 1998 by Stephen L. Moshier
  73 */
  74
  75
  76 #include <math.h>
  77
  78 #ifdef UNK
  79 static long double P[] = {
  80 -8.6863818178092187535440E-1L,
  81 -1.4683508633175792446076E1L,
  82 -6.3976888655834347413154E1L,
  83 -9.9988763777265819915721E1L,
  84 -5.0894116899623603312185E1L,
  85 };
  86 static long double Q[] = {
  87 /* 1.00000000000000000000E0L,*/
  88  2.2981886733594175366172E1L,
  89  1.4399096122250781605352E2L,
  90  3.6144079386152023162701E2L,
  91  3.9157570175111990631099E2L,
  92  1.5268235069887081006606E2L,
  93 };
  94
  95 /* tan( 3*pi/8 ) */
  96 static long double T3P8 =  2.41421356237309504880169L;
  97
  98 /* tan( pi/8 ) */
  99 static long double TP8 =  4.1421356237309504880169e-1L;
 100 #endif
 101
 102
 103 #ifdef IBMPC
 104 static unsigned short P[] = {
 105 0x8ece,0xce53,0x1266,0xde5f,0xbffe, XPD
 106 0x07e6,0xa061,0xa6bf,0xeaef,0xc002, XPD
 107 0x53ee,0xf291,0x557f,0xffe8,0xc004, XPD
 108 0xf9d6,0xeda6,0x3f3e,0xc7fa,0xc005, XPD
 109 0xb6c3,0x6abc,0x9361,0xcb93,0xc004, XPD
 110 };
 111 static unsigned short Q[] = {
 112 /*0x0000,0x0000,0x0000,0x8000,0x3fff,*/
 113 0x54d4,0x894e,0xe76e,0xb7da,0x4003, XPD
 114 0x76b9,0x7a46,0xafa2,0x8ffd,0x4006, XPD
 115 0xe3a9,0xe9c0,0x6bee,0xb4b8,0x4007, XPD
 116 0xabc1,0x50a7,0xb098,0xc3c9,0x4007, XPD
 117 0x891c,0x100d,0xae89,0x98ae,0x4006, XPD
 118 };
 119
 120 /* tan( 3*pi/8 ) = 2.41421356237309504880 */
 121 static unsigned short T3P8A[] = {0x3242,0xfcef,0x7999,0x9a82,0x4000, XPD};
 122 #define T3P8 *(long double *)T3P8A
 123
 124 /* tan( pi/8 ) = 0.41421356237309504880 */
 125 static unsigned short TP8A[] = {0x9211,0xe779,0xcccf,0xd413,0x3ffd, XPD};
 126 #define TP8 *(long double *)TP8A
 127 #endif
 128
 129 #ifdef MIEEE
 130 static unsigned long P[] = {
 131 0xbffe0000,0xde5f1266,0xce538ece,
 132 0xc0020000,0xeaefa6bf,0xa06107e6,
 133 0xc0040000,0xffe8557f,0xf29153ee,
 134 0xc0050000,0xc7fa3f3e,0xeda6f9d6,
 135 0xc0040000,0xcb939361,0x6abcb6c3,
 136 };
 137 static unsigned long Q[] = {
 138 /*0x3fff0000,0x80000000,0x00000000,*/
 139 0x40030000,0xb7dae76e,0x894e54d4,
 140 0x40060000,0x8ffdafa2,0x7a4676b9,
 141 0x40070000,0xb4b86bee,0xe9c0e3a9,
 142 0x40070000,0xc3c9b098,0x50a7abc1,
 143 0x40060000,0x98aeae89,0x100d891c,
 144 };
 145
 146 /* tan( 3*pi/8 ) = 2.41421356237309504880 */
 147 static long T3P8A[] = {0x40000000,0x9a827999,0xfcef3242};
 148 #define T3P8 *(long double *)T3P8A
 149
 150 /* tan( pi/8 ) = 0.41421356237309504880 */
 151 static long TP8A[] = {0x3ffd0000,0xd413cccf,0xe7799211};
 152 #define TP8 *(long double *)TP8A
 153 #endif
 154
 155 #ifdef ANSIPROT
 156 extern long double polevll ( long double, void *, int );
 157 extern long double p1evll ( long double, void *, int );
 158 extern long double fabsl ( long double );
 159 extern int signbitl ( long double );
 160 extern int isnanl ( long double );
 161 long double atanl ( long double );
 162 #else
 163 long double polevll(), p1evll(), fabsl(), signbitl(), isnanl();
 164 long double atanl();
 165 #endif
 166 #ifdef INFINITIES
 167 extern long double INFINITYL;
 168 #endif
 169 #ifdef NANS
 170 extern long double NANL;
 171 #endif
 172 #ifdef MINUSZERO
 173 extern long double NEGZEROL;
 174 #endif
 175
 176 long double atanl(x)
 177 long double x;
 178 {
 179 extern long double PIO2L, PIO4L;
 180 long double y, z;
 181 short sign;
 182
 183 #ifdef MINUSZERO
 184 if( x == 0.0L )
 185         return(x);
 186 #endif
 187 #ifdef INFINITIES
 188 if( x == INFINITYL )
 189         return( PIO2L );
 190 if( x == -INFINITYL )
 191         return( -PIO2L );
 192 #endif
 193 /* make argument positive and save the sign */
 194 sign = 1;
 195 if( x < 0.0L )
 196         {
 197         sign = -1;
 198         x = -x;
 199         }
 200
 201 /* range reduction */
 202 if( x > T3P8 )
 203         {
 204         y = PIO2L;
 205         x = -( 1.0L/x );
 206         }
 207
 208 else if( x > TP8 )
 209         {
 210         y = PIO4L;
 211         x = (x-1.0L)/(x+1.0L);
 212         }
 213 else
 214         y = 0.0L;
 215
 216 /* rational form in x**2 */
 217 z = x * x;
 218 y = y + ( polevll( z, P, 4 ) / p1evll( z, Q, 5 ) ) * z * x + x;
 219
 220 if( sign < 0 )
 221         y = -y;
 222
 223 return(y);
 224 }
 225 \f
 226 /*                                                      atan2   */
 227
 228
 229 extern long double PIL, PIO2L, MAXNUML;
 230
 231 #if ANSIC
 232 long double atan2l( y, x )
 233 #else
 234 long double atan2l( x, y )
 235 #endif
 236 long double x, y;
 237 {
 238 long double z, w;
 239 short code;
 240
 241 code = 0;
 242
 243 if( x < 0.0L )
 244         code = 2;
 245 if( y < 0.0L )
 246         code |= 1;
 247
 248 #ifdef NANS
 249 if( isnanl(x) )
 250         return(x);
 251 if( isnanl(y) )
 252         return(y);
 253 #endif
 254 #ifdef MINUSZERO
 255 if( y == 0.0L )
 256         {
 257         if( signbitl(y) )
 258                 {
 259                 if( x > 0.0L )
 260                         z = y;
 261                 else if( x < 0.0L )
 262                         z = -PIL;
 263                 else
 264                         {
 265                         if( signbitl(x) )
 266                                 z = -PIL;
 267                         else
 268                                 z = y;
 269                         }
 270                 }
 271         else /* y is +0 */
 272                 {
 273                 if( x == 0.0L )
 274                         {
 275                         if( signbitl(x) )
 276                                 z = PIL;
 277                         else
 278                                 z = 0.0L;
 279                         }
 280                 else if( x > 0.0L )
 281                         z = 0.0L;
 282                 else
 283                         z = PIL;
 284                 }
 285         return z;
 286         }
 287 if( x == 0.0L )
 288         {
 289         if( y > 0.0L )
 290                 z = PIO2L;
 291         else
 292                 z = -PIO2L;
 293         return z;
 294         }
 295 #endif /* MINUSZERO */
 296 #ifdef INFINITIES
 297 if( x == INFINITYL )
 298         {
 299         if( y == INFINITYL )
 300                 z = 0.25L * PIL;
 301         else if( y == -INFINITYL )
 302                 z = -0.25L * PIL;
 303         else if( y < 0.0L )
 304                 z = NEGZEROL;
 305         else
 306                 z = 0.0L;
 307         return z;
 308         }
 309 if( x == -INFINITYL )
 310         {
 311         if( y == INFINITYL )
 312                 z = 0.75L * PIL;
 313         else if( y == -INFINITYL )
 314                 z = -0.75L * PIL;
 315         else if( y >= 0.0L )
 316                 z = PIL;
 317         else
 318                 z = -PIL;
 319         return z;
 320         }
 321 if( y == INFINITYL )
 322         return( PIO2L );
 323 if( y == -INFINITYL )
 324         return( -PIO2L );
 325 #endif /* INFINITIES */
 326
 327 #ifdef INFINITIES
 328 if( x == 0.0L )
 329 #else
 330 if( fabsl(x) <= (fabsl(y) / MAXNUML) )
 331 #endif
 332         {
 333         if( code & 1 )
 334                 {
 335 #if ANSIC
 336                 return( -PIO2L );
 337 #else
 338                 return( 3.0L*PIO2L );
 339 #endif
 340                 }
 341         if( y == 0.0L )
 342                 return( 0.0L );
 343         return( PIO2L );
 344         }
 345
 346 if( y == 0.0L )
 347         {
 348         if( code & 2 )
 349                 return( PIL );
 350         return( 0.0L );
 351         }
 352
 353
 354 switch( code )
 355         {
 356         default:
 357 #if ANSIC
 358         case 0:
 359         case 1: w = 0.0L; break;
 360         case 2: w = PIL; break;
 361         case 3: w = -PIL; break;
 362 #else
 363         case 0: w = 0.0L; break;
 364         case 1: w = 2.0L * PIL; break;
 365         case 2:
 366         case 3: w = PIL; break;
 367 #endif
 368         }
 369
 370 z = w + atanl( y/x );
 371 #ifdef MINUSZERO
 372 if( z == 0.0L && y < 0.0L )
 373         z = NEGZEROL;
 374 #endif
 375 return( z );
 376 }