ack/lang/pc/libpc/sqt.c

/* $Header$ */

#include	<pc_err.h>

extern	double	_fef();
extern		_trp();

/*
	sqrt returns the square root of its floating
	point argument. Newton's method.

	calls _fef
*/

double
_sqt(arg)
double arg;
{
	double x, temp;
	int exp;
	int i;

	if(arg <= 0) {
		if(arg < 0)
			_trp(ESQT);
		return(0);
	}
	x = _fef(arg,&exp);
	/*
	while(x < 0.5) {
		x =* 2;
		exp--;
	}
	*/
	/*
	 * NOTE
	 * this wont work on 1's comp
	 */
	if(exp & 1) {
		x *= 2;
		exp--;
	}
	temp = 0.5*(1 + x);

	while(exp > 28) {
		temp *= (1<<14);
		exp -= 28;
	}
	while(exp < -28) {
		temp /= (1<<14);
		exp += 28;
	}
	if(exp >= 0)
		temp *= 1 << (exp/2);
	else
		temp /= 1 << (-exp/2);
	for(i=0; i<=4; i++)
		temp = 0.5*(temp + arg/temp);
	return(temp);
}