CocoRoboDesktop/pdf1.js/node_modules/cephes/cephes/atan.c

/*							atan.c
 *
 *	Inverse circular tangent
 *      (arctangent)
 *
 *
 *
 * SYNOPSIS:
 *
 * double x, y, atan();
 *
 * y = atan( x );
 *
 *
 *
 * DESCRIPTION:
 *
 * Returns radian angle between -pi/2 and +pi/2 whose tangent
 * is x.
 *
 * Range reduction is from three intervals into the interval
 * from zero to 0.66.  The approximant uses a rational
 * function of degree 4/5 of the form x + x**3 P(x)/Q(x).
 *
 *
 *
 * ACCURACY:
 *
 *                      Relative error:
 * arithmetic   domain     # trials      peak         rms
 *    DEC       -10, 10     50000       2.4e-17     8.3e-18
 *    IEEE      -10, 10      10^6       1.8e-16     5.0e-17
 *
 */
/*							atan2()
 *
 *	Quadrant correct inverse circular tangent
 *
 *
 *
 * SYNOPSIS:
 *
 * double x, y, z, atan2();
 *
 * z = atan2( y, x );
 *
 *
 *
 * DESCRIPTION:
 *
 * Returns radian angle whose tangent is y/x.
 * Define compile time symbol ANSIC = 1 for ANSI standard,
 * range -PI < z <= +PI, args (y,x); else ANSIC = 0 for range
 * 0 to 2PI, args (x,y).
 *
 *
 *
 * ACCURACY:
 *
 *                      Relative error:
 * arithmetic   domain     # trials      peak         rms
 *    IEEE      -10, 10      10^6       2.5e-16     6.9e-17
 * See atan.c.
 *
 */

/*							atan.c */

/*
Cephes Math Library Release 2.8:  June, 2000
Copyright 1984, 1995, 2000 by Stephen L. Moshier
*/

#include "mconf.h"

/* arctan(x)  = x + x^3 P(x^2)/Q(x^2)
   0 <= x <= 0.66
   Peak relative error = 2.6e-18  */
#ifdef UNK
static double P[5] = {
    -8.750608600031904122785E-1, -1.615753718733365076637E1,
    -7.500855792314704667340E1,  -1.228866684490136173410E2,
    -6.485021904942025371773E1,
};
static double Q[5] = {
    /* 1.000000000000000000000E0, */
    2.485846490142306297962E1, 1.650270098316988542046E2,
    4.328810604912902668951E2, 4.853903996359136964868E2,
    1.945506571482613964425E2,
};

/* tan( 3*pi/8 ) */
static double T3P8 = 2.41421356237309504880;
#endif

#ifdef DEC
static short P[20] = {
    0140140, 0001775, 0007671, 0026242, 0141201, 0041242, 0155534,
    0001715, 0141626, 0002141, 0132100, 0011625, 0141765, 0142771,
    0064055, 0150453, 0141601, 0131517, 0164507, 0062164,
};
static short Q[20] = {
    /* 0040200,0000000,0000000,0000000, */
    0041306, 0157042, 0154243, 0000742, 0042045, 0003352, 0016707,
    0150452, 0042330, 0070306, 0113425, 0170730, 0042362, 0130770,
    0116602, 0047520, 0042102, 0106367, 0156753, 0013541,
};

/* tan( 3*pi/8 ) = 2.41421356237309504880 */
static unsigned short T3P8A[] = {
    040432,
    0101171,
    0114774,
    0167462,
};
#define T3P8 *(double *)T3P8A
#endif

#ifdef IBMPC
static short P[20] = {
    0x2594, 0xa1f7, 0x007f, 0xbfec, 0x807a, 0x5b6b, 0x2854,
    0xc030, 0x0273, 0x3688, 0xc08c, 0xc052, 0xba25, 0x2d05,
    0xb8bf, 0xc05e, 0xec8e, 0xfd28, 0x3669, 0xc050,
};
static short Q[20] = {
    /* 0x0000,0x0000,0x0000,0x3ff0, */
    0x603c, 0x5b14, 0xdbc4, 0x4038, 0xfa25, 0x43b8, 0xa0dd,
    0x4064, 0xbe3b, 0xd2e2, 0x0e18, 0x407b, 0x49ea, 0x13b0,
    0x563f, 0x407e, 0x62ec, 0xfbbd, 0x519e, 0x4068,
};

/* tan( 3*pi/8 ) = 2.41421356237309504880 */
static unsigned short T3P8A[] = {0x9de6, 0x333f, 0x504f, 0x4003};
#define T3P8 *(double *)T3P8A
#endif

#ifdef MIEEE
static short P[20] = {
    0xbfec, 0x007f, 0xa1f7, 0x2594, 0xc030, 0x2854, 0x5b6b,
    0x807a, 0xc052, 0xc08c, 0x3688, 0x0273, 0xc05e, 0xb8bf,
    0x2d05, 0xba25, 0xc050, 0x3669, 0xfd28, 0xec8e,
};
static short Q[20] = {
    /* 0x3ff0,0x0000,0x0000,0x0000, */
    0x4038, 0xdbc4, 0x5b14, 0x603c, 0x4064, 0xa0dd, 0x43b8,
    0xfa25, 0x407b, 0x0e18, 0xd2e2, 0xbe3b, 0x407e, 0x563f,
    0x13b0, 0x49ea, 0x4068, 0x519e, 0xfbbd, 0x62ec,
};

