/*							sqrtl.c
 *
 *	Square root, long double precision
 *
 *
 *
 * SYNOPSIS:
 *
 * long double x, y, sqrtl();
 *
 * y = sqrtl( x );
 *
 *
 *
 * DESCRIPTION:
 *
 * Returns the square root of x.
 *
 * Range reduction involves isolating the power of two of the
 * argument and using a polynomial approximation to obtain
 * a rough value for the square root.  Then Heron's iteration
 * is used three times to converge to an accurate value.
 *
 * Note, some arithmetic coprocessors such as the 8087 and
 * 68881 produce correctly rounded square roots, which this
 * routine will not.
 *
 * ACCURACY:
 *
 *
 *                      Relative error:
 * arithmetic   domain     # trials      peak         rms
 *    IEEE      0,10        30000       8.1e-20     3.1e-20
 *
 *
 * ERROR MESSAGES:
 *
 *   message         condition      value returned
 * sqrt domain        x < 0            0.0
 *
 */

/*
Cephes Math Library Release 2.2:  December, 1990
Copyright 1984, 1990 by Stephen L. Moshier
Direct inquiries to 30 Frost Street, Cambridge, MA 02140
*/


#include <math.h>

#define SQRT2 1.4142135623730950488017E0L
#ifdef ANSIPROT
extern long double frexpl ( long double, int * );
extern long double ldexpl ( long double, int );
#else
long double frexpl(), ldexpl();
#endif

long double sqrtl(x)
long double x;
{
int e;
long double z, w;
#ifndef UNK
short *q;
#endif

if( x <= 0.0 )
	{
	if( x < 0.0 )
		mtherr( "sqrtl", DOMAIN );
	return( 0.0 );
	}
w = x;
/* separate exponent and significand */
#ifdef UNK
z = frexpl( x, &e );
#endif

/* Note, frexp and ldexp are used in order to
 * handle denormal numbers properly.
 */
#ifdef IBMPC
z = frexpl( x, &e );
q = (short *)&x; /* point to the exponent word */
q += 4;
/*
e = ((*q >> 4) & 0x0fff) - 0x3fe;
*q &= 0x000f;
*q |= 0x3fe0;
z = x;
*/
#endif
#ifdef MIEEE
z = frexpl( x, &e );
q = (short *)&x;
/*
e = ((*q >> 4) & 0x0fff) - 0x3fe;
*q &= 0x000f;
*q |= 0x3fe0;
z = x;
*/
#endif

/* approximate square root of number between 0.5 and 1
 * relative error of linear approximation = 7.47e-3
 */
/*
x = 0.4173075996388649989089L + 0.59016206709064458299663L * z;
*/

/* quadratic approximation, relative error 6.45e-4 */
x = ( -0.20440583154734771959904L  * z
     + 0.89019407351052789754347L) * z
     + 0.31356706742295303132394L;

/* adjust for odd powers of 2 */
if( (e & 1) != 0 )
	x *= SQRT2;

/* re-insert exponent */
#ifdef UNK
x = ldexpl( x, (e >> 1) );
#endif
#ifdef IBMPC
x = ldexpl( x, (e >> 1) );
/*
*q += ((e >>1) & 0x7ff) << 4;
*q &= 077777;
*/
#endif
#ifdef MIEEE
x = ldexpl( x, (e >> 1) );
/*
*q += ((e >>1) & 0x7ff) << 4;
*q &= 077777;
*/
#endif

/* Newton iterations: */
#ifdef UNK
x += w/x;
x = ldexpl( x, -1 );	/* divide by 2 */
x += w/x;
x = ldexpl( x, -1 );
x += w/x;
x = ldexpl( x, -1 );
#endif

/* Note, assume the square root cannot be denormal,
 * so it is safe to use integer exponent operations here.
 */
#ifdef IBMPC
x += w/x;
*q -= 1;
x += w/x;
*q -= 1;
x += w/x;
*q -= 1;
#endif
#ifdef MIEEE
x += w/x;
*q -= 1;
x += w/x;
*q -= 1;
x += w/x;
*q -= 1;
#endif

return(x);
}