summaryrefslogtreecommitdiff
path: root/libm/double/atan.c
diff options
context:
space:
mode:
Diffstat (limited to 'libm/double/atan.c')
-rw-r--r--libm/double/atan.c393
1 files changed, 393 insertions, 0 deletions
diff --git a/libm/double/atan.c b/libm/double/atan.c
new file mode 100644
index 000000000..f2d50768d
--- /dev/null
+++ b/libm/double/atan.c
@@ -0,0 +1,393 @@
+/* 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 );
+}