/*
	This C language subroutine computes the tangent
	of the floating point input.
	
	The Run Time Library for the C Language.

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

	Created on 4/28/90
	Updated on 5/94 by AS

	#include <math.h>	    <Header>
	double tan(double x);	    <Prototype>
	double cot(double x);	    <Prototype>
*/

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

#define eps1 9.0E-7;		    /*Smallest value*/

.SEGMENT/CODE	Code_Space_Name;
.FILE 		RTL_FILENAME;

.GLOBAL	    	_tan, _cot;

_tan:		JUMP (PC, tan_core) (DB);
		R0=R0-R0, put(R3);
		i_reg=tangent_data;

_cot:		R0=1;
		put(R3);
		i_reg=tangent_data;

tan_core:	put(R5);
		put(R7);
		put(R11);

		reads(R4,1);
		reads(R8,2);

		CALL (PC, ___lib_dtof) (DB);
		R5=PASS R0, R0=R4;
		R1=PASS R8;

		R5=PASS R5;
		IF EQ JUMP (PC, compute_modulo) (DB);
		F8=ABS F0, F2=mem(i_reg,1); /*Use absolute value of input*/
		F4=eps1;
	
		F8=MAX(F8,F4);

compute_modulo:	F4=F8*F2, F1=F0;	    /*Compute fp modulo value*/
		R2=FIX F4, F12=mem(i_reg,1);/*Rnd nearest fractional portion*/
		F4=FLOAT R2, R0=R2;	    /*Return to fp*/

compute_f:	F12=F12*F4, F2=mem(i_reg,1);/*Compute XN*C1*/
		F2=F2*F4, F12=F8-F12;	    /*Compute |X|-XN*C1, and XN*C2*/
		F8=F12-F2, F4=mem(i_reg,1); /*Compute f=(|X|-XN*C1)-XN*C2*/
		F12=ABS F8, F7=F8;
		F4=F12-F4, F12=mem(i_reg,1);/*Check for TAN(x)=x*/		    
		IF LT JUMP (PC, compute_quot);    /*Compute quotient with NUM=f DEN=1*/

compute_P:	F12=F8*F8, F4=mem(i_reg,1); /*g=f*f*/
		F4=F12*F4, F2=mem(i_reg,1); /*Compute p3*g*/
		F4=F2+F4;		    /*Compute (p3*g + p2)*/
		F4=F12*F4, F2=mem(i_reg,1); /*Compute (p3*g + p2)*g*/
		F4=F2+F4;		    /*Compute (p3*g + p2)*g + p1*/
		F4=F12*F4;		    /*Compute ((p3*g + p2)*g + p1)*g*/
		F4=F4*F8;		    /*Compute ((p3*g + p2)*g + p1)*g*f*/
		F8=F4+F8, F4=mem(i_reg,1);  /*Compute f*P(g)*/

compute_Q:	F4=F12*F4, F2=mem(i_reg,1); /*Compute sum*g*/
		F4=F2+F4;		    /*Compute sum=sum+next q*/
		F4=F12*F4, F2=mem(i_reg,1); /*Compute sum*g*/
		F4=F2+F4;		    /*Compute sum=sum+next q*/
		F4=F12*F4, F2=mem(i_reg,1); /*Compute sum*g*/
		F12=F2+F4, F7=F8;	    /*Compute sum=sum+next q*/

compute_quot:	R5=PASS R5, F11=mem(i_reg,1);	/*Test value of IFLAG*/
		IF NE JUMP (PC, cot_path);	/*Cotan route*/

tan_path:	CALL (PC, ___float_divide) (DB);
		BTST R0 BY 0; 
		IF NOT SZ F12=-F7, F7=F12;
		JUMP (PC, restore_state) (DB);
		F0=PASS F7;
		FETCH_RETURN

cot_path:	BTST R0 BY 0;			/*Test for even N*/
		CALL (PC, ___float_divide) (DB);
		IF SZ F12=PASS F7, F7=F12;
		IF NOT SZ F7=-F7;

		F0=PASS F7, FETCH_RETURN

restore_state:	CALL (PC, ___lib_ftod) (DB);
		get(R11,1);
		get(R7,2);
		get(R5,3);
		get(R3,4);
		RETURN (DB);
		RESTORE_STACK
		RESTORE_FRAME
.ENDSEG;

.SEGMENT/SPACE	Data_Space_Name;

.PRECISION=MEMORY_PRECISION;

.VAR tangent_data[13] =	 0.6366197723675834308,		/*2/PI*/
			 1.57080078125,			/*C1, almost PI/2*/
			-4.454455103380768678308E-6,	/*C2, PI/2=C1+C2*/
			 9.536743164E-7,		/*eps, TAN(eps)=eps*/
			 1.0,				/*Used in one path*/
			-0.7483634966612065149E-5,	/*P3*/
			 0.2805918241169988906E-2,	/*P2*/
			-0.1282834704095743847,		/*P1*/
			-0.2084480442203870948E-3,	/*Q3*/
			 0.2334485282206872802E-1,	/*Q2*/
			-0.4616168037429048840,		/*Q1*/
			 1.0,				/*Q0*/
			 2.0;				/*Used in divide*/

.ENDSEG;