/*
	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
	Updated on 5/94 by AS

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

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

.SEGMENT/CODE	Code_Space_Name;
.FILE 		RTL_FILENAME;

.GLOBAL		_sqrtf;

_sqrtf:		F0=PASS F4, R1=MODE1;		/*Tst for neg, hold old mode*/
		IF LT JUMP (PC, domain_error);

		IF EQ JUMP (PC, restore_state) (DB); /* see rsqrts as to why*/
		FETCH_RETURN
		F8=3.0;
		
		F2=0.5;
		BIT CLR MODE1 65536;		/*Set to 40-BIT mode*/

		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*/
		IF LT JUMP (PC, __big_input);
		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*/
		IF LT JUMP (PC, __big_input);
		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*/
		IF LT JUMP (PC, __big_input);
		F4=F4*F12;			/*F4=X3=.5*X2(3-C*X2^2)*/
		F0=F4*F0, MODE1=R1;		/*X=sqrt(Y)=Y/sqrt(Y)*/

restore_state:	RETURN (DB);
		RESTORE_STACK
                RESTORE_FRAME

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

__big_input:	JUMP (PC, restore_state) (DB);
		FETCH_RETURN
		F0=F0*F4, MODE1=R1;
		
.ENDSEG;