/*
   FUNCTION:  double sqrt (x)

   PURPOSE:
	To compute the square root of a double precision floating point number

   INPUT PARAMETERS:
	double x - The given number

   OUTPUT PARAMETERS:
	Returned - SQRT(x) if x is non-negative.  If x is negative, x is 
		   returned unchanged.

   DESCRIPTION:
	The algorithm is standard:  make an initial guess by halving the 
	exponent of the number.  Then use Newton's method to converge on the
	correct value

   WRITTEN BY:  Scott D. Roth		DATE:  March 27, 1981
*/

double sqrt (x)
    double x;
	{
	double root;
	register int i, *ptr;

	if ((root = x) > 0) {
	    ptr = &root;		/* Halve the exponent--first guess */
	    *ptr = (*ptr >> 1) + 020100;
	    for (i = 0; i < 4; ++i) {	/* Do four Newton iterations */
		root = x/root + root;
		if (*ptr != 0)		/* y /= 2	*/
		    *ptr -= 0200;
		}
	    }
	return(root);
	}