/* tan( 3*pi/8 ) = 2.41421356237309504880 */
static unsigned short T3P8A[] = {0x4003, 0x504f, 0x333f, 0x9de6};
#define T3P8 *(double *)T3P8A
#endif

#ifdef ANSIPROT
extern double polevl(double, void *, int);
extern double p1evl(double, void *, int);
extern double atan(double);
extern double fabs(double);
extern int signbit(double);
extern int isnan(double);
#else
double polevl(), p1evl(), atan(), fabs();
int signbit(), isnan();
#endif
extern double PI, PIO2, PIO4, INFINITY, NEGZERO, MAXNUM;

/* pi/2 = PIO2 + MOREBITS.  */
#ifdef DEC
#define MOREBITS 5.721188726109831840122E-18
#else
#define MOREBITS 6.123233995736765886130E-17
#endif

double atan(x) double x;
{
  double y, z;
  short sign, flag;

#ifdef MINUSZERO
  if (x == 0.0)
    return (x);
#endif
#ifdef INFINITIES
  if (x == INFINITY)
    return (PIO2);
  if (x == -INFINITY)
    return (-PIO2);
#endif
  /* make argument positive and save the sign */
  sign = 1;
  if (x < 0.0) {
    sign = -1;
    x = -x;
  }
  /* range reduction */
  flag = 0;
  if (x > T3P8) {
    y = PIO2;
    flag = 1;
    x = -(1.0 / x);
  } else if (x <= 0.66) {
    y = 0.0;
  } else {
    y = PIO4;
    flag = 2;
    x = (x - 1.0) / (x + 1.0);
  }
  z = x * x;
  z = z * polevl(z, P, 4) / p1evl(z, Q, 5);
  z = x * z + x;
  if (flag == 2)
    z += 0.5 * MOREBITS;
  else if (flag == 1)
    z += MOREBITS;
  y = y + z;
  if (sign < 0)
    y = -y;
  return (y);
}

/*							atan2	*/

#ifdef ANSIC
double atan2(y, x)
#else
double atan2(x, y)
#endif
    double x,
    y;
{
  double z, w;
  short code;

  code = 0;

#ifdef NANS
  if (isnan(x))
    return (x);
  if (isnan(y))
    return (y);
#endif
#ifdef MINUSZERO
  if (y == 0.0) {
    if (signbit(y)) {
      if (x > 0.0)
        z = y;
      else if (x < 0.0)
        z = -PI;
      else {
        if (signbit(x))
          z = -PI;
        else
          z = y;
      }
    } else /* y is +0 */
    {
      if (x == 0.0) {
        if (signbit(x))
          z = PI;
        else
          z = 0.0;
      } else if (x > 0.0)
        z = 0.0;
      else
        z = PI;
    }
    return z;
  }
  if (x == 0.0) {
    if (y > 0.0)
      z = PIO2;
    else
      z = -PIO2;
    return z;
  }
#endif /* MINUSZERO */
#ifdef INFINITIES
  if (x == INFINITY) {
    if (y == INFINITY)
      z = 0.25 * PI;
    else if (y == -INFINITY)
      z = -0.25 * PI;
    else if (y < 0.0)
      z = NEGZERO;
    else
      z = 0.0;
    return z;
  }
  if (x == -INFINITY) {
    if (y == INFINITY)
      z = 0.75 * PI;
    else if (y <= -INFINITY)
      z = -0.75 * PI;
    else if (y >= 0.0)
      z = PI;
    else
      z = -PI;
    return z;
  }
  if (y == INFINITY)
    return (PIO2);
  if (y == -INFINITY)
    return (-PIO2);
#endif

  if (x < 0.0)
    code = 2;
  if (y < 0.0)
    code |= 1;

#ifdef INFINITIES
  if (x == 0.0)
#else
  if (fabs(x) <= (fabs(y) / MAXNUM))
#endif
  {
    if (code & 1) {
#if ANSIC
      return (-PIO2);
#else
      return (3.0 * PIO2);
#endif
    }
    if (y == 0.0)
      return (0.0);
    return (PIO2);
  }

  if (y == 0.0) {
    if (code & 2)
      return (PI);
    return (0.0);
  }

  switch (code) {
#if ANSIC
  default:
  case 0:
  case 1:
    w = 0.0;
    break;
  case 2:
    w = PI;
    break;
  case 3:
    w = -PI;
    break;
#else
  default:
  case 0:
    w = 0.0;
    break;
  case 1:
    w = 2.0 * PI;
    break;
  case 2:
  case 3:
    w = PI;
    break;
#endif
  }

  z = w + atan(y / x);
#ifdef MINUSZERO
  if (z == 0.0 && y < 0)
    z = NEGZERO;
#endif
  return (z);
}