/*
	This C language subroutine computes the 
	reciprocal square root of its double 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>
	double rsqrt(double x);			<Prototype>
*/

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

.SEGMENT/CODE	Code_Space_Name;
.FILE 		RTL_FILENAME;

.GLOBAL		_rsqrt;

/*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.*/

_rsqrt:		dm_ptr=MODE1;
		BIT CLR MODE1 65536;		/*Set 40-BIT mode*/
		CALL (PC, ___lib_dtof) (DB);		/*convert dbl to float*/
		R0=PASS R4, put(R3);
		R1=PASS R8;

		F0=PASS F0, FETCH_RETURN
		IF EQ JUMP (PC, restore_state);
		F8=3.0;
		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, F12=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, F12=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*/		
		CALL (PC, ___lib_ftod) (DB);
		F4=F2*F4, F12=F8-F12;		/*F4=.5*X2, F12=3-C*X2^2*/
		F0=F4*F12, MODE1=dm_ptr;	/*F4=X3=.5*X2(3-C*X2^2)*/

restore_state:	get(R3,1);
		RETURN (DB);
		RESTORE_STACK
		RESTORE_FRAME
.ENDSEG;