/* @(#)z_atangent.c 1.0 98/08/13 */
/******************************************************************
 * The following routines are coded directly from the algorithms
 * and coefficients given in "Software Manual for the Elementary
 * Functions" by William J. Cody, Jr. and William Waite, Prentice
 * Hall, 1980.
 ******************************************************************/

/*
FUNCTION
        <<atan>>, <<atanf>>, <<atan2>>, <<atan2f>>, <<atangent>>, <<atangentf>>---arc tangent

INDEX
   atan2
INDEX
   atan2f
INDEX
   atan
INDEX
   atanf

SYNOPSIS
        #include <math.h>
        double atan(double <[x]>);
        float atan(float <[x]>);
        double atan2(double <[y]>,double <[x]>);
        float atan2f(float <[y]>,float <[x]>);

DESCRIPTION

<<atan2>> computes the inverse tangent (arc tangent) of y / x.

<<atan2f>> is identical to <<atan2>>, save that it operates on <<floats>>.

<<atan>> computes the inverse tangent (arc tangent) of the input value.

<<atanf>> is identical to <<atan>>, save that it operates on <<floats>>.

RETURNS
@ifnottex
<<atan>> returns a value in radians, in the range of -pi/2 to pi/2.
<<atan2>> returns a value in radians, in the range of -pi/2 to pi/2.
@end ifnottex
@tex
<<atan>> returns a value in radians, in the range of $-\pi/2$ to $\pi/2$.
<<atan2>> returns a value in radians, in the range of $-\pi/2$ to $\pi/2$.
@end tex

PORTABILITY
<<atan>> is ANSI C.  <<atanf>> is an extension.
<<atan2>> is ANSI C.  <<atan2f>> is an extension.

*/

/******************************************************************
 * Arctangent
 *
 * Input:
 *   x - floating point value
 *
 * Output:
 *   arctangent of x
 *
 * Description:
 *   This routine calculates arctangents.
 *
 *****************************************************************/
#include <float.h>
#include "fdlibm.h"
#include "zmath.h"

#ifndef _DOUBLE_IS_32BITS

static const double ROOT3 = 1.73205080756887729353;
static const double a[] = { 0.0, 0.52359877559829887308, 1.57079632679489661923,
                     1.04719755119659774615 };
static const double q[] = { 0.41066306682575781263e+2,
                     0.86157349597130242515e+2,
                     0.59578436142597344465e+2,
                     0.15024001160028576121e+2 };
static const double p[] = { -0.13688768894191926929e+2,
                     -0.20505855195861651981e+2,
                     -0.84946240351320683534e+1,
                     -0.83758299368150059274 };

double
_DEFUN (atangent, (double, double, double, int),
        double x _AND
        double v _AND
        double u _AND
        int arctan2)
{
  double f, g, R, P, Q, A, res;
  int N;
  int branch = 0;
  int expv, expu;

  /* Preparation for calculating arctan2. */
  if (arctan2)
    {
      if (u == 0.0)
        if (v == 0.0)
          {
            errno = ERANGE;
            return (z_notanum.d);
          }
        else
          {
            branch = 1;
            res = __PI_OVER_TWO;
          }

      if (!branch)
        {
          int e;
          /* Get the exponent values of the inputs. */
          g = frexp (v, &expv);
          g = frexp (u, &expu);

          /* See if a divide will overflow. */
          e = expv - expu;
          if (e > DBL_MAX_EXP)
            {
               branch = 1;
               res = __PI_OVER_TWO;
            }

          /* Also check for underflow. */
          else if (e < DBL_MIN_EXP)
            {
               branch = 2;
               res = 0.0;
            }
         }
    }

  if (!branch)
    {
      if (arctan2)
        f = fabs (v / u);
      else
        f = fabs (x);

      if (f > 1.0)
        {
          f = 1.0 / f;
          N = 2;
        }
      else
        N = 0;

      if (f > (2.0 - ROOT3))
        {
          A = ROOT3 - 1.0;
          f = (((A * f - 0.5) - 0.5) + f) / (ROOT3 + f);
          N++;
        }

      /* Check for values that are too small. */
      if (-z_rooteps < f && f < z_rooteps)
        res = f;

      /* Calculate the Taylor series. */
      else
        {
          g = f * f;
          P = (((p[3] * g + p[2]) * g + p[1]) * g + p[0]) * g;
          Q = (((g + q[3]) * g + q[2]) * g + q[1]) * g + q[0];
          R = P / Q;

          res = f + f * R;
        }

      if (N > 1)
        res = -res;

      res += a[N];
    }

  if (arctan2)
    {
      if (u < 0.0)
        res = __PI - res;
      if (v < 0.0)
        res = -res;
    }
  else if (x < 0.0)
    {
      res = -res;
    }

  return (res);
}

#endif /* _DOUBLE_IS_32BITS */