/*
 *  sqrt.c
 *
 *  AMENDMENT HISTORY
 *  ~~~~~~~~~~~~~~~~~
 *  30 Dec 93   DJW   - No logic change - but changed to use 'float' variable
 *                      and functions instead of 'double' ones.
 */

#define __LIBRARY__

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

/* initial minimum mean-squared-error linear approximation constants */
/* N.B. will be rounded to 'float' accuracy bu=y compiler */
#define C1  0.582905362174731941742E0
#define C0  0.424749790911252694093E0

float sqrtf _LIB_F1_(float, x)
{
    int exponent;
    float y;

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

    if(x == 0)
        return 0;

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

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

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

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

    if(exponent & 1)
        y *= (float)M_SQRT2;

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