/* ====================================================================
 * ====================================================================
 *
 * Module: asin.c
 * $Revision: 1.1.1.2 $
 * $Date: 2002/10/29 08:06:00 $
 * $Author: blythe $
 * $Source: /root/leakn64/depot/rf/sw/n64os20l/apps/spin/asin.c,v $
 *
 * Revision history:
 *  09-Jun-93 - Original Version
 *
 * Description:	source code for asin function
 *
 * ====================================================================
 * ====================================================================
 */


#ifdef _CALL_MATHERR
#include <stdio.h>
#include <math.h>
#include <errno.h>
#endif

#include "libm.h"
#include "gu.h"

#pragma weak asin = __asin
#define	asin __asin

#ifndef ABS
#define ABS(x)		((x) < 0 ? -(x) : (x))
#endif

/* coefficients for rational approximation of asin on +/- 0.5    */

static const du	P[] =
{
{0x00000000,	0x00000000},
{0x400cf4d7,	0x1166375d},
{0xc019dde2,	0xfd680d32},
{0x400c15ee,	0x4cc68a6a},
{0xbfe20f7e,	0xdc1c40fe},
{0x3f780cd5,	0x52bc78fd},
};

static const du	Q[] =
{
{0x4035b7a1,	0x4d0ca925},
{0xc0484954,	0xef63fb64},
{0x40428d6e,	0x183c02d2},
{0xc026104e,	0x748fbc12},
{0x3ff00000,	0x00000000},
};

/* coefficients for rational approximation of asin on +/- sqrt(2 - sqrt(3))/2 */

static const du	P2[] =
{
{0x00000000,	0x00000000},
{0xc0097a02,	0x2e7a0e13},
{0x4009aeb0,	0x0736a7c7},
{0xbfe37279,	0x9e195a2e},
};

static const du	Q2[] =
{
{0xc0331b81,	0xa2db8f8a},
{0x403bdc31,	0x8eb7aeb5},
{0xc0262173,	0xdc9ece8a},
{0x3ff00000,	0x00000000},
};

static const du	one =
{0x3ff00000,	0x00000000};

static const du	m_one =
{0xbff00000,	0x00000000};

static const du	piby4 =
{0x3fe921fb,	0x54442d18};

static const du	piby2 =
{0x3ff921fb,	0x54442d18};

static const du	m_piby4 =
{0xbfe921fb,	0x54442d18};

static const du	m_piby2 =
{0xbff921fb,	0x54442d18};

static const du	root3by2 =
{0x3febb67a,	0xe8584caa};

#ifdef _CALL_MATHERR
static struct exception	exstruct;
#endif


/* ====================================================================
 *
 * FunctionName		asin
 *
 * Description		computes arcsine of arg
 *
 * ====================================================================
 */

double
  asin( double x )
{
double	xsq, num, denom, result;
double	y, ysq;
double	q, absx;
int	ix, xpt;

	/* extract exponent of x for some quick screening */

	ix = *(int *)&x;
	xpt = (ix >> 20);
	xpt &= 0x7ff;

	if ( xpt < 0x3fe )
	{
		/* |x| < 0.5    */

		if ( xpt >= 0x3e3 )
		{
			/* |x| >= 2^(-28)   */

			xsq = x*x;
			num = ((((P[5].d*xsq + P[4].d)*xsq + P[3].d)*xsq + P[2].d)*xsq +
				P[1].d);

			denom = ((((xsq + Q[3].d)*xsq + Q[2].d)*xsq + Q[1].d)*xsq +
				Q[0].d);

			result = x + x*(xsq*num)/denom;

			return ( result );
		}

		return ( x );
	}
	
	if ( xpt < 0x3ff )
	{
		/*  |x| < 1.0   */

		absx = ABS(x);

		if ( absx < root3by2.d )
		{
			/* |x| < sqrt(3)/2    */

			xsq = x*x;
			y = xsq + xsq - one.d;
			xsq = y*y;
	
			num = ((((P[5].d*xsq + P[4].d)*xsq + P[3].d)*xsq + P[2].d)*xsq +
				P[1].d);
	
			denom = ((((xsq + Q[3].d)*xsq + Q[2].d)*xsq + Q[1].d)*xsq +
				Q[0].d);
	
			result = y + y*(xsq*num)/denom;
	
			if ( ix > 0 )
				return ( piby4.d + 0.5*result );
			else
				return ( m_piby4.d - 0.5*result );
		}
	
		ysq = 0.5*(one.d - absx);
		y = (double)sqrtf((float)ysq);

		num = ((P2[3].d*ysq + P2[2].d)*ysq + P2[1].d);
		denom = (((ysq + Q2[2].d)*ysq + Q2[1].d)*ysq + Q2[0].d);

		q = y + y*(ysq*num)/denom;

		if ( ix > 0 )
			return ( piby2.d - (q + q) );
		else
			return ( m_piby2.d + (q + q) );
	}

	if ( x != x )
	{
		/* x is a NaN; return a quiet NaN */

#ifdef _CALL_MATHERR

		exstruct.type = DOMAIN;
		exstruct.name = "asin";
		exstruct.arg1 = x;
		exstruct.retval = __libm_qnan_d;

		if ( matherr( &exstruct ) == 0 )
		{
			fprintf(stderr, "domain error in asin\n");
			*__errnoaddr = EDOM;
		}

		return ( exstruct.retval );
#else

#ifdef _IP_NAN_SETS_ERRNO

		*__errnoaddr = EDOM;
#endif

		/* return ( __libm_qnan_d ); */
		return ( m_piby2.d + (q + q) );
#endif
	}

	if ( x == one.d )
		return ( piby2.d );

	if ( x == m_one.d )
		return ( m_piby2.d );

	/*  |x| > 1.0; return a NaN */

#ifdef _CALL_MATHERR

	exstruct.type = DOMAIN;
	exstruct.name = "asin";
	exstruct.arg1 = x;
	exstruct.retval = __libm_qnan_d;

	if ( matherr( &exstruct ) == 0 )
	{
		fprintf(stderr, "domain error in asin\n");
		*__errnoaddr = EDOM;
	}

	return ( exstruct.retval );
#else

#ifdef _IP_NAN_SETS_ERRNO

	*__errnoaddr = EDOM;
#endif

	/* return ( __libm_qnan_d ); */
	return ( m_piby2.d );
#endif
}