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)
201         res = __PI - res;
202       if (v < 0.0)
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 */
214