diff options
author | Eric Andersen <andersen@codepoet.org> | 2001-11-22 14:04:29 +0000 |
---|---|---|
committer | Eric Andersen <andersen@codepoet.org> | 2001-11-22 14:04:29 +0000 |
commit | 7ce331c01ce6eb7b3f5c715a38a24359da9c6ee2 (patch) | |
tree | 3a7e8476e868ae15f4da1b7ce26b2db6f434468c /libm/ldouble/atanl.c | |
parent | c117dd5fb183afb1a4790a6f6110d88704be6bf8 (diff) |
Totally rework the math library, this time based on the MacOs X
math library (which is itself based on the math lib from FreeBSD).
-Erik
Diffstat (limited to 'libm/ldouble/atanl.c')
-rw-r--r-- | libm/ldouble/atanl.c | 376 |
1 files changed, 0 insertions, 376 deletions
diff --git a/libm/ldouble/atanl.c b/libm/ldouble/atanl.c deleted file mode 100644 index 9e6d9af3c..000000000 --- a/libm/ldouble/atanl.c +++ /dev/null @@ -1,376 +0,0 @@ -/* 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 ); -} |