/* DEC	*InternalRateOfReturn(irr, flows, nflow)
 *
 * ARGUMENT
 *	DEC	*irr, **flows;
 *	int	nflow;
 *
 * DESCRIPTION
 *	Calculates the internal	rate of	return of a series of cash flows.
 *   The internal rate of return, rounded to wGMIntrPrec decimal places,
 *   is	stored in irr.
 *
 * SIDE	EFFECTS
 *	None.
 *
 * RETURNS
 *	irr if successful, otherwise GM_NULL.
 *
 * POSSIBLE ERRORS
 *	GM_NULLPOINTER
 *	GM_ARGVAL
 *
 * AUTHOR
 *  Jared Levy
 *   Copyright (C) 1988-1990 Greenleaf Software	Inc.  All rights reserved.
 *
 * MODIFICATIONS
 *
 *
 */

#include <stdio.h>
#include "gmsystem.h"
static void	_NPVaux(DEC *, DEC **, int, DEC	*);

DEC	*InternalRateOfReturn(irr, flows, nflow)
DEC	*irr, **flows;
int	nflow;
{
	mbool	wfPos, wfNeg, wfDone;
	int	i, j, places;
	DEC	dx1, *x1=&dx1, dx2, *x2=&dx2;
	DEC	dy1, *y1=&dy1, dy2, *y2=&dy2;
	DEC	dxn, *xn=&dxn, dyn, *yn=&dyn;
	DEC	dtemp, *temp=&dtemp, dtemp2, *temp2=&dtemp2;

	_MacStart(GM_IRR);

	if (nflow<2)  {
		_MacErr(GM_ARGVAL);
		_MacRet(GM_NULL);
	}

	if (!flows||!irr)  {
		_MacErr(GM_NULLPOINTER);
		_MacRet(GM_NULL);
	}

	/* check for null pointers, at least one pos & neg cash	flow */
	wfPos=FALSE;
	wfNeg=FALSE;
	for (i=0;i<nflow;i++)  {
		_MacInVarD(flows[i]);
		wfPos=wfPos||_MacIsDecP(flows[i]);
		wfNeg=wfNeg||_MacIsDecN(flows[i]);
	}
	if (!wfPos||!wfNeg)  {
		_MacErr(GM_ARGVAL);
		_MacRet(GM_NULL);
	}

	places = wIntrPrec;
	if (places<0)
		places=0;
	if (places>18)
		places=18;

/* find	bounds */
/*
 * The initial 2 guesses for IRR have been changed.  The old code started
 * with	guesses	of 2% &	2.2%, and tried	to extrapolate from there.  I
 * have	changed	it to guess 2% and -2%.
 */
	/* load	initial	guesses	*/
	_MacDCopy(x1, pGMIRRGuess);
	if (_MacIsDecZ(x1))
		(void) MakeDecimalFromLong(x2, 1L, 1);
	_MacDCopy( x2, x1 );
	_MacDChgs( x2 );

	_NPVaux(y1, flows, nflow, x1);
	_NPVaux(y2, flows, nflow, x2);

	/* use secant method to	find bounds */
	i=0;
	while ((!(_MacIsDecN(y1)^_MacIsDecN(y2))) && i<15)  {
		i++;
		j = CompareDecimal(y1, y2);
		if (j==0)  {				/* |y1|=|y2| */
			(void) _AddDec80Bit(x2,	x1, x2);
			(void) _MulDec80Bit(x2,	x2, &decPoint5);
			_NPVaux(y2, flows, nflow, x2);
		}
		else  {
			/* compute new value */
			(void) _SubDec80Bit(temp, x2, x1);
			(void) _SubDec80Bit(temp2, y2, y1);
			(void) _DivDec80Bit(temp, temp,	temp2);
			(void) _MulDec80Bit(temp, temp,	&decOnePoint6);
			if ((j==1)^(_MacIsDecN(y1))) {	/* |y1|>|y2| */
				(void) _MulDec80Bit(temp, temp,	y2);
				(void) _SubDec80Bit(x1,	x2, temp);
				_NPVaux(y1, flows, nflow, x1);
			}
			else  {				/* |y2|>|y1| */
				(void) _MulDec80Bit(temp, temp,	y1);
				(void) _SubDec80Bit(x2,	x1, temp);
				_NPVaux(y2, flows, nflow, x2);
			}
		}
	}

	/* find	bounds,	method 2 */
	if (i>=15)  {
		/* load	initial	guesses	*/
		_MacDCopy(x1, pGMIRRGuess);
		if (_MacIsDecZ(x1))
			MakeDecimalFromLong(x2,	1L, 1);
		else  {
			_MacDCopy(x2, x1);
			x2->dc.id++;
			(void) _SubDec80Bit(x2,	x1, x2);
		}
		_NPVaux(y1, flows, nflow, x1);
		_NPVaux(y2, flows, nflow, x2);

		while (!(_MacIsDecN(y1)^_MacIsDecN(y2)))  {
			j = CompareDecimal(y1, y2);
			if ((j==1)^(_MacIsDecN(y1))) {	/* |y1|>|y2| */
				_AddDec80Bit(temp, x2, x2);
				_SubDec80Bit(x2, temp, x1);
				if ((CompareDecimal(x2,	&decHundred) ==	1)
				 || (CompareDecimal(x2,	&decMinus50) ==	-1)) {
					 _MacErr(GM_ARGVAL);
					 _MacRet(GM_NULL);
					 }
				_NPVaux(y2, flows, nflow, x2);
				 }
			else  {				/* |y2|>|y1| */
				_AddDec80Bit(temp, x1, x1);
				_SubDec80Bit(x1, temp, x2);
				if ((CompareDecimal(x1,	&decHundred) ==	1)
				 || (CompareDecimal(x1,	&decMinus50) ==	-1)) {
					 _MacErr(GM_ARGVAL);
					 _MacRet(GM_NULL);
				 }
				_NPVaux(y1, flows, nflow, x1);
			}
		}
	}


/* If this stage is reached, bounds have been found */
/* Now use bisection to	find the IRR between those bounds */
	wfDone=FALSE;
	while (!wfDone)	 {
/* calculate new position */
		(void) _AddDec80Bit(xn,	x2, x1);
	       _HalveUnsArr(xn->dc.sl, 5);
	       if (xn->dc.id !=	places)
		       _ScaleDec80Bit(xn, xn, places);

/* calculate corresponding y */
		_NPVaux(yn, flows, nflow, xn);

/* check for approximate or exact zero */
		if (_MacIsDecZ(yn))  {
			_MacDCopy(irr,xn);
			_MacRet(irr);
		}

/* adjust appropriate boundary */
		if (_MacIsDecN(y1)^_MacIsDecN(yn))  {
			wfDone=(CompareDecimal(x2,xn)==0);
			_MacDCopy(x2,xn);
			_MacDCopy(y2,yn);
		}
		else  {
			wfDone=(CompareDecimal(x1,xn)==0);
			_MacDCopy(x1,xn);
			_MacDCopy(y1,yn);
		}

	}

	_MacDCopy(irr, xn);
	_MacRet(xn);
}

static void	_NPVaux(net, flows, nflow, intr)
DEC	*net, **flows, *intr;
int	nflow;
{
	int	i;
	DEC	dtemp, *temp=&dtemp, dooopi, *ooopi=&dooopi;
	DEC	dpow, *pow=&dpow;

	_MacDZero(net);
	_MacDCopy(temp,	intr);
	if (temp->dc.id	<= GM_MAXID-2)
		temp->dc.id+=2;
	else
		(void) _DivUnsArrByPwrOf10(temp->dc.sl,	5, 2);
	(void) _AddDec80Bit(temp, temp,	&decOne);
	(void) _DivDec80Bit(ooopi, &decOne, temp);

	_MacDCopy(pow, (&decOne));

	for (i=0;i<nflow;i++)  {
		(void) _MulDec80Bit(temp, flows[i], pow);
		(void) _AddDec80Bit(net, net, temp);
		(void) _MulDec80Bit(pow, pow, ooopi);
	}

}