/* void	_ATanDec80Bit(pDst,pSrc)
 *
 * ARGUMENT
 *	DEC	*pDst;
 *	DEC	*pSrc;
 *
 * DESCRIPTION
 *	Sets pDst = the	arctangent (in radians)	of pSrc,
 *  which always lies between -pi/2 and	pi/2.
 *
 * SIDE	EFFECTS
 *	None.
 *
 * RETURNS
 *	None. (always successful)
 *
 * AUTHOR
 *  Jared Levy		Aug 7, 1987
 *   Copyright (C) 1987-1990 Greenleaf Software	Inc.  All rights reserved.
 *
 * ALGORITHM
 *	For x <	0, atan(x) = - atan(-x)
 *	For x >	1  atan(x) = pi	/ 2 - acot(x) =	pi / 2 - atan(1	/ x)
 *	atan(x)	= pi / 6 + atan	[ (x - 1/3^.5) / (1 + x/3^.5) ]
 *	atan(x)	= x - x^3/3 + x^5/5 - ...
 *
 * MODIFICATIONS
 *
 */

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

void	_ATanDec80Bit(pDst,pSrc)
DEC	*pDst;
DEC	*pSrc;
{
	int	isn=0, isrecip=0;
	DEC	*nsrc, dnsrc, *x, dx, *xsq, dxsq, *xpow, dxpow;
	DEC	*term, dterm, *fact, dfact;
	DEC	*num, dnum, *den, dden,	*temp, dtemp;

	nsrc = &dnsrc;
	_MacDCopy(nsrc,	pSrc);
	if (_MacIsDecN(nsrc))  {
		isn = 1;
		_MacDChgs(nsrc);
		}

	if (CompareDecimal(nsrc,&decOne) == 1)	{
		isrecip	= 1;
		(void) _DivRndDec80Bit(nsrc, &decOne, nsrc, 23);
		}

	(void) _SubDec80Bit(num	= &dnum, nsrc, &decReciprocalOfSqrt3);
	(void) _MulDec80Bit(temp = &dtemp, &decReciprocalOfSqrt3, nsrc);
	(void) _AddDec80Bit(den	= &dden, temp, &decOne);
	(void) _DivRndDec80Bit(x = &dx,	num, den, 23);

	(void) _MulDec80Bit(xsq	= &dxsq, x , x);
	(void) _MulDec80Bit(xpow = &dxpow, xsq,	x);
	_MacDCopy(pDst,	x);
	(void) ConvLongToDecimal(fact =	&dfact,	3L);
	(void) _DivRndDec80Bit(term = &dterm, xpow, fact, 23);
	while (!_MacIsDecZ(term)) {
		if ((fact->dc.sl[0] % 4) == 1)
			(void) _AddDec80Bit(pDst, pDst,	term);
		else
			(void) _SubDec80Bit(pDst, pDst,	term);
		fact->dc.sl[0]+= 2;
		(void) _MulDec80Bit(xpow, xpow,	xsq);
		(void) _DivRndDec80Bit(term, xpow, fact, 23);
		}

	(void) _AddDec80Bit(pDst, &decPiOver6, pDst);

	if (isrecip)
		(void) _SubDec80Bit(pDst, &decPiOver2, pDst);

	if (isn)
		_MacDChgs(pDst);
}