/*
	This C language subroutine computes the Root Mean
	Square (RMS) of its input.

	The Run Time Library for the C Language.

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

	Created on 04/19/90 
	Updated on 5/94 by AS

	#include <stats.h>			<Header>
	float rms(float in[], int length);	<Prototype>

*/

#include "lib_glob.h"
#include "sts_glob.h"

.SEGMENT/CODE	Code_Space_Name;
.FILE 		RTL_FILENAME;

.GLOBAL	    	_rms;
.GLOBAL		___rmsD, ___rmsP;

_rms:
___rmsD:	JUMP (PC, rms_core) (DB);
		pm_mdf=0;
		pm_ptr=rmsd;

___rmsP:	JUMP (PC, rms_core) (DB);
		pm_mdf=0;
		pm_ptr=rmsp;

rms_core:	CALL (pm_mdf, pm_ptr) (DB);	/* Call routine to compute sum	*/
		put(R7);
		put(R11);
	    
		CALL (PC, ___float_divide);	/* Setup done in subroutine	*/

		F2=0.5;

		F4=RSQRTS F7;			/*Fetch seed*/
		F12=F4*F4, F11=F7;		/*F12=X0^2*/
		F12=F12*F11;			/*F12=C*X0^2*/
		F4=F2*F4, F12=F8-F12;		/*F4=.5*X0, F10=3-C*X0^2*/
		F4=F4*F12;			/*F4=X1=.5*X0(3-C*X0^2)*/
		F12=F4*F4;			/*F12=X1^2*/
		F12=F12*F11;			/*F12=C*X1^2*/
		F4=F2*F4, F12=F8-F12;		/*F4=.5*X1, F10=3-C*X1^2*/
		F4=F4*F12;			/*F4=X2=.5*X1(3-C*X1^2)*/
		F12=F4*F4;			/*F12=X2^2*/
		F12=F12*F11;			/*F12=C*X2^2*/
		F4=F2*F4, F12=F8-F12;		/*F4=.5*X2, F10=3-C*X2^2*/
		F4=F4*F12;			/*F4=X3=.5*X2(3-C*X2^2)*/
		F0=F4*F11, FETCH_RETURN;	/*X=sqrt(Y)=Y/sqrt(Y)*/

restore_state:	get(R11,1);
		get(R7,2);
		RETURN (DB);
		RESTORE_STACK
		RESTORE_FRAME

error_return:	JUMP (PC, restore_state) (DB);
		R0=0;
		NOP;			    	/*Set error if desired		*/

rmsd:		dm_ptr=R4;			/*Point to input		*/
		F12=FLOAT R8, dm_lnt=dm_0;	/*Convert to FP			*/
		R7=R7-R7, R2=R8;		/*Set first input to 0		*/

		R2=R2-1, F0=DM(dm_ptr,dm_1);
		IF LE JUMP (PC, error_return);	/*Must have > 1 point*/
		LCNTR=R2, DO sum_sqrsd UNTIL LCE;
		    F8=F0*F0, F0=DM(dm_ptr,dm_1);
sum_sqrsd:	    F7=F8+F7;
		RTS (DB);
		F11=2.0;
		F8=3.0;

rmsp:		pm_ptr=R4;			/*Point to input		*/
		F12=FLOAT R8, pm_lnt=pm_0;	/*Convert to FP			*/
		R7=R7-R7, R2=R8;		/*Set first input to 0		*/

		R2=R2-1, F0=PM(pm_ptr,pm_1);
		IF LE JUMP (PC, error_return);	/*Must have > 1 point*/
		LCNTR=R2, DO sum_sqrsp UNTIL LCE;
		    F8=F0*F0, F0=PM(pm_ptr,pm_1);
sum_sqrsp:	    F7=F8+F7;
		RTS (DB);
		F11=2.0;
		F8=3.0;

.ENDSEG;