/************************************************************************
 *									*
 *				N O T I C E				*
 *									*
 *			Copyright Abandoned, 1987, Fred Fish		*
 *									*
 *	This previously copyrighted work has been placed into the	*
 *	public domain by the author (Fred Fish) and may be freely used	*
 *	for any purpose, private or commercial.  I would appreciate	*
 *	it, as a courtesy, if this notice is left in all copies and	*
 *	derivative works.  Thank you, and enjoy...			*
 *									*
 *	The author makes no warranty of any kind with respect to this	*
 *	product and explicitly disclaims any implied warranties of	*
 *	merchantability or fitness for any particular purpose.		*
 *									*
 ************************************************************************
 */


/*
 *  FUNCTION
 *
 *	sqrt   double precision square root
 *
 *  KEY WORDS
 *
 *	sqrt
 *	machine independent routines
 *	math libraries
 *
 *  DESCRIPTION
 *
 *	Returns double precision square root of double precision
 *	floating point argument.
 *
 *  USAGE
 *
 *	double sqrt (x)
 *	double x;
 *
 *  REFERENCES
 *
 *	Fortran IV-PLUS user's guide, Digital Equipment Corp. pp B-7.
 *
 *	Computer Approximations, J.F. Hart et al, John Wiley & Sons,
 *	1968, pp. 89-96.
 *
 *  RESTRICTIONS
 *
 *	The relative error is 10**(-30.1) after three applications
 *	of Heron's iteration for the square root.
 *
 *	However, this assumes exact arithmetic in the iterations
 *	and initial approximation.  Additional errors may occur
 *	due to truncation, rounding, or machine precision limits.
 *	
 *  PROGRAMMER
 *
 *	Fred Fish
 *
 *  INTERNALS
 *
 *	Computes square root by:
 *
 *	  (1)	Range reduction of argument to [0.5,1.0]
 *		by application of identity:
 *
 *		sqrt(x)  =  2**(k/2) * sqrt(x * 2**(-k))
 *
 *		k is the exponent when x is written as
 *		a mantissa times a power of 2 (m * 2**k).
 *		It is assumed that the mantissa is
 *		already normalized (0.5 =< m < 1.0).
 *
 *	  (2)	An approximation to sqrt(m) is obtained
 *		from:
 *
 *		u = sqrt(m) = (P0 + P1*m) / (Q0 + Q1*m)
 *
 *		P0 = 0.594604482
 *		P1 = 2.54164041
 *		Q0 = 2.13725758
 *		Q1 = 1.0
 *
 *		(coefficients from HART table #350 pg 193)
 *
 *	  (3)	Three applications of Heron's iteration are
 *		performed using:
 *
 *		y[n+1] = 0.5 * (y[n] + (m/y[n]))
 *
 *		where y[0] = u = approx sqrt(m)
 *
 *	  (4)	If the value of k was odd then y is either
 *		multiplied by the square root of two or
 *		divided by the square root of two for k positive
 *		or negative respectively.  This rescales y
 *		by multiplying by 2**frac(k/2).
 *
 *	  (5)	Finally, y is rescaled by int(k/2) which
 *		is equivalent to multiplication by 2**int(k/2).
 *
 *		The result of steps 4 and 5 is that the value
 *		of y between 0.5 and 1.0 has been rescaled by
 *		2**(k/2) which removes the original rescaling
 *		done prior to finding the mantissa square root.
 *
 *  NOTES
 *
 *	The Convergent Technologies compiler optimizes division
 *	by powers of two to "arithmetic shift right" instructions.
 *	This is ok for positive integers but gives different
 *	results than other C compilers for negative integers.
 *	For example, (-1)/2 is -1 on the CT box, 0 on every other
 *	machine I tried.
 *
 */
/*
 *  MODIFICATIONS
 *
 *	This routine modified to use polynomial, instead of rational,
 *	approximation with coefficients  
 *
 *		P0  0.22906994529e+00
 *		P1  0.13006690496e+01
 *		P2 -0.90932104982e+00
 *		P3  0.50104207633e+00
 *		P4 -0.12146838249e+00
 *
 *	as given by Hart (op. cit.) in SQRT 0132.  This approximation
 *	gives around 5 digits correct.  Two applications of Heron's
 *	approximation will give more then practically achievable
 *	13-15 decimal digits.  Multiplications by powers of 2 are
 *	replaced by explicit calls to ldexp.
 *
 *	Michal Jaegermann, 24 October 1990
 */
 
#include <stdio.h>
#include <pmluser.h>
#include "pml.h"
 
#ifdef OLD
#define P0 0.594604482			/* Approximation coeff	*/
#define P1 2.54164041			/* Approximation coeff	*/
#define Q0 2.13725758			/* Approximation coeff	*/
#define Q1 1.0				/* Approximation coeff	*/

#define ITERATIONS 3			/* Number of iterations	*/

#endif

#define P0  0.22906994529e+00		/* Hart SQRT 0132 */
#define P1  0.13006690496e+01
#define P2 -0.90932104982e+00
#define P3  0.50104207633e+00
#define P4 -0.12146838249e+00

static char funcname[] = "sqrt";


double sqrt (x)
double x;
{
#ifdef OLD
    auto int k;
    register int bugfix;
    register int kmod2;
    register int count;
    auto int exponent;
#else
#ifdef __MSHORT__
    auto short k;
#else
    auto int k;
#endif
#endif
    auto double m;
    auto double u;
#ifdef OLD
    auto double y;
#endif
    auto struct exception xcpt;
    extern double frexp ();
    extern double ldexp ();
 
    DBUG_ENTER ("sqrt");
    DBUG_3 ("sqrtin", "arg %le", x);
    if (x == 0.0) {
	xcpt.retval = 0.0;
    } else if (x < 0.0) {
	xcpt.type = DOMAIN;
	xcpt.name = funcname;
	xcpt.arg1 = x;
	if (!matherr (&xcpt)) {
	    fprintf (stderr, "%s: DOMAIN error\n", funcname);
	    errno = EDOM;
	    xcpt.retval = 0.0;
	}
    } else {
	m = frexp (x, &k);
#ifdef OLD
	u = (P0 + (P1 * m)) / (Q0 + (Q1 * m));
	for (count = 0, y = u; count < ITERATIONS; count++) {
	    y = 0.5 * (y + (m / y));
	}
	if ((kmod2 = (k % 2)) < 0) {
	    y /= SQRT2;
	} else if (kmod2 > 0) {
	    y *= SQRT2;
	}
	bugfix = 2;
	xcpt.retval = ldexp (y, k/bugfix);
#else
	u = (((P4 * m + P3) * m + P2) * m + P1) * m + P0;
	u = ldexp((u + (m / u)), -1);	/* Heron's iteration */
	u += m / u;			/* and a part of the second one */
	if (k & 1) {
	    u *= SQRT2;
	}
	/* 
	 * here we rely on the fact that -3/2 and (-3 >> 1)
	 * do give different results
	 */
	xcpt.retval = ldexp (u, (k >> 1) - 1);
#endif
    }
    DBUG_3 ("sqrtout", "result %le", xcpt.retval);
    DBUG_RETURN (xcpt.retval);
}