/* 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 "mconf.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); }