/*
 *  sqrt.c
 *
 *  AMENDMENT HISTORY
 *  ~~~~~~~~~~~~~~~~~
 *  30 Dec 93   DJW   - No logic change - but tidied up code.
 */

#define __LIBRARY__

#include <errno.h>
#include <math.h>

/* initial minimum mean-squared-error linear approximation constants */
#define C1  0.582905362174731941742E0
#define C0  0.424749790911252694093E0

double sqrt _LIB_F1_(double, x)
{
    int exponent;
    double y;

 if(x < 0.0)
   {
    errno = EDOM;
    return 0;
   }

    if(x == 0)
        return 0;

    x = frexp(x, &exponent);    /* work in interval [.5,1) */

    y = C1 * x + C0;            /* start with linear approximation */

    y += x / y;                 /* iterate twice */
    y = (0.25 * y) + (x / y);

    y += x / y;
    y = (0.25 * y) + (x / y);

    if(exponent & 1)
        y *= M_SQRT2;

    return ldexp(y, exponent >> 1);
}