/* arcsin.c */

/* For |x| < 0.65 uses Abramowitz & Stegun continued fraction formula
* 4.4.44  and matrix recurrence A&S 3.10.1 for evaluation starting 
* with second term to cope with a1 = x whereas higher terms involve x*x.
* For |x| >= 0.65 use A&S series 4.4.41  
*  S E Peterson  3/92    */

#define __LIBRARY__

#include <math.h>
#include <stdlib.h>
#include <errno.h>
#include <string.h>
#include <stdio.h>
#include <float.h>

static double confrac _LIB_F1_(double,x)
{
	char msg[70], msg2[40]; 
	double     x2, res, old, diff;
	double	   an, An, Bn;
	double	   a1,a2,b1,b2;
	int n, bn;
	(void)strcpy(msg, "No convergence - increase no. of iterations or value of epsilon");
	(void)strcpy(msg2, "An or Bn exceeds numeric limits");
	a1=x;
	a2=0.0;
	b1=1.0;
	b2=1.0;

	n=2;  
	old = HUGE_VAL;			 		
	x2 = x*x;
	while (n < 50){  		/* control for no. of iterations */
		if (n%2 == 0)	    /* if n is even increment an,else stays same */
			an = (double)n *(double)(1-n) * x2;
		bn = 2*n - 1;
		An = (a1* (double)bn)+(a2*an);
		Bn = (b1* (double)bn)+(b2*an);
		if (fabs(An) >= HUGE_VAL || fabs(Bn) >= HUGE_VAL){
			return (HUGE_VAL);
		}	
		res = An/Bn;
		diff = old - res;
		if ( fabs(diff) < DBL_EPSILON ){
			return (res);
		}
		else{
			old=res;
			a2=a1;			/* update matrix */
			a1=An;
			b2=b1;
			b1=Bn;
			n++;	
		}
	}
	(void)printf("n = %6d\n",n); 
	(void)printf("%s\n", msg); 
	return (DBL_EPSILON);
}

static double series _LIB_F1_(double,z)
{
	static double coeff[19] = {1.0,
					8.33333333333333333333E-2,
					1.87500000000000000000E-2,
					5.58035714285714285714E-3,
					1.89887152777777777778E-3,
					6.99129971590909090909E-4,
					2.71136944110576923077E-4,
					1.09100341796875000000E-4,
					4.51242222505457261029E-5,
					1.90656436117071854441E-5,
					8.19368731407892136347E-6,
					3.57056927421818608823E-6,
					1.57402595505118370056E-6,
					7.00688192241445735649E-7,
					3.14533061665033215079E-7,
					1.42216292935641362302E-7, 
					6.47111067761133282064E-8,
					2.96094097811711825281E-8,
					1.36154380562817937676E-8
					};

	double ans;
	int ncf, i;
	
	ncf = 9;		/*  9 coeffs exhausts acuracy of MFFP - 
				*   increase this when true double available */
	
	ans = coeff[ncf];
	if (ans < 1.0E-15)
		ans = 0.0;
	for ( i = (ncf-1); i >= 0;i--){
		ans = coeff[i] + (z * ans);
		if (ans < 1.0E-15)	
			ans = 0.0;
	}
		return (ans);
}
													
double asin(x)
double x;

{
	double ans, z;
	if (x < -1.0 || x > 1.0){
		errno = EDOM;
		return (-HUGE_VAL);
	}	
	if (fabs(x) < 0.65){
		ans = confrac(x);
		if (ans == HUGE_VAL)
			return (HUGE_VAL);
		ans *= sqrt(1.0 - x*x);
		return(ans);
	}	
	else{
		z = 1.0  - fabs(x); 
		ans = series(z);
		ans = PID2 - (sqrt(2.0 * z)) * ans;
		if (x < 0.0)
			ans = -ans;
		return (ans);
	}	
}