diff options
Diffstat (limited to 'libm/ldouble/atanl.c')
| -rw-r--r-- | libm/ldouble/atanl.c | 376 |
1 files changed, 376 insertions, 0 deletions
diff --git a/libm/ldouble/atanl.c b/libm/ldouble/atanl.c new file mode 100644 index 000000000..9e6d9af3c --- /dev/null +++ b/libm/ldouble/atanl.c @@ -0,0 +1,376 @@ +/* atanl.c + * + * Inverse circular tangent, long double precision + * (arctangent) + * + * + * + * SYNOPSIS: + * + * long double x, y, atanl(); + * + * y = atanl( x ); + * + * + * + * DESCRIPTION: + * + * Returns radian angle between -pi/2 and +pi/2 whose tangent + * is x. + * + * Range reduction is from four intervals into the interval + * from zero to tan( pi/8 ). The approximant uses a rational + * function of degree 3/4 of the form x + x**3 P(x)/Q(x). + * + * + * + * ACCURACY: + * + * Relative error: + * arithmetic domain # trials peak rms + * IEEE -10, 10 150000 1.3e-19 3.0e-20 + * + */ + /* atan2l() + * + * Quadrant correct inverse circular tangent, + * long double precision + * + * + * + * SYNOPSIS: + * + * long double x, y, z, atan2l(); + * + * z = atan2l( 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 60000 1.7e-19 3.2e-20 + * See atan.c. + * + */ + +/* atan.c */ + + +/* +Cephes Math Library Release 2.7: May, 1998 +Copyright 1984, 1990, 1998 by Stephen L. Moshier +*/ + + +#include <math.h> + +#ifdef UNK +static long double P[] = { +-8.6863818178092187535440E-1L, +-1.4683508633175792446076E1L, +-6.3976888655834347413154E1L, +-9.9988763777265819915721E1L, +-5.0894116899623603312185E1L, +}; +static long double Q[] = { +/* 1.00000000000000000000E0L,*/ + 2.2981886733594175366172E1L, + 1.4399096122250781605352E2L, + 3.6144079386152023162701E2L, + 3.9157570175111990631099E2L, + 1.5268235069887081006606E2L, +}; + +/* tan( 3*pi/8 ) */ +static long double T3P8 = 2.41421356237309504880169L; + +/* tan( pi/8 ) */ +static long double TP8 = 4.1421356237309504880169e-1L; +#endif + + +#ifdef IBMPC +static unsigned short P[] = { +0x8ece,0xce53,0x1266,0xde5f,0xbffe, XPD +0x07e6,0xa061,0xa6bf,0xeaef,0xc002, XPD +0x53ee,0xf291,0x557f,0xffe8,0xc004, XPD +0xf9d6,0xeda6,0x3f3e,0xc7fa,0xc005, XPD +0xb6c3,0x6abc,0x9361,0xcb93,0xc004, XPD +}; +static unsigned short Q[] = { +/*0x0000,0x0000,0x0000,0x8000,0x3fff,*/ +0x54d4,0x894e,0xe76e,0xb7da,0x4003, XPD +0x76b9,0x7a46,0xafa2,0x8ffd,0x4006, XPD +0xe3a9,0xe9c0,0x6bee,0xb4b8,0x4007, XPD +0xabc1,0x50a7,0xb098,0xc3c9,0x4007, XPD +0x891c,0x100d,0xae89,0x98ae,0x4006, XPD +}; + +/* tan( 3*pi/8 ) = 2.41421356237309504880 */ +static unsigned short T3P8A[] = {0x3242,0xfcef,0x7999,0x9a82,0x4000, XPD}; +#define T3P8 *(long double *)T3P8A + +/* tan( pi/8 ) = 0.41421356237309504880 */ +static unsigned short TP8A[] = {0x9211,0xe779,0xcccf,0xd413,0x3ffd, XPD}; +#define TP8 *(long double *)TP8A +#endif + +#ifdef MIEEE +static unsigned long P[] = { +0xbffe0000,0xde5f1266,0xce538ece, +0xc0020000,0xeaefa6bf,0xa06107e6, +0xc0040000,0xffe8557f,0xf29153ee, +0xc0050000,0xc7fa3f3e,0xeda6f9d6, +0xc0040000,0xcb939361,0x6abcb6c3, +}; +static unsigned long Q[] = { +/*0x3fff0000,0x80000000,0x00000000,*/ +0x40030000,0xb7dae76e,0x894e54d4, +0x40060000,0x8ffdafa2,0x7a4676b9, +0x40070000,0xb4b86bee,0xe9c0e3a9, +0x40070000,0xc3c9b098,0x50a7abc1, +0x40060000,0x98aeae89,0x100d891c, +}; + +/* tan( 3*pi/8 ) = 2.41421356237309504880 */ +static long T3P8A[] = {0x40000000,0x9a827999,0xfcef3242}; +#define T3P8 *(long double *)T3P8A + +/* tan( pi/8 ) = 0.41421356237309504880 */ +static long TP8A[] = {0x3ffd0000,0xd413cccf,0xe7799211}; +#define TP8 *(long double *)TP8A +#endif + +#ifdef ANSIPROT +extern long double polevll ( long double, void *, int ); +extern long double p1evll ( long double, void *, int ); +extern long double fabsl ( long double ); +extern int signbitl ( long double ); +extern int isnanl ( long double ); +long double atanl ( long double ); +#else +long double polevll(), p1evll(), fabsl(), signbitl(), isnanl(); +long double atanl(); +#endif +#ifdef INFINITIES +extern long double INFINITYL; +#endif +#ifdef NANS +extern long double NANL; +#endif +#ifdef MINUSZERO +extern long double NEGZEROL; +#endif + +long double atanl(x) +long double x; +{ +extern long double PIO2L, PIO4L; +long double y, z; +short sign; + +#ifdef MINUSZERO +if( x == 0.0L ) + return(x); +#endif +#ifdef INFINITIES +if( x == INFINITYL ) + return( PIO2L ); +if( x == -INFINITYL ) + return( -PIO2L ); +#endif +/* make argument positive and save the sign */ +sign = 1; +if( x < 0.0L ) + { + sign = -1; + x = -x; + } + +/* range reduction */ +if( x > T3P8 ) + { + y = PIO2L; + x = -( 1.0L/x ); + } + +else if( x > TP8 ) + { + y = PIO4L; + x = (x-1.0L)/(x+1.0L); + } +else + y = 0.0L; + +/* rational form in x**2 */ +z = x * x; +y = y + ( polevll( z, P, 4 ) / p1evll( z, Q, 5 ) ) * z * x + x; + +if( sign < 0 ) + y = -y; + +return(y); +} + +/* atan2 */ + + +extern long double PIL, PIO2L, MAXNUML; + +#if ANSIC +long double atan2l( y, x ) +#else +long double atan2l( x, y ) +#endif +long double x, y; +{ +long double z, w; +short code; + +code = 0; + +if( x < 0.0L ) + code = 2; +if( y < 0.0L ) + code |= 1; + +#ifdef NANS +if( isnanl(x) ) + return(x); +if( isnanl(y) ) + return(y); +#endif +#ifdef MINUSZERO +if( y == 0.0L ) + { + if( signbitl(y) ) + { + if( x > 0.0L ) + z = y; + else if( x < 0.0L ) + z = -PIL; + else + { + if( signbitl(x) ) + z = -PIL; + else + z = y; + } + } + else /* y is +0 */ + { + if( x == 0.0L ) + { + if( signbitl(x) ) + z = PIL; + else + z = 0.0L; + } + else if( x > 0.0L ) + z = 0.0L; + else + z = PIL; + } + return z; + } +if( x == 0.0L ) + { + if( y > 0.0L ) + z = PIO2L; + else + z = -PIO2L; + return z; + } +#endif /* MINUSZERO */ +#ifdef INFINITIES +if( x == INFINITYL ) + { + if( y == INFINITYL ) + z = 0.25L * PIL; + else if( y == -INFINITYL ) + z = -0.25L * PIL; + else if( y < 0.0L ) + z = NEGZEROL; + else + z = 0.0L; + return z; + } +if( x == -INFINITYL ) + { + if( y == INFINITYL ) + z = 0.75L * PIL; + else if( y == -INFINITYL ) + z = -0.75L * PIL; + else if( y >= 0.0L ) + z = PIL; + else + z = -PIL; + return z; + } +if( y == INFINITYL ) + return( PIO2L ); +if( y == -INFINITYL ) + return( -PIO2L ); +#endif /* INFINITIES */ + +#ifdef INFINITIES +if( x == 0.0L ) +#else +if( fabsl(x) <= (fabsl(y) / MAXNUML) ) +#endif + { + if( code & 1 ) + { +#if ANSIC + return( -PIO2L ); +#else + return( 3.0L*PIO2L ); +#endif + } + if( y == 0.0L ) + return( 0.0L ); + return( PIO2L ); + } + +if( y == 0.0L ) + { + if( code & 2 ) + return( PIL ); + return( 0.0L ); + } + + +switch( code ) + { + default: +#if ANSIC + case 0: + case 1: w = 0.0L; break; + case 2: w = PIL; break; + case 3: w = -PIL; break; +#else + case 0: w = 0.0L; break; + case 1: w = 2.0L * PIL; break; + case 2: + case 3: w = PIL; break; +#endif + } + +z = w + atanl( y/x ); +#ifdef MINUSZERO +if( z == 0.0L && y < 0.0L ) + z = NEGZEROL; +#endif +return( z ); +} |