/*
	This C language subroutine computes the ArcSine
	and ArcCosine of its floating point input.

	The Run Time Library for the C Language.

	Gordon A. Sterling	(617) 461 - 3076
	DSP Development Tools Engineering

	Created on 5/7/90
	Updated on 5/94 by AS

	#include <math.h>		    <Header>
	float acosf(float x);		    <Prototype>
	float asinf(float x);		    <Prototype>
*/

#include "lib_glob.h"
#include "mth_glob.h"
#include "asin.h"

.SEGMENT/CODE	Code_Space_Name;
.FILE 		RTL_FILENAME;

.GLOBAL	    	_asinf, _acosf;

.EXTERN	    	_sqrtf;

_acosf:		JUMP (PC, arc_core) (DB);
		R0=R0-R0, put(R6);
		R0=R0+1, put(R7);	/*Set flag for acos*/

_asinf:		put(R6);
		R0=R0-R0, put(R7);	/*Clear flag for asin*/

arc_core:	R2=R2-R2, put(R10);
		R2=R2+1, put(R11);
		R1=PASS R4, put(R13);	/*Save flag	*/
		F14=ABS F4, put(R14);	/*Y=ABS(X)	*/
		
		F8=0.5;
		COMP(F8,F14), R11=R0;	/*In case we jump, we need flag in R11 */
		IF GT JUMP (PC, small_input)(DB);
		R13=PASS R0;
		i_reg=arc_data;		/*Point to data array*/

		R2=R2-R0, R6=R1;	/*i=1-Flag, save r1 temp in r6*/

		F13=1.0;
		COMP(F13,F14), R11=R2;	/* Hold i in R11*/
		IF LT JUMP (PC, domain_error);
		F8=F13-F14, R10=dm_M1;
		F14=SCALB F8 BY R10, R13=R0; /*Divide by 2, Hold R0*/
		R10=-R10, CALLER_HOLD(R0)
		R4=PASS R14, CALLER_SWAP
		RTLCALL (PC, _sqrtf) (DB);
		SAVE_OLD_FRAME(R0)
label_1:	SAVE_RET_ADDR

		F14=SCALB F0 BY R10, F8=F14; 	/*Multiply by 2*/
		F14=-F14, R1=R6;		/*Place x input into R1 again, sqrtf */
						/*might have trashed it	*/
		
compute_R:	R2=PASS R11, F6=mem(p3,i_reg); /*Place i back into R2 */
		F7=F6*F8, F6=mem(p2,i_reg); /*Compute p3*g	*/
		F7=F7+F6;		    /*Compute p3*g+p2	*/
		F7=F7*F8, F6=mem(p1,i_reg); /*Compute (p3*g+p2)*g*/
		F7=F7+F6;		    /* (p3*g+p2)*g+p1 	*/
		F7=F7*F8;		    /* g*P(g)=[(p3*g+p2)*g+p1]*g*/
		F12=PASS F8, F6=mem(q2,i_reg);
		F12=F12+F6;		    /* Compute g+q2 	*/
		F12=F12*F8, F6=mem(q1,i_reg); /* (g+q2)*g 	*/
		F12=F12+F6;		    /* [(g+q2)*g+q1]	*/
		CALL (PC,___float_divide) (DB);
		F12=F12*F8, F6=mem(q0,i_reg);
		F12=F12+F6, F11=mem(two,i_reg);/* Q(g)=[(g+q2)*g+q1]+q0*/

		F7=F7*F14;		    /* Y*R(g) 		*/
		F7=F7+F14;	    	    /* Y+Y*R(g) 	*/

test_flag:	m_reg=R2;		    /*Load i in MOD reg	*/
		MODIFY(i_reg, a0);	    /*Point to a0	*/
		R13=PASS R13;
		IF NE JUMP (PC, test_x_input);

    		F14=mem(m_reg,i_reg);
		F7=F7+F14;		    /*Compute ai+result	*/
		F0=F7+F14;		    /*Compute (ai+result)+ai */
		JUMP (PC, restore_state) (DB);
		F1=PASS F1;		    /*Test x input*/
		IF LT F0=-F0;

test_x_input:	F1=PASS F1;		    /*Test x input*/
		IF GE JUMP (PC, compute_a);
		MODIFY(i_reg,db0-a0);	    /* Point to db0 */
		F14=mem(m_reg,i_reg);
		JUMP (PC, restore_state) (DB);
		F7=F7+F14;
		F0=F7+F14;

compute_a:	F14=mem(m_reg,i_reg);
		F7=F14-F7;
		F0=F7+F14;
		
restore_state:	FETCH_RETURN
		get(R14,1);
		get(R13,2);
		get(R11,3);
		get(R10,4);
		get(R7,5);
		get(R6,6);
		RETURN (DB);
		RESTORE_STACK
		RESTORE_FRAME

small_input:	R2=PASS R0, F4=mem(eps,i_reg);
		COMP(F4,F14), F7=F14;
		IF GE JUMP (PC, test_flag), R13=PASS R0; /*setup R13 in test_flag */
		JUMP (PC, compute_R) (DB);
		F8=F14*F14;		/*g = Y * Y */
		nop;

domain_error:	JUMP (PC, restore_state) (DB), R0=R0-R0;
		ram_ireg=_errno;
		rammem(ram_ireg,ram_0)=EDOM;

.ENDSEG;

.SEGMENT/SPACE	Data_Space_Name;

.PRECISION=MEMORY_PRECISION;

.VAR	arc_data[13] =	 0.0,			    /*eps*/
			-0.27516555290596E+1,	    /*p1*/
			 0.29058762374859E+1,	    /*p2*/
			-0.59450144193246,	    /*p3*/
			-0.16509933202424E+2,	    /*q0*/
			 0.24864728969164E+2,	    /*q1*/
			-0.10333867072113E+2,	    /*q2*/
			 0.10000000000000E+1,	    /*q3*/
			 0.0,			    /*a0*/
			 0.78539816339744830962,    /*a1*/
			 1.57079632679489661923,    /*db0*/
			 0.78539816339744830962,    /*b1*/
			 2.0;			    /*Two*/

.ENDSEG;