atanll.c

128位长双精度型数字运算包

/*							atanl.c
 *
 *	Inverse circular tangent, 128-bit long double precision
 *      (arctangent)
 *
 *
 *
 * SYNOPSIS:
 *
 * long double x, y, atanl();
 *
 * y = atanl( 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
 *    IEEE      -10, 10    100,000      2.6e-34     6.5e-35
 *
 */
/*							atan2l()
 *
 *	Quadrant correct inverse circular tangent,
 *	long double precision
 *
 *
 *
 * SYNOPSIS:
 *
 * long double x, y, z, atan2l();
 *
 * z = atan2l( 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    100,000      3.2e-34      5.9e-35
 * See atan.c.
 *
 */

/*							atan.c */


/*
Cephes Math Library Release 2.2:  December, 1990
Copyright 1984, 1990 by Stephen L. Moshier
Direct inquiries to 30 Frost Street, Cambridge, MA 02140
*/


#include "mconf.h"

/* arctan(x) = x + x^3 P(x^2)
 * Theoretical peak relative error = 3.0e-36
 * relative peak error spread = 6.6e-8
 */
static long double P[9] = {
-6.635810778635296712545011270011752799963E-4L,
-8.768423468036849091777415076702113400070E-1L,
-2.548067867495502632615671450650071218995E1L,
-2.497759878476618348858065206895055957104E2L,
-1.148164399808514330375280133523543970854E3L,
-2.792272753241044941703278827346430350236E3L,
-3.696264445691821235400930243493001671932E3L,
-2.514829758941713674909996882101723647996E3L,
-6.880597774405940432145577545328795037141E2L
};
static long double Q[8] = {
/* 1.000000000000000000000000000000000000000E0L, */
 3.566239794444800849656497338030115886153E1L,
 4.308348370818927353321556740027020068897E2L,
 2.494680540950601626662048893678584497900E3L,
 7.928572347062145288093560392463784743935E3L,
 1.458510242529987155225086911411015961174E4L,
 1.547394317752562611786521896296215170819E4L,
 8.782996876218210302516194604424986107121E3L,
 2.064179332321782129643673263598686441900E3L
};

/* tan( 3*pi/8 ) */
static long double T3P8 = 2.414213562373095048801688724209698078569672L;

/* tan( pi/8 ) */
static long double TP8 = 0.414213562373095048801688724209698078569672L;




long double atanl(x)
long double x;
{
extern long double PIO2L, PIO4L;
long double y, z;
long double polevll(), p1evll();
short sign;

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

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

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

/* rational form in x**2 */
z = x * x;
y = y + ( polevll( z, P, 8 ) / p1evll( z, Q, 8 ) ) * z * x + x;

if( sign < 0 )
	y = -y;

return(y);
}

/*							atan2	*/


extern long double PIL, PIO2L;

#if ANSIC
long double atan2l( y, x )
#else
long double atan2l( x, y )
#endif
long double x, y;
{
long double z, w;
short code;
long double atanl();


code = 0;
w = 0.0L;

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

if( x == 0.0L )
	{
	if( code & 1 )
		{
#if ANSIC
		return( -PIO2L );
#else
		return( 3.0L*PIO2L );
#endif
		}
	if( y == 0.0L )
		return( 0.0L );
	return( PIO2L );
	}

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


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

z = atanl( y/x );

return( w + z );
}