Subversion Repositories or1k

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

ANSI_SYNOPSIS

        #include <math.h>

        double atan(double <[x]>);

        float atan(float <[x]>);

        double atan2(double <[y]>,double <[x]>);

        float atan2f(float <[y]>,float <[x]>);

TRAD_SYNOPSIS

        #include <math.h>

        double atan2(<[y]>,<[x]>);

        double <[y]>;

        double <[x]>;

        float atan2f(<[y]>,<[x]>);

        float <[y]>;

        float <[x]>;

        #include <math.h>

        double atan(<[x]>);

        double <[x]>;

        float atanf(<[x]>);

        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

@ifinfo

<<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 ifinfo

@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)

        {

          /* Get the exponent values of the inputs. */

          g = frexp (v, &expv);

          g = frexp (u, &expu);

          /* See if a divide will overflow. */

          if ((expv - expu) > DBL_MAX_EXP + 1023)

            {

               branch = 1;

               res = __PI_OVER_TWO;

            }

          /* Also check for underflow. */

          else if ((expv - expu) < -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 && branch)

    {

      if (u < 0.0 || branch == 2)

        res = __PI - res;

      if (v < 0.0 || branch == 1)

        res = -res;

    }

  else if (x < 0.0)

    {

      res = -res;

    }

  return (res);

}

#endif /* _DOUBLE_IS_32BITS */