/*	qsqrta.c		*/
/* Square root check routine, done by long division. */
/* Copyright (C) 1984-1988 by Stephen L. Moshier. */


#include "qhead.h"
#include "mconf.h"

extern QELT qzero[], qone[];
static QELT rndbit[NQ+1] = {0};
static QELT rndset = 0;
#ifdef ANSIPROT
int addm( QELT *, QELT * );
int subm( QELT *, QELT * );
int cmpm( QELT *, QELT * );
int qmovz( QELT *, QELT * );
int shup1( QELT * ), shdn1( QELT * );
int mdnorm( QELT *, int );
#else
int mtherr(), qmovz(), shdn1(), shup1(), addm(), cmpm(), subm(), mdnorm();
#endif

int qsqrt( x, y )
QELT *x, *y;
{
QELT temp[NQ+1], num[NQ+1], sq[NQ+1], xx[NQ+1];
int i, j, k, n;
long m, exp;

if( rndset == 0 )
	{
	qclear( rndbit );
	rndbit[NQ] = 0;
	rndbit[NQ-1] = 1;
	rndset = 1;
	}
/* Check for arg <= 0 */
if( x[0] != 0 )
	{
	if( qcmp( qzero, x ) == 0 )
		{
		qclear(y);
		qneg(y);
		}
	else
		{
		mtherr( "qsqrt", DOMAIN );
		qclear(y);
		}
	return 0;
	}
if( x[1] == 0 )
	{
	qclear(y);
	return 0;
	}

/* Bring in the arg and renormalize if it is denormal. */
qmovz( x, xx );
m = xx[1]; /* local long word exponent */
/*
if( m == 0 )
	m -= normlz( xx );
*/
/* Divide exponent by 2 */
m -= EXPONE - 1;
exp = (QELT) ( (m / 2) + (EXPONE - 1) );

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

qclear( sq );
sq[NQ] = 0;
qclear( num );
num[NQ] = 0;
n = WORDSIZE / 2;

for( j=0; j<2*NQ-6; j++ )
	{
/* bring in next word of arg */
	if( j <= NQ-3 )
		num[NQ] = xx[j+3];
/* Do additional bit on last outer loop, for roundoff. */
	if( j == 2*NQ-7 )
		n = WORDSIZE / 2  + 1;
	for( i=0; i<n; i++ )
		{
/* Next 2 bits of arg */
		shup1( num );
		shup1( num );
/* Shift up answer */
		shup1( sq );
/* Make trial divisor */
		for( k=0; k<=NQ; k++ )
			temp[k] = sq[k];
		shup1( temp );
		addm( rndbit, temp );
/* Subtract and insert answer bit if it goes in */
		if( cmpm( temp, num ) <= 0 )
			{
			subm( temp, num );
			sq[NQ-1] |= 1;
			}
		}
	}

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

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

/* Renormalize and round off. */
mdnorm( sq, k );
qmov( sq, y );
return 0;
}