/*
	This C language subroutine computes the Sine 
	or Cosine values of a 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>
	float sinf(float x);	    <Prototype>

	50 Cycles
*/

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

.PRECISION=MACHINE_PRECISION;

.SEGMENT/CODE	Code_Space_Name;
.FILE 		RTL_FILENAME;

.GLOBAL	    	_sinf, _cosf;

_cosf:		JUMP (PC, reg_save) (DB);
		R0=R0-R0, put(R3);  	  	/*Clear flag*/
		i_reg=sine_data;

_sinf:		F0=1.0;		    		/*Set flag*/
		R12=R12-R12, put(R3);		/*Used later after modulo*/

reg_save:	put(R5);
		R1=PASS R4, put(R7);

		F0=PASS F0, R5=MODE1;		/*Read input value*/ 
		F7=1.0;			    	/*Sign flag is set to 1*/

		IF NE JUMP (PC, ___sin_path) (DB); /*Sine requested*/
		BIT CLR MODE1 65536;	     	/*Set to 40-bit mode*/
		F8=ABS F1;		    	/*Use absolute value of input*/
		
		F2=1.57079632679489661923;  /* and add PI/2*/
		JUMP (PC, compute_modulo) (DB); /*Follow sin code from here!*/
		F4=F8+F2, F2=mem(i_reg,1);
		F12=0.5;		    	/*Used later after modulo for cos*/

___sin_path:	i_reg=sine_data;
		F8=ABS F1, F2=mem(i_reg,1);	
		F1=PASS F1, F4=F8;
		IF LT F7=-F7;		    /*If input was negative, invert sign*/

compute_modulo:	F4=F4*F2;	    	    /*Compute fp modulo value*/
		R2=FIX F4;		    /*Round nearest fractional portion*/
		BTST R2 BY 0;		    /*Test for odd number*/
		IF NOT SZ F7=-F7;	    /*Invert sign if odd modulo*/
		F4=FLOAT R2;		    /*Return to fp*/
		F4=F4-F12, F2=mem(i_reg,1); /*Add cos adjust if necessary, F3=XN*/

compute_f:	F12=F2*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;		    /*Need magnitude for test*/
		F4=F12-F4, F12=F8;	    /*Check for sin(x)=x*/		    
		IF LT JUMP (PC, compute_sign); /*Return with result in F1*/

compute_R:	F12=F12*F12, F4=mem(i_reg,1);
		LCNTR=6, DO compute_poly UNTIL LCE;
		    F4=F12*F4, F2=mem(i_reg,1);	/*Compute sum*g*/
compute_poly:	    F4=F2+F4;			/*Compute sum=sum+next r*/
		F4=F12*F4;			/*Final multiply by g*/
	    	F4=F4*F8;		    /*Compute f*R*/
		F12=F4+F8;		    /*Compute Result=f+f*R*/
compute_sign:	F0=F12*F7, FETCH_RETURN	    /*Restore sign of result*/

restore_state:	F0=RND F0, MODE1=R5;
		get(R7,1);
		get(R5,2);
		get(R3,3);
		RETURN (DB);
		RESTORE_STACK
		RESTORE_FRAME
.ENDSEG;

.SEGMENT/SPACE	Data_Space_Name;

.PRECISION=MEMORY_PRECISION;

.VAR sine_data[11] =	 0.31830988618379067154,	/*1/PI*/ 
			 3.14160156250000000000,	/*C1, almost PI*/
			-8.908910206761537356617E-6,	/*C2, PI=C1+C2*/
			 9.536743164E-7,		/*eps, sin(eps)=eps*/
			-0.737066277507114174E-12,	/*R7*/
			 0.160478446323816900E-9,	/*R6*/
			-0.250518708834705760E-7,	/*R5*/
			 0.275573164212926457E-5,	/*R4*/
			-0.198412698232225068E-3,	/*R3*/
			 0.833333333327592139E-2,	/*R2*/
			-0.166666666666659653;		/*R1*/

.ENDSEG;