/* void	_SqrtDec80Bit(pDst,pSrc)
 *
 * ARGUMENT
 *	DEC	*pDst;
 *	DEC	*pSrc;
 *
 * DESCRIPTION
 *	Sets pDst = square root	of pSrc;
 *
 * SIDE	EFFECTS
 *	None.
 *
 * RETURNS
 *	None.  Always succeeds because CALLER GUARANTEES a valid 64-bit
 *  value for square root operation.
 *
 * AUTHOR
 *  Jared Levy		Aug 7, 1987
 *   Copyright (C) 1987-1990 Greenleaf Software	Inc.  All rights reserved.
 *
 * MODIFICATIONS
 *
 */

#include <stdio.h>
#include "gm.h"
#include "gmsystem.h"

void	_SqrtDec80Bit(pDst,pSrc)
DEC	*pDst, *pSrc;
{
	int	m, precision;
	DEC	*temp, dtemp, *nsrc, dnsrc;

	if(_MacIsDecZ(pSrc))  {
		_MacDZero(pDst);
		return;
	}

	temp = &dtemp;
	nsrc = &dnsrc;

/* adjust pSrc to range	.01 < nsrc < 1 */
	_MacDCopy(nsrc,pSrc);
	m = MagnitudeOfDecimal(nsrc);

/* calculate first guess */
	if ((m & 0X0001) == 0)	{
		m+=2;
		nsrc->dc.id+= m;
		if (nsrc->dc.id	> GM_IMAXID)  {
		   _DivUnsArrByPwrOf10(nsrc->dc.sl, 5, nsrc->dc.id-GM_IMAXID);
		   nsrc->dc.id = GM_IMAXID;
		    }
			/* decTwoPoint59 = 2.59	*/
		(void) _MulDec80Bit(temp, &decTwoPoint59, nsrc);
			/* decPoint0819	 =  .0819 */
		(void) _AddDec80Bit(pDst, &decPoint0819, temp);	 }
	else {
		m++;
		nsrc->dc.id+= m;
			/* decPoint819 = .819 */
		(void) _MulDec80Bit(temp, &decPoint819,	nsrc);
			/* decPoint259 = .259 */
		(void) _AddDec80Bit(pDst, &decPoint259,	temp);	}

/* do iterative	calculation */
	precision = 3;
	(void) _ScaleDec80Bit(pDst, pDst, 4);
	while (precision < 18)	{
		precision = 2 *	precision - 2;
		(void) _DivRndDec80Bit(temp, nsrc, pDst, precision);
		(void) _AddDec80Bit(pDst, pDst,	temp);
		_HalveUnsArr(pDst->dc.sl, 5);
	}

/* last	iteration */
	(void) _DivDec80Bit(temp, nsrc,	pDst);
	(void) _AddDec80Bit(pDst, pDst,	temp);
	if (((pDst->dc.sl[0] & 1) == 0)	|| (pDst->dc.msd > 6553))
		_HalveUnsArr(pDst->dc.sl, 5);
	else
		(void) _MulDec80Bit(pDst, pDst,	&decPoint5);

/* recalibrate to correct magnitude */
	pDst->dc.id -= (m/2);

}