/*
	This C language subroutine computes the value of 
	X raised to the Y

	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 pow(double x, double y);	    <Prototype>
*/

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

.SEGMENT/CODE	Code_Space_Name;
.FILE RTL_FILENAME;

.GLOBAL	    _pow;

_pow:		put(R3);
		put(R5);
		put(R6);
		put(R7);
		put(R9);
		put(R10);
		put(R11);
		put(R13);
		put(R14);
		put(R15);
		
		reads(R4,1);
		reads(R8,2);
		reads(R12,3);

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

		CALL (PC, ___lib_dtof) (DB);
		R0=PASS R12, F5=F0;
		reads(R1,4);

		F0=PASS F5, F1=F0;		/*Test for x<=0*/
		m_reg=dm_1;

		IF LT JUMP (PC, x_neg_error);	/*Report error*/
		IF EQ JUMP (PC, check_y);	/*Test y input*/
determine_m:	R2=LOGB F0;			/*Get exponent of x input*/
		R2=R2+1;			/*Reduce to proper format*/
		R3=-R2;				/*Used to produce f*/
determine_g:	F15=SCALB F0 BY R3;		/* .5 <= g < 1*/
determine_p:	i_reg=a1_values;
		R14=R14-R14, R10=i_reg;
		R13=9;
		R14=R14+1, F12=mem(9,i_reg);	/*Get A1(9)*/
		COMP(F12,F15), F12=mem(5,i_reg);/*A1(9) - g*/
		F11=mem(13,i_reg);
		IF GE F12=F11;		        /*Use A(13) next*/
		IF GE R14=R13;		        /*IF (g<=A1(9)) p=9*/
		R9=R10+R14;				
		i_reg=R9;
		R13=4;
		COMP(F12,F15), F12=mem(2,i_reg);/*A1(p+4) - g*/
		F11=mem(6,i_reg);
		IF GE F12=F11;
		IF GE R14=R14+R13;	        /*IF (g<=A1(p+4)) p=p+4*/
		R13=2;
		COMP(F12,F15);		        /*A1(p+2) - g*/
		IF GE R14=R14+R13;	        /*IF (g<=A1(p+4)) p=p+2*/
determine_z:	R14=R14+1, R4=R14;
		m_reg=R14;
		i_reg=a1_values;
		R11=ASHIFT R14 BY -1;		/*Comp (p+1)/2	*/
		F12=mem(m_reg,i_reg);		/*Fetch A1(p+1) */
		i_reg=a2_values;	        /*Correction array*/
		F0=F12+F15, m_reg=R11;	        /*g + A1(p+1)*/
		F14=F15-F12, F11=mem(m_reg,i_reg);
		CALL (PC, ___float_divide) (DB);
		F7=F14-F11, F12=F0;		/*[g-A1(p+1)]-A2((p+1)/2)*/
		F11=2.0;

		F7=F7+F7;			/* z = z + z*/
determine_R:	i_reg=power_array;
		F8=F7*F7, F9=mem(p2,i_reg);	/* v = z * z */
		F10=F8*F9, F9=mem(p1,i_reg);	/* p2*v */
		F10=F10+F9;			/* p2*v + p1 */
		F10=F10*F8;			/* (p2*v + p1) * v */
		F10=F10*F7, F9=mem(K,i_reg);	/* (p2*v + p1) * v * z */
determine_u2:	F11=F10*F9;			/* K*R */
		F11=F10+F11;			/* K + K*R */
		F9=F9*F7;			/* z*K */
		F9=F9+F11;			/* R + z*K */
		F9=F9+F7;			/* (R + z*K) + z*/
determine_u1:	R3=16;
		R2=R2*R3 (SSI);			/* m*16 */
		R2=R2-R4;			/* m*16-p */
		R3=-4;
		F2=FLOAT R2 BY R3;		/*FLOAT(m*16-p)*.0625*/
determine_w:	R4=4;				/*Used in reduce*/
		R8=FIX F1 BY R4;
		F8=FLOAT R8 BY R3;		/* y1=REDUCE(Y) */
		F7=F1-F8;			/* y2 = y - y1 */
		F15=F9*F1;			/* U2*Y */
		F14=F2*F7;			/* U1*Y2 */
		F15=F14+F15;			/* W = U2*Y + U1*Y2*/
		R14=FIX F15 BY R4;
		F14=FLOAT R14 BY R3;		/* W1=REDUCE(W) */
		F13=F15-F14;			/* W2 = W - W1 */
		F12=F2*F8;			/* U1*Y1 */
		F12=F14+F12;			/* W = W1 + U1*Y1 */
		R14=FIX F12 BY R4;		
		F14=FLOAT R14 BY R3;		/* W1 = REDUCE(W) */
		F11=F12-F14;			/* W-W1 */
		F13=F11+F13;			/* W2 = W2 + (W-W1) */
		R12=FIX F13 BY R4;
		F12=FLOAT R12 BY R3;		/* W = REDUCE(W2) */
		F10=F12+F14;
		R10=FIX F10 BY R4;		/* IW1 = INT(16*(W1+W)) */
		F13=F13-F12, R9=mem(bigx,i_reg);/* W2 = W2 - W */
		COMP(R10,R9), R9=mem(smallx,i_reg);/* Test for overflow */
		IF GE JUMP (PC, overflow);
		COMP(R10,R9);			/* Test for underflow */
		IF LE JUMP (PC, underflow);
