From 7ce331c01ce6eb7b3f5c715a38a24359da9c6ee2 Mon Sep 17 00:00:00 2001 From: Eric Andersen Date: Thu, 22 Nov 2001 14:04:29 +0000 Subject: 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 --- libm/double/atan.c | 393 ----------------------------------------------------- 1 file changed, 393 deletions(-) delete mode 100644 libm/double/atan.c (limited to 'libm/double/atan.c') diff --git a/libm/double/atan.c b/libm/double/atan.c deleted file mode 100644 index f2d50768d..000000000 --- a/libm/double/atan.c +++ /dev/null @@ -1,393 +0,0 @@ -/* 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 - -/* 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 ); -} -- cgit v1.2.3