/*
	This C language subroutine computes the 
	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 4/28/90

	#include <math.h>			<Header>
	double sqrt(double x);			<Prototype>
*/

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

.SEGMENT/CODE	Code_Space_Name;
.FILE 		RTL_FILENAME;

.GLOBAL		_sqrt;

_sqrt:		put(R3);

		reads(R4,1);
		reads(R8,2);

		R12=MODE1;
		BIT CLR MODE1 65536;		/*Set to 40-BIT mode*/
		CALL (PC, ___lib_dtof) (DB);	/*Convert dbl to float*/
		R0=PASS R4;
		R1=PASS R8;

		F0=PASS F0;			/*Test for negative input*/
		IF LT JUMP (PC, domain_error);	/*Return with 0 and EDOM*/

		IF EQ JUMP (PC, restore_state) (DB); /* see rsqrts for reason */
		FETCH_RETURN
		F8=3.0;
		F2=0.5;

		F4=RSQRTS F0, R1=R12;		/*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, R12=R1;		/*F4=X3=.5*X2(3-C*X2^2)*/
		F4=F4*F0, F2=F0;		/*X=sqrt(Y)=Y/sqrt(Y)*/

cnvrt_to_dbl:	CALL (PC, ___lib_ftod) (DB);		/*Convert to double*/
		F2=PASS F2, F0=F4;		/*Test for zero input*/
		IF EQ F0=PASS F2;		/*Pass zero to result*/

restore_state:	MODE1=R12;			/*Restore old mode*/
		get(R3,1);
		RETURN (DB);
		RESTORE_STACK
		RESTORE_FRAME

domain_error:	ram_ireg=_errno;
		JUMP (PC, cnvrt_to_dbl) (DB);
		rammem(ram_ireg,ram_0)=EDOM;
		R0=R0-R0, FETCH_RETURN
.ENDSEG;