diff options
Diffstat (limited to 'libm/double/atan.c')
| -rw-r--r-- | libm/double/atan.c | 393 |
1 files changed, 393 insertions, 0 deletions
diff --git a/libm/double/atan.c b/libm/double/atan.c new file mode 100644 index 000000000..f2d50768d --- /dev/null +++ b/libm/double/atan.c @@ -0,0 +1,393 @@ +/* atan.c + * + * Inverse circular tangent + * (arctangent) + * + * + * + * SYNOPSIS: + * + * double x, y, atan(); + * + * y = atan( x ); + * + * + * + * DESCRIPTION: + * + * Returns radian angle between -pi/2 and +pi/2 whose tangent + * is x. + * + * Range reduction is from three intervals into the interval + * from zero to 0.66. The approximant uses a rational + * function of degree 4/5 of the form x + x**3 P(x)/Q(x). + * + * + * + * ACCURACY: + * + * Relative error: + * arithmetic domain # trials peak rms + * DEC -10, 10 50000 2.4e-17 8.3e-18 + * IEEE -10, 10 10^6 1.8e-16 5.0e-17 + * + */ + /* atan2() + * + * Quadrant correct inverse circular tangent + * + * + * + * SYNOPSIS: + * + * double x, y, z, atan2(); + * + * z = atan2( y, x ); + * + * + * + * DESCRIPTION: + * + * Returns radian angle whose tangent is y/x. + * Define compile time symbol ANSIC = 1 for ANSI standard, + * range -PI < z <= +PI, args (y,x); else ANSIC = 0 for range + * 0 to 2PI, args (x,y). + * + * + * + * ACCURACY: + * + * Relative error: + * arithmetic domain # trials peak rms + * IEEE -10, 10 10^6 2.5e-16 6.9e-17 + * See atan.c. + * + */ + +/* atan.c */ + + +/* +Cephes Math Library Release 2.8: June, 2000 +Copyright 1984, 1995, 2000 by Stephen L. Moshier +*/ + + +#include <math.h> + +/* arctan(x) = x + x^3 P(x^2)/Q(x^2) + 0 <= x <= 0.66 + Peak relative error = 2.6e-18 */ +#ifdef UNK +static double P[5] = { +-8.750608600031904122785E-1, +-1.615753718733365076637E1, +-7.500855792314704667340E1, +-1.228866684490136173410E2, +-6.485021904942025371773E1, +}; +static double Q[5] = { +/* 1.000000000000000000000E0, */ + 2.485846490142306297962E1, + 1.650270098316988542046E2, + 4.328810604912902668951E2, + 4.853903996359136964868E2, + 1.945506571482613964425E2, +}; + +/* tan( 3*pi/8 ) */ +static double T3P8 = 2.41421356237309504880; +#endif + +#ifdef DEC +static short P[20] = { +0140140,0001775,0007671,0026242, +0141201,0041242,0155534,0001715, +0141626,0002141,0132100,0011625, +0141765,0142771,0064055,0150453, +0141601,0131517,0164507,0062164, +}; +static short Q[20] = { +/* 0040200,0000000,0000000,0000000, */ +0041306,0157042,0154243,0000742, +0042045,0003352,0016707,0150452, +0042330,0070306,0113425,0170730, +0042362,0130770,0116602,0047520, +0042102,0106367,0156753,0013541, +}; + +/* tan( 3*pi/8 ) = 2.41421356237309504880 */ +static unsigned short T3P8A[] = {040432,0101171,0114774,0167462,}; +#define T3P8 *(double *)T3P8A +#endif + +#ifdef IBMPC +static short P[20] = { +0x2594,0xa1f7,0x007f,0xbfec, +0x807a,0x5b6b,0x2854,0xc030, +0x0273,0x3688,0xc08c,0xc052, +0xba25,0x2d05,0xb8bf,0xc05e, +0xec8e,0xfd28,0x3669,0xc050, +}; +static short Q[20] = { +/* 0x0000,0x0000,0x0000,0x3ff0, */ +0x603c,0x5b14,0xdbc4,0x4038, +0xfa25,0x43b8,0xa0dd,0x4064, +0xbe3b,0xd2e2,0x0e18,0x407b, +0x49ea,0x13b0,0x563f,0x407e, +0x62ec,0xfbbd,0x519e,0x4068, +}; + +/* tan( 3*pi/8 ) = 2.41421356237309504880 */ +static unsigned short T3P8A[] = {0x9de6,0x333f,0x504f,0x4003}; +#define T3P8 *(double *)T3P8A +#endif + +#ifdef MIEEE +static short P[20] = { +0xbfec,0x007f,0xa1f7,0x2594, +0xc030,0x2854,0x5b6b,0x807a, +0xc052,0xc08c,0x3688,0x0273, +0xc05e,0xb8bf,0x2d05,0xba25, +0xc050,0x3669,0xfd28,0xec8e, +}; +static short Q[20] = { +/* 0x3ff0,0x0000,0x0000,0x0000, */ +0x4038,0xdbc4,0x5b14,0x603c, +0x4064,0xa0dd,0x43b8,0xfa25, +0x407b,0x0e18,0xd2e2,0xbe3b, +0x407e,0x563f,0x13b0,0x49ea, +0x4068,0x519e,0xfbbd,0x62ec, +}; + +/* tan( 3*pi/8 ) = 2.41421356237309504880 */ +static unsigned short T3P8A[] = { +0x4003,0x504f,0x333f,0x9de6 +}; +#define T3P8 *(double *)T3P8A +#endif + +#ifdef ANSIPROT +extern double polevl ( double, void *, int ); +extern double p1evl ( double, void *, int ); +extern double atan ( double ); +extern double fabs ( double ); +extern int signbit ( double ); +extern int isnan ( double ); +#else +double polevl(), p1evl(), atan(), fabs(); +//int signbit(), isnan(); +#endif +extern double PI, PIO2, PIO4, INFINITY, NEGZERO, MAXNUM; + +/* pi/2 = PIO2 + MOREBITS. */ +#ifdef DEC +#define MOREBITS 5.721188726109831840122E-18 +#else +#define MOREBITS 6.123233995736765886130E-17 +#endif + + +double atan(x) +double x; +{ +double y, z; +short sign, flag; + +#ifdef MINUSZERO +if( x == 0.0 ) + return(x); +#endif +#ifdef INFINITIES +if(x == INFINITY) + return(PIO2); +if(x == -INFINITY) + return(-PIO2); +#endif +/* make argument positive and save the sign */ +sign = 1; +if( x < 0.0 ) + { + sign = -1; + x = -x; + } +/* range reduction */ +flag = 0; +if( x > T3P8 ) + { + y = PIO2; + flag = 1; + x = -( 1.0/x ); + } +else if( x <= 0.66 ) + { + y = 0.0; + } +else + { + y = PIO4; + flag = 2; + x = (x-1.0)/(x+1.0); + } +z = x * x; +z = z * polevl( z, P, 4 ) / p1evl( z, Q, 5 ); +z = x * z + x; +if( flag == 2 ) + z += 0.5 * MOREBITS; +else if( flag == 1 ) + z += MOREBITS; +y = y + z; +if( sign < 0 ) + y = -y; +return(y); +} + +/* atan2 */ + +#ifdef ANSIC +double atan2( y, x ) +#else +double atan2( x, y ) +#endif +double x, y; +{ +double z, w; +short code; + +code = 0; + +#ifdef NANS +if( isnan(x) ) + return(x); +if( isnan(y) ) + return(y); +#endif +#ifdef MINUSZERO +if( y == 0.0 ) + { + if( signbit(y) ) + { + if( x > 0.0 ) + z = y; + else if( x < 0.0 ) + z = -PI; + else + { + if( signbit(x) ) + z = -PI; + else + z = y; + } + } + else /* y is +0 */ + { + if( x == 0.0 ) + { + if( signbit(x) ) + z = PI; + else + z = 0.0; + } + else if( x > 0.0 ) + z = 0.0; + else + z = PI; + } + return z; + } +if( x == 0.0 ) + { + if( y > 0.0 ) + z = PIO2; + else + z = -PIO2; + return z; + } +#endif /* MINUSZERO */ +#ifdef INFINITIES +if( x == INFINITY ) + { + if( y == INFINITY ) + z = 0.25 * PI; + else if( y == -INFINITY ) + z = -0.25 * PI; + else if( y < 0.0 ) + z = NEGZERO; + else + z = 0.0; + return z; + } +if( x == -INFINITY ) + { + if( y == INFINITY ) + z = 0.75 * PI; + else if( y <= -INFINITY ) + z = -0.75 * PI; + else if( y >= 0.0 ) + z = PI; + else + z = -PI; + return z; + } +if( y == INFINITY ) + return( PIO2 ); +if( y == -INFINITY ) + return( -PIO2 ); +#endif + +if( x < 0.0 ) + code = 2; +if( y < 0.0 ) + code |= 1; + +#ifdef INFINITIES +if( x == 0.0 ) +#else +if( fabs(x) <= (fabs(y) / MAXNUM) ) +#endif + { + if( code & 1 ) + { +#if ANSIC + return( -PIO2 ); +#else + return( 3.0*PIO2 ); +#endif + } + if( y == 0.0 ) + return( 0.0 ); + return( PIO2 ); + } + +if( y == 0.0 ) + { + if( code & 2 ) + return( PI ); + return( 0.0 ); + } + + +switch( code ) + { +#if ANSIC + default: + case 0: + case 1: w = 0.0; break; + case 2: w = PI; break; + case 3: w = -PI; break; +#else + default: + case 0: w = 0.0; break; + case 1: w = 2.0 * PI; break; + case 2: + case 3: w = PI; break; +#endif + } + +z = w + atan( y/x ); +#ifdef MINUSZERO +if( z == 0.0 && y < 0 ) + z = NEGZERO; +#endif +return( z ); +} |