{
    Subroutine to compute the Sine or Cosine values of
    a floating point input.

	    Y=SIN(X) or
	    Y=COS(X)

	Calling Registers
    F0 = Input Value X=[6E-20, 6E20] 
    l_reg=0

	Result Registers
    F0 = Sine (or Cosine) of input Y=[-1,1]

	Altered Registers
    F0, F2, F4, F7, F8, F12
    i_reg

	Computation Time
    38 Cycles

    Version 0.03    7/4/90	Gordon A. Sterling
}

#include "asm_glob.h"

.SEGMENT/PM	Assembly_Library_Code_Space;

.PRECISION=MACHINE_PRECISION;

#define half_PI    1.57079632679489661923

.GLOBAL	cosine, sine;

cosine:		i_reg=sine_data;	    {Load pointer to data}
		F8=ABS F0;		    {Use absolute value of input}
		F12=0.5;		    {Used later after modulo}
		F2=1.57079632679489661923;  { and add PI/2}
		JUMP compute_modulo (DB);   {Follow sin code from here!}
		F4=F8+F2, F2=mem(i_reg,1);   
		F7=1.0;			    {Sign flag is set to 1}

sine:		i_reg=sine_data;	    {Load pointer to data}
		F7=1.0;			    {Assume a positive sign}
		F12=0.0;		    {Used later after modulo}
		F8=ABS F0, F2=mem(i_reg,1);	
		F0=PASS F0, 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, F4=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 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}
		RTS (DB), F4=F4*F8;	    {Compute f*R}
		F12=F4+F8;		    {Compute Result=f+f*R}
compute_sign:	F0=F12*F7;		    {Restore sign of result}
		RTS;			    {This return only for sin(eps)=eps path}

.ENDSEG;

.SEGMENT/SPACE Assembly_Library_Data_Space;

.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;