/*
	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

	#include <math.h>		<Header>
	float atanf(float x);		<Prototype>
	float atan2f(float x, float y);	<Prototype>
*/

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

.SEGMENT/CODE	Code_Space_Name;
.FILE 		RTL_FILENAME;

.GLOBAL	    	_atanf, _atan2f;

_atanf:		JUMP (PC, save_stack) (DB);
		R0=R0-R0, put(R7);
		put(R11);

_atan2f:	R0=R0-R0, put(R7);
		R0=R0+1, put(R11);

save_stack:	put(R5);
		R5=PASS R0, put(R3);	/*Test for atan2 routine*/ 

		IF EQ JUMP (PC, get_f) (DB); /*Go straight to atan routine*/
		R0=PASS R4;
		R3=PASS R8;		/*Test for 0 divisor*/

do_divide:	IF NE JUMP (PC, overflow_tst); /*Test for overflow division*/
		R4=PASS R4;		/*Test for dividend also 0*/
		IF NE JUMP (PC, overflow); /*Same path as overflow*/

input_error:	JUMP (PC, restore_state) (DB);
		F0=0.0;
		FETCH_RETURN

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

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

		F0=PASS F7;

get_f:		F11=2.0;		/*For divides*/
		F0=ABS F0, R1=pm_0;	/*Get abs val of X or X/Y, N=0*/
		F12=1.0;		
		COMP(F0,F12);
		IF LE JUMP (PC, tst_f);	

		CALL (PC, ___float_divide) (DB);
		F12=PASS F0, F7=F12;	/*Set N=1.0, D=f*/
		R1=2;			/*Set N to two*/

		F0=PASS F7;		/*F7=N*R0*R1*R2*R3*/

tst_f:		F12=two_minus_sqrt_3;
		COMP(F0,F12);
		IF LE JUMP (PC, tst_for_eps);	
		R1=R1+1;		/*N=N+1*/
		F8=sqrt_3_minus_1;
		F7=F8*F0;		/*F7=A*f*/
		F8=0.5;
		F7=F7-F8;		/*F7=(A*f-0.5)*/
		F7=F7-F8;		/*F7=((A*f-0.5)-0.5)*/
		F7=F7+F0;		/*F7=((A*F-0.5)-0.5)+f)*/
		CALL (PC, ___float_divide) (DB);
		F8=sqrt_3;
		F12=F8+F0;		/*F12=sqrt(3)+f*/
		
		F0=PASS F7;

tst_for_eps:	F8=ABS F0;		/*Test for small input*/
		F12=eps;
		COMP(F8,F12);
		IF LE JUMP (PC, tst_N);
		F8=F0*F0;		/*g = f*f*/
		F2=p1;
		F7=F8*F2;		/*P(g)=p1*g*/
		F2=p0;
		F7=F7+F2;		/*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, F8=F0;	/*Q(g) = (g + q1) * g, F8=f*/
		CALL (PC, ___float_divide) (DB); /* Trashes F0 ^^^^*/
		F2=q0;			
		F12=F12+F2;		/*Q(g) = (g + q1) * g + q0 */

		F7=F7*F8;		/* f*R 			*/
		F0=F7+F8;		/* f+f*R		*/

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

		F12=pi_over_3;
adjust_by_AN:	F12=F0+F12;
		R5=PASS R5;
		IF EQ JUMP (PC, tst_sign_x);
tst_sign_y:	F0=pi;
		F3=PASS F3;		/* was r8 (2nd parameter) */
		IF LT F12=F0-F12;	/* Result = pi - result*/
tst_sign_x:	F4=PASS F4;
		IF LT F12=-F12;

		F0=PASS F12, FETCH_RETURN
				
restore_state:	get(R3,1);
		get(R5,2);
		get(R11,3);
		get(R7,4);
		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;