/*	esqrt.c		*/
/* Square root check routine, done by long division. */
/* by Stephen L. Moshier. */


#include "mconf.h"
#include "ehead.h"

void emovi(), emovo(), ecleaz(), emdnorm(), eaddm(), esubm();
int ecmpm();
static unsigned short rndbit[NI] = {0,0,0,0,0,0,0,1,0};

void esqrt( x, y )
short *x, *y;
{
unsigned short temp[NI], num[NI], sq[NI], xx[NI];
int i, j, k, n;
long m, exp;

/* Check for arg <= 0 */
i = ecmp( x, ezero );
if( i <= 0 )
	{
	eclear(y);
	if( i < 0 )
		mtherr( "esqrt", DOMAIN );
	return;
	}

/* Bring in the arg and renormalize if it is denormal. */
emovi( x, xx );
m = (long )xx[1]; /* local long word exponent */
if( m == 0 )
	m -= enormlz( xx );

/* Divide exponent by 2 */
m -= 0x3ffe;
exp = (unsigned short )( (m / 2) + 0x3ffe );

/* Adjust if exponent odd */
if( (m & 1) != 0 )
	{
	if( m > 0 )
		exp += 1;
	eshdn1( xx );
	}

ecleaz( sq );
ecleaz( num );
n = 8; /* get 8 bits of result per inner loop */

for( j=0; j<(2*NE-2); j++ )
	{
/* bring in next word of arg */
	if( j < NE )
		num[NI-1] = xx[j+3];
/* Do additional bit on last outer loop, for roundoff. */
	if( j == (2*NE-3) )
		n = 9;
	for( i=0; i<n; i++ )
		{
/* Next 2 bits of arg */
		eshup1( num );
		eshup1( num );
/* Shift up answer */
		eshup1( sq );
/* Make trial divisor */
		for( k=0; k<NI; k++ )
			temp[k] = sq[k];
		eshup1( temp );
		eaddm( rndbit, temp );
/* Subtract and insert answer bit if it goes in */
		if( ecmpm( temp, num ) <= 0 )
			{
			esubm( temp, num );
			sq[NI-2] |= 1;
			}
		}
	}

exp -= 1; /* Adjust for extra, roundoff loop done. */

/* Sticky bit = 1 if the remainder is nonzero. */
k = 0;
for( i=3; i<NI; i++ )
	k |= (int )num[i];

/* Renormalize and round off. */
emdnorm( sq, k, 0, exp, 64 );
emovo( sq, y );
}