newlib-cygwin/newlib/libm/mathfp/sf_atangent.c

141 lines
2.9 KiB
C

/* @(#)z_atangentf.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.
******************************************************************/
/******************************************************************
* Arctangent
*
* Input:
* x - floating point value
*
* Output:
* arctangent of x
*
* Description:
* This routine calculates arctangents.
*
*****************************************************************/
#include <float.h>
#include "fdlibm.h"
#include "zmath.h"
static const float ROOT3 = 1.732050807;
static const float a[] = { 0.0, 0.523598775, 1.570796326,
1.047197551 };
static const float q[] = { 0.1412500740e+1 };
static const float p[] = { -0.4708325141, -0.5090958253e-1 };
float
_DEFUN (atangentf, (float, float, float, int),
float x _AND
float v _AND
float u _AND
int arctan2)
{
float 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_f.f);
}
else
{
branch = 1;
res = __PI_OVER_TWO;
}
if (!branch)
{
int e;
/* Get the exponent values of the inputs. */
g = frexpf (v, &expv);
g = frexpf (u, &expu);
/* See if a divide will overflow. */
e = expv - expu;
if (e > FLT_MAX_EXP)
{
branch = 1;
res = __PI_OVER_TWO;
}
/* Also check for underflow. */
else if (e < FLT_MIN_EXP)
{
branch = 2;
res = 0.0;
}
}
}
if (!branch)
{
if (arctan2)
f = fabsf (v / u);
else
f = fabsf (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 && f < z_rooteps_f)
res = f;
/* Calculate the Taylor series. */
else
{
g = f * f;
P = (p[1] * g + p[0]) * g;
Q = 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);
}