/*							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 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
 *    DEC       0, 10        9500       2.8e-17     8.6e-18
 *    DEC       0,MAXNUM     9500       4.2e-17     4.3e-18
 *    IEEE      -10, 10     30000       3.1e-16     7.9e-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     30000       3.6e-16     7.1e-17
 * See atan.c.
 *
 */

/*							atan.c */


/*
Cephes Math Library Release 2.1:  January, 1989
Copyright 1984, 1987, 1989 by Stephen L. Moshier
Direct inquiries to 30 Frost Street, Cambridge, MA 02140
*/


#include "mconf.h"

#ifdef UNK
static double P[] = {
-8.40980878064499716001E-1,
-8.83860837023772394279E0,
-2.18476213081316705724E1,
-1.48307050340438946993E1
};
static double Q[] = {
/* 1.00000000000000000000E0,*/
 1.54974124675307267552E1,
 6.27906555762653017263E1,
 9.22381329856214406485E1,
 4.44921151021319438465E1,
};

/* tan( 3*pi/8 ) */
static double T3P8 = 2.41421356237309504880;

/* tan( pi/8 ) */
static double TP8 = 0.41421356237309504880;
#endif

#ifdef DEC
static short P[] = {
0140127,0045205,0153731,0030027,
0141015,0065360,0116105,0025210,
0141256,0143755,0127056,0105716,
0141155,0045221,0056235,0072437
};
static short Q[] = {
/*0040200,0000000,0000000,0000000,*/
0041167,0172546,0143212,0126224,
0041573,0024641,0116611,0153210,
0041670,0074754,0110422,0127624,
0041461,0173755,0002566,0014374
};

/* tan( 3*pi/8 ) = 2.41421356237309504880 */
static short T3P8A[] = {040432,0101171,0114774,0167462,};
#define T3P8 *(double *)T3P8A

/* tan( pi/8 ) = 0.41421356237309504880 */
static short TP8A[] = {037724,011714,0147747,074622,};
#define TP8 *(double *)TP8A
#endif

#ifdef IBMPC
static short P[] = {
0x2603,0xbafb,0xe950,0xbfea,
0xa551,0x1388,0xad5e,0xc021,
0xd17a,0xb5c5,0xd8fd,0xc035,
0xaea4,0x2b93,0xa952,0xc02d
};
static short Q[] = {
/*0x0000,0x0000,0x0000,0x3ff0,*/
0x5592,0xd8d1,0xfeac,0x402e,
0x3ad1,0x33b1,0x6534,0x404f,
0x55f2,0x9222,0x0f3d,0x4057,
0xc31f,0xa0ae,0x3efd,0x4046
};

/* tan( 3*pi/8 ) = 2.41421356237309504880 */
static short T3P8A[] = {0x9de6,0x333f,0x504f,0x4003};
#define T3P8 *(double *)T3P8A

/* tan( pi/8 ) = 0.41421356237309504880 */
static short TP8A[] = {0xef32,0x99fc,0x8279,0x3fda};
#define TP8 *(double *)TP8A
#endif

#ifdef MIEEE
static short P[] = {
0xbfea,0xe950,0xbafb,0x2603,
0xc021,0xad5e,0x1388,0xa551,
0xc035,0xd8fd,0xb5c5,0xd17a,
0xc02d,0xa952,0x2b93,0xaea4
};
static short Q[] = {
0x402e,0xfeac,0xd8d1,0x5592,
0x404f,0x6534,0x33b1,0x3ad1,
0x4057,0x0f3d,0x9222,0x55f2,
0x4046,0x3efd,0xa0ae,0xc31f,
};

/* tan( 3*pi/8 ) = 2.41421356237309504880 */
static short T3P8A[] = {
0x4003,0x504f,0x333f,0x9de6
};
#define T3P8 *(double *)T3P8A

/* tan( pi/8 ) = 0.41421356237309504880 */
static short TP8A[] = {
0x3fda,0x8279,0x99fc,0xef32
};
#define TP8 *(double *)TP8A
#endif

double atan(x)
double x;
{
extern double PIO2, PIO4;
double y,z,p,q;
double polevl(), p1evl();
short sign;

/* make argument positive and save the sign */
sign = 1;
if( x < 0.0 )
	{
	sign = -1;
	x = -x;
	}

/* range reduction */
if( x > T3P8 )
	{
	y = PIO2;
	x = -( 1.0/x );
	}

else if( x > TP8 )
	{
	y = PIO4;
	x = (x-1.0)/(x+1.0);
	}
else
	y = 0.0;

/* rational form in x**2 */
z = x * x;
y = y + ( polevl( z, P, 3 ) / p1evl( z, Q, 4 ) ) * z * x + x;

if( sign < 0 )
	y = -y;

return(y);
}

/*							atan2	*/

#ifndef GENERIC
static char fname[] = {"atan2"};
#endif /* GENERIC */

extern double PI, PIO2;

#if ANSIC
double atan2( y, x )
#else
double atan2( x, y )
#endif
double x, y;
{
double z, w;
short code;
double atan();


code = 0;

if( x < 0.0 )
	code = 2;
if( y < 0.0 )
	code |= 1;

if( x == 0.0 )
	{
	if( code & 1 )
		{
#if ANSIC
		return( -PIO2 );
#else
		return( 3.0*PIO2 );
#endif
		}
	if( y == 0.0 )
#ifndef GENERIC
		mtherr (fname, DOMAIN);
#endif /* GENERIC */
		return( 0.0 );
	return( PIO2 );
	}

if( y == 0.0 )
	{
	if( code & 2 )
		return( PI );
	return( 0.0 );
	}


switch( code )
	{
#if ANSIC
	case 0:
	case 1: w = 0.0; break;
	case 2: w = PI; break;
	case 3: w = -PI; break;
#else
	case 0: w = 0.0; break;
	case 1: w = 2.0 * PI; break;
	case 2:
	case 3: w = PI; break;
#endif
	}

z = atan( y/x );

return( w + z );
}