diff options
Diffstat (limited to 'newlib/libm/mathfp/s_atangent.c')
| -rw-r--r-- | newlib/libm/mathfp/s_atangent.c | 213 |
1 files changed, 213 insertions, 0 deletions
diff --git a/newlib/libm/mathfp/s_atangent.c b/newlib/libm/mathfp/s_atangent.c new file mode 100644 index 00000000000..c6f3c9bd6bc --- /dev/null +++ b/newlib/libm/mathfp/s_atangent.c @@ -0,0 +1,213 @@ + +/* @(#)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) + { + int e; + /* Get the exponent values of the inputs. */ + g = frexp (v, &expv); + g = frexp (u, &expu); + + /* See if a divide will overflow. */ + e = expv - expu; + if (e > DBL_MAX_EXP) + { + branch = 1; + res = __PI_OVER_TWO; + } + + /* Also check for underflow. */ + else if (e < 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) + { + 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 */ |