flow_to_a:	F13=PASS F13;			/* W2 must be <=0 */
		IF LE JUMP (PC, determine_mp);
		F8=0.0625;
		F13=F13-F8;
		R10=R10+1;
determine_mp:	R10=PASS R10, R8=dm_1;
		IF LT R8=R8-R8;			/* I=0 if IWI < 0 */
		R6=ABS R10;
		R7=ASHIFT R6 BY -4;
		R6=PASS R10;
		IF LT R7=-R7;
		R7=R7+R8;			/* m' = IW1/16 + I */
		R6=ASHIFT R7 BY 4;		/* m'*16 */
		R6=R6-R10, F5=mem(q5,i_reg);	/* p'  = 16*m' - IW1 */
determine_Z:	F4=F5*F13, F5=mem(q4,i_reg);	/* q5*W2 */
		F4=F4+F5, F5=mem(q3,i_reg);	/* q5*W2 + q4 */
		F4=F4*F13;			/* (q5*W2+q4)*W2 */
		F4=F4+F5, F5=mem(q2,i_reg);	/* (q5*W2+q4)*W2+q3 */
		F4=F4*F13;			/* ((q5*W2+q4)*W2+q3)*W2 */
		F4=F4+F5, F5=mem(q1,i_reg);	/* ((q5*W2+q4)*W2+q3)*W2+q2 */
		F4=F4*F13;			/* (((q5*W2+q4)*W2+q3)*W2+q2)*W2 */
		R6=R6+1;			/* Compute p'+1 */
		F4=F4+F5, m_reg=R6;		/* Q(W2)=(((q5*W2+q4)*W2+q3)*W2+q2)*W2+q1*/
		i_reg=a1_values;
		F4=F4*F13, F5=mem(m_reg,i_reg);	/* Z = W2*Q(W2), Fetch A1(p'+1) */
		F4=F4*F5;			/* A1(p'+1)*Z */
		F4=F4+F5;			/* Z = A1(p'+1) + A1(p'+1)*Z */
		F0=SCALB F4 BY R7, FETCH_RETURN	/* Result = ADX(Z,m') */
						
restore_state:	CALL (PC, ___lib_ftod) (DB);
		get(R15,1);
		get(R14,2);
		get(R13,3);
		get(R11,4);
		get(R10,5);
		get(R9,6);
		get(R7,7);
		get(R6,8);
		get(R5,9);
		get(R3,10);
		RETURN (DB);
		RESTORE_STACK
		RESTORE_FRAME

underflow:	ram_ireg=_errno;
		JUMP (PC, restore_state) (DB);
		rammem(ram_ireg, ram_0)=ERANGE;
		R0=R0-R0, FETCH_RETURN

x_neg_error:	ram_ireg=_errno;
		JUMP (PC, restore_state) (DB);
		rammem(ram_ireg, ram_0)=EDOM;
		R0=R0-R0, FETCH_RETURN

check_y:	F1=PASS F1, FETCH_RETURN		/*If y>0 return x*/
		IF GT JUMP (PC, restore_state);

overflow:	ram_ireg=_errno;
		JUMP (PC, restore_state) (DB);
		rammem(ram_ireg, ram_0)=ERANGE;
		R0=R0-R0, FETCH_RETURN

.ENDSEG;

.SEGMENT/SPACE	Data_Space_Name;

.PRECISION=MEMORY_PRECISION;

.VAR	a1_values[18] = 0,
			0x3F800000,
			0x3F75257D,
			0x3F6AC0C6,
			0x3F60CCDE,
			0x3F5744FC,
			0x3F4E248C,
			0x3F45672A ,
			0x3F3D08A3 ,
			0x3F3504F3 ,
			0x3F2D583E ,
			0x3F25FED6 ,
			0x3F1EF532 ,
			0x3F1837F0 ,
			0x3F11C3D3 ,
			0x3F0B95C1 ,
			0x3F05AAC3 ,
			0x3F000000;

.VAR	a2_values[9] =	0,
			0x31A92436,
			0x336C2A95,
			0x31A8FC24,
			0x331F580C,
			0x336A42A1,
			0x32C12342,
			0x32E75624,
			0x32CF9891;
			
.VAR	power_array[10]= 0.833333286245E-1,	/* p1 */
			0.125064850052E-1,	/* p2 */
			0.693147180556341,	/* q1 */
			0.240226506144710,	/* q2 */
			0.555040488130765E-1,	/* q3 */
			0.961620659583789E-2,	/* q4 */
			0.130525515942810E-2,	/* q5 */
			0.44269504088896340736,	/* K */
			2032,			/* bigx */
		       -2015;			/* smallx */
.ENDSEG;