text: 
  /* @(#)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 */
  