toolchain/gcc/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 ANSI_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 TRAD_SYNOPSIS
  31         #include <math.h>
  32         double atan2(<[y]>,<[x]>);
  33         double <[y]>;
  34         double <[x]>;
  35
  36         float atan2f(<[y]>,<[x]>);
  37         float <[y]>;
  38         float <[x]>;
  39
  40         #include <math.h>
  41         double atan(<[x]>);
  42         double <[x]>;
  43
  44         float atanf(<[x]>);
  45         float <[x]>;
  46
  47 DESCRIPTION
  48
  49 <<atan2>> computes the inverse tangent (arc tangent) of y / x.
  50
  51 <<atan2f>> is identical to <<atan2>>, save that it operates on <<floats>>.
  52
  53 <<atan>> computes the inverse tangent (arc tangent) of the input value.
  54
  55 <<atanf>> is identical to <<atan>>, save that it operates on <<floats>>.
  56
  57 RETURNS
  58 @ifnottex
  59 <<atan>> returns a value in radians, in the range of -pi/2 to pi/2.
  60 <<atan2>> returns a value in radians, in the range of -pi/2 to pi/2.
  61 @end ifnottex
  62 @tex
  63 <<atan>> returns a value in radians, in the range of $-\pi/2$ to $\pi/2$.
  64 <<atan2>> returns a value in radians, in the range of $-\pi/2$ to $\pi/2$.
  65 @end tex
  66
  67 PORTABILITY
  68 <<atan>> is ANSI C.  <<atanf>> is an extension.
  69 <<atan2>> is ANSI C.  <<atan2f>> is an extension.
  70
  71 */
  72
  73 /******************************************************************
  74  * Arctangent
  75  *
  76  * Input:
  77  *   x - floating point value
  78  *
  79  * Output:
  80  *   arctangent of x
  81  *
  82  * Description:
  83  *   This routine calculates arctangents.
  84  *
  85  *****************************************************************/
  86 #include <float.h>
  87 #include "fdlibm.h"
  88 #include "zmath.h"
  89
  90 #ifndef _DOUBLE_IS_32BITS
  91
  92 static const double ROOT3 = 1.73205080756887729353;
  93 static const double a[] = { 0.0, 0.52359877559829887308, 1.57079632679489661923,
  94                      1.04719755119659774615 };
  95 static const double q[] = { 0.41066306682575781263e+2,
  96                      0.86157349597130242515e+2,
  97                      0.59578436142597344465e+2,
  98                      0.15024001160028576121e+2 };
  99 static const double p[] = { -0.13688768894191926929e+2,
 100                      -0.20505855195861651981e+2,
 101                      -0.84946240351320683534e+1,
 102                      -0.83758299368150059274 };
 103
 104 double
 105 _DEFUN (atangent, (double, double, double, int),
 106         double x _AND
 107         double v _AND
 108         double u _AND
 109         int arctan2)
 110 {
 111   double f, g, R, P, Q, A, res;
 112   int N;
 113   int branch = 0;
 114   int expv, expu;
 115
 116   /* Preparation for calculating arctan2. */
 117   if (arctan2)
 118     {
 119       if (u == 0.0)
 120         if (v == 0.0)
 121           {
 122             errno = ERANGE;
 123             return (z_notanum.d);
 124           }
 125         else
 126           {
 127             branch = 1;
 128             res = __PI_OVER_TWO;
 129           }
 130
 131       if (!branch)
 132         {
 133           int e;
 134           /* Get the exponent values of the inputs. */
 135           g = frexp (v, &expv);
 136           g = frexp (u, &expu);
 137
 138           /* See if a divide will overflow. */
 139           e = expv - expu;
 140           if (e > DBL_MAX_EXP)
 141             {
 142                branch = 1;
 143                res = __PI_OVER_TWO;
 144             }
 145
 146           /* Also check for underflow. */
 147           else if (e < DBL_MIN_EXP)
 148             {
 149                branch = 2;
 150                res = 0.0;
 151             }
 152          }
 153     }
 154
 155   if (!branch)
 156     {
 157       if (arctan2)
 158         f = fabs (v / u);
 159       else
 160         f = fabs (x);
 161
 162       if (f > 1.0)
 163         {
 164           f = 1.0 / f;
 165           N = 2;
 166         }
 167       else
 168         N = 0;
 169
 170       if (f > (2.0 - ROOT3))
 171         {
 172           A = ROOT3 - 1.0;
 173           f = (((A * f - 0.5) - 0.5) + f) / (ROOT3 + f);
 174           N++;
 175         }
 176
 177       /* Check for values that are too small. */
 178       if (-z_rooteps < f && f < z_rooteps)
 179         res = f;
 180
 181       /* Calculate the Taylor series. */
 182       else
 183         {
 184           g = f * f;
 185           P = (((p[3] * g + p[2]) * g + p[1]) * g + p[0]) * g;
 186           Q = (((g + q[3]) * g + q[2]) * g + q[1]) * g + q[0];
 187           R = P / Q;
 188
 189           res = f + f * R;
 190         }
 191
 192       if (N > 1)
 193         res = -res;
 194
 195       res += a[N];
 196     }
 197
 198   if (arctan2)
 199     {
 200       if (u < 0.0 || branch == 2)
 201         res = __PI - res;
 202       if (v < 0.0 || branch == 1)
 203         res = -res;
 204     }
 205   else if (x < 0.0)
 206     {
 207       res = -res;
 208     }
 209
 210   return (res);
 211 }
 212
 213 #endif /* _DOUBLE_IS_32BITS */