Quelle  atanl.c

/* atanl.c
 *
 * Inverse circular tangent, 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    150000       1.3e-19     3.0e-20
 *
 */
/* 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     60000       1.7e-19     3.2e-20
 * See atan.c.
 *
 */

/* atan.c */


/*
Cephes Math Library Release 2.7:  May, 1998
Copyright 1984, 1990, 1998 by Stephen L. Moshier
*/


#include "mconf.h"

#ifdef UNK
static long double P[] = {
-8.6863818178092187535440E-1L,
-1.4683508633175792446076E1L,
-6.3976888655834347413154E1L,
-9.9988763777265819915721E1L,
-5.0894116899623603312185E1L,
};
static long double Q[] = {
/* 1.00000000000000000000E0L,*/
 2.2981886733594175366172E1L,
 1.4399096122250781605352E2L,
 3.6144079386152023162701E2L,
 3.9157570175111990631099E2L,
 1.5268235069887081006606E2L,
};

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

/* tan( pi/8 ) */
static long double TP8 =  4.1421356237309504880169e-1L;
#endif


#ifdef IBMPC
static unsigned short P[] = {
0x8ece,0xce53,0x1266,0xde5f,0xbffe, XPD
0x07e6,0xa061,0xa6bf,0xeaef,0xc002, XPD
0x53ee,0xf291,0x557f,0xffe8,0xc004, XPD
0xf9d6,0xeda6,0x3f3e,0xc7fa,0xc005, XPD
0xb6c3,0x6abc,0x9361,0xcb93,0xc004, XPD
};
static unsigned short Q[] = {
/*0x0000,0x0000,0x0000,0x8000,0x3fff,*/
0x54d4,0x894e,0xe76e,0xb7da,0x4003, XPD
0x76b9,0x7a46,0xafa2,0x8ffd,0x4006, XPD
0xe3a9,0xe9c0,0x6bee,0xb4b8,0x4007, XPD
0xabc1,0x50a7,0xb098,0xc3c9,0x4007, XPD
0x891c,0x100d,0xae89,0x98ae,0x4006, XPD
};

/* tan( 3*pi/8 ) = 2.41421356237309504880 */
static unsigned short T3P8A[] = {0x3242,0xfcef,0x7999,0x9a82,0x4000, XPD};
#define T3P8 *(long double *)T3P8A

/* tan( pi/8 ) = 0.41421356237309504880 */
static unsigned short TP8A[] = {0x9211,0xe779,0xcccf,0xd413,0x3ffd, XPD};
#define TP8 *(long double *)TP8A
#endif

#ifdef MIEEE
static unsigned long P[] = {
0xbffe0000,0xde5f1266,0xce538ece,
0xc0020000,0xeaefa6bf,0xa06107e6,
0xc0040000,0xffe8557f,0xf29153ee,
0xc0050000,0xc7fa3f3e,0xeda6f9d6,
0xc0040000,0xcb939361,0x6abcb6c3,
};
static unsigned long Q[] = {
/*0x3fff0000,0x80000000,0x00000000,*/
0x40030000,0xb7dae76e,0x894e54d4,
0x40060000,0x8ffdafa2,0x7a4676b9,
0x40070000,0xb4b86bee,0xe9c0e3a9,
0x40070000,0xc3c9b098,0x50a7abc1,
0x40060000,0x98aeae89,0x100d891c,
};

/* tan( 3*pi/8 ) = 2.41421356237309504880 */
static long T3P8A[] = {0x40000000,0x9a827999,0xfcef3242};
#define T3P8 *(long double *)T3P8A

/* tan( pi/8 ) = 0.41421356237309504880 */
static long TP8A[] = {0x3ffd0000,0xd413cccf,0xe7799211};
#define TP8 *(long double *)TP8A
#endif

#ifdef ANSIPROT
extern long double polevll ( long double, void *, int );
extern long double p1evll ( long double, void *, int );
extern long double fabsl ( long double );
extern int signbitl ( long double );
extern int isnanl ( long double );
long double atanl ( long double );
#else
long double polevll(), p1evll(), fabsl(), signbitl(), isnanl();
long double atanl();
#endif
#ifdef INFINITIES
extern long double INFINITYL;
#endif
#ifdef NANS
extern long double NANL;
#endif
#ifdef MINUSZERO
extern long double NEGZEROL;
#endif

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

#ifdef MINUSZERO
if( x == 0.0L )
 return(x);
#endif
#ifdef INFINITIES
if( x == INFINITYL )
 return( PIO2L );
if( x == -INFINITYL )
 return( -PIO2L );
#endif
/* 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, 4 ) / p1evll( z, Q, 5 ) ) * z * x + x;

if( sign < 0 )
 y = -y;

return(y);
}

/* atan2 */


extern long double PIL, PIO2L, MAXNUML;

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

code = 0;

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

#ifdef NANS
if( isnanl(x) )
 return(x);
if( isnanl(y) )
 return(y);
#endif
#ifdef MINUSZERO
if( y == 0.0L )
 {
 if( signbitl(y) )
  {
  if( x > 0.0L )
   z = y;
  else if( x < 0.0L )
   z = -PIL;
  else
   {
   if( signbitl(x) )
    z = -PIL;
   else
    z = y;
   }
  }
 else /* y is +0 */
  {
  if( x == 0.0L )
   {
   if( signbitl(x) )
    z = PIL;
   else
    z = 0.0L;
   }
  else if( x > 0.0L )
   z = 0.0L;
  else
   z = PIL;
  }
 return z;
 }
if( x == 0.0L )
 {
 if( y > 0.0L )
  z = PIO2L;
 else
  z = -PIO2L;
 return z;
 }
#endif /* MINUSZERO */
#ifdef INFINITIES
if( x == INFINITYL )
 {
 if( y == INFINITYL )
  z = 0.25L * PIL;
 else if( y == -INFINITYL )
  z = -0.25L * PIL;
 else if( y < 0.0L )
  z = NEGZEROL;
 else
  z = 0.0L;
 return z;
 }
if( x == -INFINITYL )
 {
 if( y == INFINITYL )
  z = 0.75L * PIL;
 else if( y == -INFINITYL )
  z = -0.75L * PIL;
 else if( y >= 0.0L )
  z = PIL;
 else
  z = -PIL;
 return z;
 }
if( y == INFINITYL )
 return( PIO2L );
if( y == -INFINITYL )
 return( -PIO2L );
#endif /* INFINITIES */

#ifdef INFINITIES
if( x == 0.0L )
#else
if( fabsl(x) <= (fabsl(y) / MAXNUML) )
#endif
 {
 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 )
 {
 default:
#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 = w + atanl( y/x );
#ifdef MINUSZERO
if( z == 0.0L && y < 0.0L )
 z = NEGZEROL;
#endif
return( z );
}