/* int	_LnDec80Bit(pDst,pSrc)
 * ARGUMENT
 *	pDst is	a DEC pointer to the destination.
 *	pSrc is	a DEC pointer to the source.
 *
 * DESCRIPTION
 *	Takes the natural logarithm of pSrc, storing result in pDst.
 *
 * SIDE	EFFECTS
 *	Possible errors:
 *	GM_NULLPOINTER
 *	GM_IMAG	if pSrc	<= 0
 *	GM_OVERFLOW if pSrc >= 0 and overflow occurs
 *
 * RETURNS
 *	GM_SUCCESS if successful, error	code otherwise.
 *
 * ALGORTIHM
 *    pSrc = 10^m * 2^(-d) * x		.7 <= x	< 1.4
 *    z	= (x - 1) / (x + 1)
 *    ln pSrc =	m * ln 10 - d *	ln 2 + ln x
 *    ln x = 2 * SUM [z^(2 * n + 1) / (2 * n + 1)]
 *	n = 0, 1, 2, 3,...
 *
 * AUTHOR
 *  Jared Levy		4 / 8 /	1987
 *   Copyright (C) 1987-1990 Greenleaf Software	Inc.  All rights reserved.
 *
 * MODIFICATIONS
 *
 */

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

int	_LnDec80Bit(pDst,pSrc)
DEC	*pDst, *pSrc;
{

	int	m;
	DEC	*pm, dm, *nsrc,	dnsrc, *one, *num, dnum;
	DEC	*den, dden, *z,	dz, *zsq, dzsq,	*pdt, dterm, *fact, dfact;
	DEC	*powz, dpowz;

	if (!_MacIsDecP(pSrc))
		return(GM_IMAG);

/* find	order of magnitude and normalize to range .1 < nsrc < 1	*/
	m = MagnitudeOfDecimal(pSrc) + 1;
	nsrc = &dnsrc;
	_MacDCopy(nsrc,	pSrc);
	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;
	}
/* calculate m * ln(10)	*/
	pm = ConvLongToDecimal(&dm, (long) m);
	(void) _MulDec80Bit(pDst, &decLn10, pm);

/* normalize to	.7 < nsrc <= 1.4,  adjust by factors of	log 2 */
	while (CompareDecimal(nsrc, &decPoint7)	<= 0)  {
			/* double nsrc */
		(void) _AddDec80Bit(nsrc, nsrc,	nsrc);
			/* subtract a ln 2 */
		(void) _SubDec80Bit(pDst, pDst,	&decLn2);
		}

/* calculate z = (nsrc - 1) / (nsrc + 1) */
	one = &decOne;
	num = &dnum;
	(void) _SubDec80Bit(&dnum, nsrc, one);
	den = &dden;
	(void) _AddDec80Bit(&dden, nsrc, one);
	z = &dz;
	(void) _DivRndDec80Bit(&dz, num, den, 23);

/* calculate log nsrc as an expansion of z */
	zsq = &dzsq;
	(void) _MulDec80Bit(&dzsq, z, z);
	fact = ConvLongToDecimal(&dfact, 3L);
	powz = &dpowz;
	/* first term is 2 * z */
	(void) _AddDec80Bit(powz, z, z);
	(void) _AddDec80Bit(pDst, pDst,	powz);
	pdt = &dterm;
	do {
			/* powz	= z^(2n+1) */
		(void) _MulDec80Bit(powz, powz,	zsq);
			/* term	= powz/(2n+1) */
		(void) _DivRndDec80Bit(&dterm, powz, fact,23);
			/* then	add new	term */
		(void) _AddDec80Bit(pDst, pDst,	pdt);
		/* and increase	factorial count	*/
		fact->dc.sl[0] += 2;
		}  while (!(_MacIsDecZ(pdt)));

	return(GM_SUCCESS);
}