/*
	This C language subroutine computes the 
	reciprocal square root of its floating point input.

	The Run Time Library for the C Language.

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

	Created on 10/18/90
	Updated on 5/94 by AS

	#include <math.h>		<Header>
	float rsqrtf(float x);		<Prototype>
*/

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

.SEGMENT/CODE	Code_Space_Name;
.FILE RTL_FILENAME;

.GLOBAL		_rsqrtf;

/*Based on a customer complaint.  When x=0, RSQRTS sets flags off b'cse it is
trying to do a 1/(sqrt(x)) which becomes 1/0.  According to IEEE standards,
'0' is part of the domain of numbers for a sqrt to handle an b'cse of this
no errors should be reported b'cse of a 0 input.  A check has been added at
the top to take care of this.*/

_rsqrtf:	R0=PASS R4, R1=MODE1;
		IF EQ JUMP (PC, restore_state) (DB);
		F8=3.0;
		FETCH_RETURN

		BIT CLR MODE1 65536;		/*Set to 40-BIT mode*/
		F2=0.5;

		F4=RSQRTS F0;			/*Fetch seed*/
		F12=F4*F4;			/*F12=X0^2*/
		F12=F12*F0;			/*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*F0;			/*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*F0;			/*F12=C*X2^2*/
		F4=F2*F4, F12=F8-F12;		/*F4=.5*X2, F10=3-C*X2^2*/
		F4=F4*F12, MODE1=R1;		/*F4=X3=.5*X2(3-C*X2^2)*/

		F0=RND F4; 

restore_state:	RETURN (DB);
		RESTORE_STACK
		RESTORE_FRAME
.ENDSEG;