newlib/libm/mathfp/s_atangent.c

   1
   2 /* @(#)z_atangent.c 1.0 98/08/13 */
   3 /******************************************************************
   4  * The following routines are coded directly from the algorithms
   5  * and coefficients given in "Software Manual for the Elementary
   6  * Functions" by William J. Cody, Jr. and William Waite, Prentice
   7  * Hall, 1980.
   8  ******************************************************************/
   9
  10 /*
  11 FUNCTION
  12         <<atan>>, <<atanf>>, <<atan2>>, <<atan2f>>, <<atangent>>, <<atangentf>>---arc tangent
  13
  14 INDEX
  15    atan2
  16 INDEX
  17    atan2f
  18 INDEX
  19    atan
  20 INDEX
  21    atanf
  22
  23 SYNOPSIS
  24         #include <math.h>
  25         double atan(double <[x]>);
  26         float atan(float <[x]>);
  27         double atan2(double <[y]>,double <[x]>);
  28         float atan2f(float <[y]>,float <[x]>);
  29
  30 DESCRIPTION
  31
  32 <<atan2>> computes the inverse tangent (arc tangent) of y / x.
  33
  34 <<atan2f>> is identical to <<atan2>>, save that it operates on <<floats>>.
  35
  36 <<atan>> computes the inverse tangent (arc tangent) of the input value.
  37
  38 <<atanf>> is identical to <<atan>>, save that it operates on <<floats>>.
  39
  40 RETURNS
  41 @ifnottex
  42 <<atan>> returns a value in radians, in the range of -pi/2 to pi/2.
  43 <<atan2>> returns a value in radians, in the range of -pi/2 to pi/2.
  44 @end ifnottex
  45 @tex
  46 <<atan>> returns a value in radians, in the range of $-\pi/2$ to $\pi/2$.
  47 <<atan2>> returns a value in radians, in the range of $-\pi/2$ to $\pi/2$.
  48 @end tex
  49
  50 PORTABILITY
  51 <<atan>> is ANSI C.  <<atanf>> is an extension.
  52 <<atan2>> is ANSI C.  <<atan2f>> is an extension.
  53
  54 */
  55
  56 /******************************************************************
  57  * Arctangent
  58  *
  59  * Input:
  60  *   x - floating point value
  61  *
  62  * Output:
  63  *   arctangent of x
  64  *
  65  * Description:
  66  *   This routine calculates arctangents.
  67  *
  68  *****************************************************************/
  69 #include <float.h>
  70 #include "fdlibm.h"
  71 #include "zmath.h"
  72
  73 #ifndef _DOUBLE_IS_32BITS
  74
  75 static const double ROOT3 = 1.73205080756887729353;
  76 static const double a[] = { 0.0, 0.52359877559829887308, 1.57079632679489661923,
  77                      1.04719755119659774615 };
  78 static const double q[] = { 0.41066306682575781263e+2,
  79                      0.86157349597130242515e+2,
  80                      0.59578436142597344465e+2,
  81                      0.15024001160028576121e+2 };
  82 static const double p[] = { -0.13688768894191926929e+2,
  83                      -0.20505855195861651981e+2,
  84                      -0.84946240351320683534e+1,
  85                      -0.83758299368150059274 };
  86
  87 double
  88 atangent (double x,
  89         double v,
  90         double u,
  91         int arctan2)
  92 {
  93   double f, g, R, P, Q, A, res;
  94   int N;
  95   int branch = 0;
  96   int expv, expu;
  97
  98   /* Preparation for calculating arctan2. */
  99   if (arctan2)
 100     {
 101       if (u == 0.0)
 102         if (v == 0.0)
 103           {
 104             errno = ERANGE;
 105             return (z_notanum.d);
 106           }
 107         else
 108           {
 109             branch = 1;
 110             res = __PI_OVER_TWO;
 111           }
 112
 113       if (!branch)
 114         {
 115           int e;
 116           /* Get the exponent values of the inputs. */
 117           g = frexp (v, &expv);
 118           g = frexp (u, &expu);
 119
 120           /* See if a divide will overflow. */
 121           e = expv - expu;
 122           if (e > DBL_MAX_EXP)
 123             {
 124                branch = 1;
 125                res = __PI_OVER_TWO;
 126             }
 127
 128           /* Also check for underflow. */
 129           else if (e < DBL_MIN_EXP)
 130             {
 131                branch = 2;
 132                res = 0.0;
 133             }
 134          }
 135     }
 136
 137   if (!branch)
 138     {
 139       if (arctan2)
 140         f = fabs (v / u);
 141       else
 142         f = fabs (x);
 143
 144       if (f > 1.0)
 145         {
 146           f = 1.0 / f;
 147           N = 2;
 148         }
 149       else
 150         N = 0;
 151
 152       if (f > (2.0 - ROOT3))
 153         {
 154           A = ROOT3 - 1.0;
 155           f = (((A * f - 0.5) - 0.5) + f) / (ROOT3 + f);
 156           N++;
 157         }
 158
 159       /* Check for values that are too small. */
 160       if (-z_rooteps < f && f < z_rooteps)
 161         res = f;
 162
 163       /* Calculate the Taylor series. */
 164       else
 165         {
 166           g = f * f;
 167           P = (((p[3] * g + p[2]) * g + p[1]) * g + p[0]) * g;
 168           Q = (((g + q[3]) * g + q[2]) * g + q[1]) * g + q[0];
 169           R = P / Q;
 170
 171           res = f + f * R;
 172         }
 173
 174       if (N > 1)
 175         res = -res;
 176
 177       res += a[N];
 178     }
 179
 180   if (arctan2)
 181     {
 182       if (u < 0.0)
 183         res = __PI - res;
 184       if (v < 0.0)
 185         res = -res;
 186     }
 187   else if (x < 0.0)
 188     {
 189       res = -res;
 190     }
 191
 192   return (res);
 193 }
 194
 195 #endif /* _DOUBLE_IS_32BITS */