/****************************************************************************/
/*  SQRT   v3.70                                                           */
/*  Copyright (c) 1993-2001 Texas Instruments Incorporated                  */
/****************************************************************************/
#include <math.h>
#include <errno.h>
#include "values.h"

/***************************************************************************/
/* SQRT() - Square Root                                                    */
/*   Computes square root of x using a Newton-Raphson approximation for    */
/*   sqrt(1/x).  Initial value x0 = .75 * 2 ^ -(e/2), where x = a * 2 ^ e. */
/*   This is the algorithm from page 11-30 of the TMS320C3x User's Guide.  */
/***************************************************************************/
#define ITERATIONS 5

double sqrt(double x)
{
    double x0;          /* estimate */
    int exp;
    int i;

    /************************************************************************/
    /* Check to see if the input is not in the function's domain.           */
    /************************************************************************/
    if (x <= 0.0)
    {
      if (x < 0.0) errno = EDOM;
      return (0.0);
    }

    /************************************************************************/
    /* initial estimate = .75 * 2 ^ -(exp/2)                                */
    /************************************************************************/
    exp = ( (*((unsigned long *) &x) >> 23) & 0xFF) - 0x7F;
    x0  = ldexp(0.75, -exp / 2);

    /************************************************************************/
    /* Refine estimate                                                      */
    /************************************************************************/
    i = ITERATIONS;
    do
      x0 *= (1.5 - x * 0.5 * x0 * x0);
    while (--i);

    /************************************************************************/
    /* sqrt(x) = x * sqrt(1/x)                                              */
    /************************************************************************/
    return (x0 * x);
}