/*
	This C language subroutine computes that arctangent
	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 6/8/90
	Updated on 5/94 by AS

	Fixed previous 3.0 bug on 7/94 by AS in calculating Q(g) at
        the end of test_for_eps:.

	#include <math.h>		    <Header>
	double atan(double x);		    <Prototype>
	double atan2(double x, double y);   <Prototype>
*/

#include "lib_glob.h"
#include "mth_glob.h"
#include "atan.h"

.SEGMENT/CODE	Code_Space_Name;
.FILE 		RTL_FILENAME;

.EXTERN		___lib_dtof, ___lib_ftod;

.PRECISION=MACHINE_PRECISION;

.GLOBAL	    	_atan, _atan2;

_atan:		JUMP (PC, save_stack) (DB);
		R0=R0-R0, put(R3);
		put(R7);

_atan2:		R0=R0-R0, put(R3);
		R0=R0+1, put(R7);

save_stack:	put(R10);
		put(R11);
		put(R13);
		put(R14);
		put(R15);
		
		reads(R4,1);
		reads(R8,2);
		reads(R12,3);

		CALL (PC, ___lib_dtof) (DB);
		R10=PASS R0, R0=R4;
		R1=PASS R8, dm_ptr=R10;	/* dm_ptr holds flag for fcn */

		R14=R14-R14, F15=F0;
		R10=PASS R10, R13=R15;	/*Test for atan2 routine*/
		IF EQ JUMP (PC, get_f);	/*Go straight to atan routine*/

do_divide:	CALL (PC, ___lib_dtof) (DB);
	    	R0=PASS R12;
		reads(R1,4);

		R14=PASS R0;		/*Test for 0 divisor*/
		IF NE JUMP (PC, overflow_tst);/*Test for overflow division*/
		R15=PASS R15;		/*Test for dividend also 0*/
		IF NE JUMP (PC, overflow);	/*Same path as overflow*/

error:		JUMP (PC, restore_state) (DB);
		F0=0;			/*Set return value???*/
		FETCH_RETURN

overflow_tst:	R10=LOGB F15;		/*Get exponent of X*/
		R1=LOGB F14;		/*Get exponent of Y*/
		R10=R10-R1;		/*Compute result exponent*/
		R1=124;			/*Max exponent - 3*/
		COMP(R10,R1);
		IF GE JUMP (PC, overflow);	/*Quotient will overflow*/
		R1=-R1;
		COMP(R10,R1);
		IF LE JUMP (PC, underflow);	/*Quotient will underflow*/

do_division:	CALL (PC, ___float_divide) (DB);
		F11=2.0;		/*Needed for divide*/
		F12=PASS F14, F7=F15;	/*Set N=X, D=Y, compute divide*/

		F13=PASS F7;

get_f:		F11=2.0;		/* Used in divide	*/
		F13=ABS F13;		/*Get absolute value of X or X/Y*/
		F12=1.0;		
		COMP(F13,F12);
		F10=0;			/*Set N to zero*/
		IF LE JUMP (PC, tst_f);	

		CALL (PC, ___float_divide) (DB);
		F12=PASS F13, F7=F12;	/*Set N=1.0, D=f*/
		NOP;			/*Dead cycle*/

		F13=PASS F7, F10=F11;	/*F7=N*R0*R1*R2*R3 Set N=2.0*/

tst_f:		R10=FIX F10;		/*Convert N to integer*/
		F12=two_minus_sqrt_3;
		COMP(F13,F12);
		IF LT JUMP (PC, tst_for_eps);	
		R10=R10+1;		/*N=N+1*/
		F1=0.5;
		F8=sqrt_3_minus_1;
		F7=F8*F13;		/*F7=A*f*/
		F7=F7-F1;		/*F7=(A*f-0.5)*/
		F7=F7-F1;		/*F7=((A*f-0.5)-0.5)*/
		F7=F7+F13;		/*F7=((A*F-0.5)-0.5)+f)*/
		CALL (PC, ___float_divide) (DB);
		F8=sqrt_3;
		F12=F8+F13;		/*F12=sqrt(3)+f		*/
		
		F13=PASS F7;

tst_for_eps:	F1=ABS F13;		/*Test for small input	*/
		F12=eps;
		COMP(F1,F12);
		IF LE JUMP (PC, tst_N);
		F8=F13*F13;		/*g = f*f		*/
		F12=p1;
		F7=F8*F12;		/*P(g)=p1*g		*/
		F12=p0;
		F7=F7+F12;		/*P(g) = p1*g + p0 	*/
		F7=F7*F8;		/*g*P(g) = (p1*g+p0)*g	*/
		F12=q1;
		F12=F8+F12;		/*Q(g) = (g + q1)	*/
		F12=F12*F8;		/*Q(g) = (g + q1) * g	*/
		CALL (PC, ___float_divide) (DB);
		F2=q0;			
		F12=F12+F2;		/*Q(g) = (g + q1)*g+q0	*/

		F7=F7*F13;		/* f*R 			*/
		F13=F7+F13;		/* f+f*R 		*/

tst_N:		R8=1;
		COMP(R10,R8);
		IF GT F13=-F13;		/*If N>1 negate result	*/
		R10=PASS R10;
		IF EQ JUMP (PC, adjust_by_AN) (DB);
		F12=0;
		R10=R10-1;
		IF EQ JUMP (PC, adjust_by_AN) (DB);
		F12=pi_over_6;
		R10=R10-1;
		IF EQ JUMP (PC, adjust_by_AN) (DB);
		F12=pi_over_2;
		R10=R10-1;

		F12=pi_over_3;
adjust_by_AN:	F12=F13+F12, R0=dm_ptr;
		R0=PASS R0;
		IF EQ JUMP (PC, tst_sign_x);
tst_sign_y:	F13=pi;
		F14=PASS F14;
		IF LT F12=F13-F12;	/* Result = pi - result*/
tst_sign_x:	F15=PASS F15;
		IF LT F12=-F12;

		F0=PASS F12, FETCH_RETURN
				
restore_state:	CALL (PC, ___lib_ftod) (DB);
		get(R15,1);
		get(R14,2);

		get(R13,3);
		get(R11,4);
		get(R10,5);
		get(R7,6);
		get(R3,7);
		RETURN (DB);
		RESTORE_STACK
		RESTORE_FRAME

overflow:	JUMP (PC, tst_sign_x) (DB);
		NOP;
		F12=pi_over_2;
		
underflow:	JUMP (PC, tst_sign_y) (DB);
		NOP;
		F12=0;

.ENDSEG;