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