/*							asin.c
 *
 *	Inverse circular sine
 *
 *
 *
 * SYNOPSIS:
 *
 * double x, y, asin();
 *
 * y = asin( x );
 *
 *
 *
 * DESCRIPTION:
 *
 * Returns radian angle between -pi/2 and +pi/2 whose sine is x.
 *
 * A rational function of the form x + x**3 P(x**2)/Q(x**2)
 * is used for |x| in the interval [0, 0.5].  If |x| > 0.5 it is
 * transformed by the identity
 *
 *    asin(x) = pi/2 - 2 asin( sqrt( (1-x)/2 ) ).
 *
 *
 * ACCURACY:
 *
 *                      Relative error:
 * arithmetic   domain     # trials      peak         rms
 *    DEC      -1, 1        20000       5.7e-17     1.2e-17
 *    IEEE     -1, 1        30000       4.8e-16     8.4e-17
 *
 *
 * ERROR MESSAGES:
 *
 *   message         condition      value returned
 * asin domain        |x| > 1           0.0
 *
 */
/*							acos()
 *
 *	Inverse circular cosine
 *
 *
 *
 * SYNOPSIS:
 *
 * double x, y, acos();
 *
 * y = acos( x );
 *
 *
 *
 * DESCRIPTION:
 *
 * Returns radian angle between -pi/2 and +pi/2 whose cosine
 * is x.
 *
 * Analytically, acos(x) = pi/2 - asin(x).  However if |x| is
 * near 1, there is cancellation error in subtracting asin(x)
 * from pi/2.  Hence if x < -0.5,
 *
 *    acos(x) =	 pi - 2.0 * asin( sqrt((1+x)/2) );
 *
 * or if x > +0.5,
 *
 *    acos(x) =	 2.0 * asin(  sqrt((1-x)/2) ).
 *
 *
 * ACCURACY:
 *
 *                      Relative error:
 * arithmetic   domain     # trials      peak         rms
 *    DEC       -1, 1       13000       3.2e-17     7.8e-18
 *    IEEE      -1, 1       30000       2.7e-16     7.1e-17
 *
 *
 * ERROR MESSAGES:
 *
 *   message         condition      value returned
 * asin domain        |x| > 1           0.0
 */

/*							asin.c	*/

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

#include "mconf.h"

#ifdef UNK
static double P[] = {
-6.96822599948686174217E-1,
 1.01598286089872099722E1,
-3.97340771391578757294E1,
 5.72912144709846496134E1,
-2.74148200465925708020E1
};
static double Q[] = {
/* 1.00000000000000000000E0,*/
-2.38368245005177488242E1,
 1.51095072703128995631E2,
-3.82340216045978957023E2,
 4.17767300951716199422E2,
-1.64488920279555473283E2
};
#endif

#ifdef DEC
static short P[] = {
0140062,0061367,0042744,0127464,
0041042,0107250,0070611,0005470,
0141436,0167661,0165345,0131200,
0041545,0025064,0020124,0000411,
0141333,0050615,0026056,0134351
};
static short Q[] = {
/*0040200,0000000,0000000,0000000,*/
0141276,0130721,0005461,0134336,
0042027,0014126,0127506,0127155,
0142277,0025614,0031413,0103353,
0042320,0161066,0165346,0163707,
0142044,0076451,0160443,0005274
};
#endif

#ifdef IBMPC
static short P[] = {
0x95e7,0xe8bc,0x4c5e,0xbfe6,
0x2167,0x0e31,0x51d5,0x4024,
0xb650,0x3d5c,0xddf6,0xc043,
0x8021,0x840a,0xa546,0x404c,
0xd71d,0xa585,0x6a31,0xc03b
};
static short Q[] = {
/*0x0000,0x0000,0x0000,0x3ff0,*/
0x371c,0x2166,0xd63a,0xc037,
0xd5ce,0xd5e8,0xe30a,0x4062,
0x70dd,0x8661,0xe571,0xc077,
0xdcf9,0xdd5c,0x1c46,0x407a,
0x6158,0x3c24,0x8fa5,0xc064
};
#endif

#ifdef MIEEE
static short P[] = {
0xbfe6,0x4c5e,0xe8bc,0x95e7,
0x4024,0x51d5,0x0e31,0x2167,
0xc043,0xddf6,0x3d5c,0xb650,
0x404c,0xa546,0x840a,0x8021,
0xc03b,0x6a31,0xa585,0xd71d
};
static short Q[] = {
0xc037,0xd63a,0x2166,0x371c,
0x4062,0xe30a,0xd5e8,0xd5ce,
0xc077,0xe571,0x8661,0x70dd,
0x407a,0x1c46,0xdd5c,0xdcf9,
0xc064,0x8fa5,0x3c24,0x6158
};
#endif

double asin(x)
double x;
{
double a, p, z, zz;
double sqrt(), polevl(), p1evl();
short sign, flag;
extern double PIO2;

if( x > 0 )
	{
	sign = 1;
	a = x;
	}
else
	{
	sign = -1;
	a = -x;
	}

if( a > 1.0 )
	{
	mtherr( "asin", DOMAIN );
	return( 0.0 );
	}

if( a < 1.0e-7 )
	{
	z = a;
	goto done;
	}

if( a > 0.5 )
	{
	zz = 0.5 -a;
	zz = (zz + 0.5)/2.0;
	z = sqrt( zz );
	flag = 1;
	}
else
	{
	z = a;
	zz = z * z;
	flag = 0;
	}

p = zz * polevl( zz, P, 4)/p1evl( zz, Q, 5);
z = z * p + z;
if( flag != 0 )
	{
	z = z + z;
	z = PIO2 - z;
	}
done:
if( sign < 0 )
	z = -z;
return(z);
}


extern double PIO2, PI;

double acos(x)
double x;
{
double asin(), sqrt();

if( x < -1.0 )
	goto domerr;

if( x < -0.5) 
	return( PI - 2.0 * asin( sqrt(0.5*(1.0+x)) ) );

if( x > 1.0 )
	{
domerr:	mtherr( "acos", DOMAIN );
	return( 0.0 );
	}

if( x > 0.5 )
	return( 2.0 * asin(  sqrt(0.5*(1.0-x) ) ) );

return( PIO2 - asin(x) );
}