summaryrefslogtreecommitdiff
path: root/libm/double/atan.c
diff options
context:
space:
mode:
authorEric Andersen <andersen@codepoet.org>2001-11-22 14:04:29 +0000
committerEric Andersen <andersen@codepoet.org>2001-11-22 14:04:29 +0000
commit7ce331c01ce6eb7b3f5c715a38a24359da9c6ee2 (patch)
tree3a7e8476e868ae15f4da1b7ce26b2db6f434468c /libm/double/atan.c
parentc117dd5fb183afb1a4790a6f6110d88704be6bf8 (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/double/atan.c')
-rw-r--r--libm/double/atan.c393
1 files changed, 0 insertions, 393 deletions
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 <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 );
-}