/*
 *				a r s i c o . c
 */

/*)LIBRARY
*/

#ifdef	DOCUMENTATION
title	arsico	arc sin and arc cosine functions
index	arc sin function
index	arc cosine function
usage
	.s
	double x, f, asin();
	.br
	double g, acos();
	.br
	f = asin(x);
	.br
	g = acos(x);
	.s
description
	.s
	asin returns arcsin of x.
	.br
	acos returns arccosine of x.
	.s
diagnostics
	.s
	If the absolute value of the argument is greater than 1 the
	error message 'asin, acos arg bad', followed by the value of the
	argument, is written to stderr.
	The algorithm returns a result of HUGE. 
	.s
internal
	.s
	Algorithm from Cody and Waite pp. 174-193.
	.s
author
	.s
	Hamish Ross.
	.s
date
	.s
	8-Jan-85
#endif

#include <math.h>

static int flag;
double asin(x)			
double x;
{
    double _asico();

    flag = 0;
    return(_asico(x));
}

double acos(x)
double x;
{
    double _asico();

    flag = 1;
    return(_asico(x));
}

static double p1 = -0.273684945241642560e2;
static double p2 =  0.572082278778917314e2;
static double p3 = -0.396888629975048773e2;
static double p4 =  0.101525222338064636e2;
static double p5 = -0.696745734473506464e0;
static double q0 = -0.164210967144985608e3;
static double q1 =  0.417144302482604126e3;
static double q2 = -0.381863033617501493e3;
static double q3 =  0.150952708410306047e3;
static double q4 = -0.238238591536702388e2;
static double a1 =  PI_BY_4;
static double b[] = {PI_BY_2, PI_BY_4};
static double _asico(x)
double x;
{
    double g, p, q, y;
    double sqrt();
    int i;

    y = x < 0.0 ? -x : x;
    if (y < ROOT_EPS)
	i = flag;
    else {
	if (y > 0.5) {
	    i = 1 - flag;
	    if (y > 1.0) {
		cmemsg(FP_ARSC, &x);
		return(HUGE);
	    }
	    g = ((0.5 - y) + 0.5) / 2.0;
	    y = -2.0 * sqrt(g);
	}
	else {
	    i = flag;
	    g = y * y;
	}
	p = ((((p5*g + p4) * g + p3) * g + p2) * g + p1) * g;
	q = ((((g + q4) * g + q3) * g + q2) * g + q1) * g + q0;
	y = y + y * (p / q);
    }
    if (flag == 0) {
	if (i != 0)
	    y = (a1 + y) + a1;
	if (x < 0.0)
	    return(-y);
	else
	    return(y);
    }
    else {
	if (x < 0.0)
	    return((b[i] + y) + b[i]);
	else {
	    if(i != 0)
		return((a1 - y) + a1);
	    else
		return(-y);
	}
    }
}