/* DEC	*_InterestAux(fun, intper, fractper, intr, pv, pmt, fv,	begend)
 *
 * ARGUMENT
 *	int	fun, intper, begend;
 *	DEC	*fractper, *intr, *pv, *pmt, *fv;
 *
 * DESCRIPTION
 *	Given four variables involved in compound interest, solves for the
 *  interest.  The variables are the number of periods nper, the percentage
 *  interest rate per period intr, the present value pv, the periodic payment
 *  pmt, and the future	value fv.  begend specifies whether payments take
 *  place at the beginning or the end of each month.
 *	This auxillary function	does the work of CompoundInterest,
 *  CompoundInterestSimple, and	CompoundInterestCompound for calculating intr.
 *
 * SIDE	EFFECTS
 *	Changes	value of intr.
 *
 * RETURNS
 *	In case	of success, intr is returned.
 *  If an error	occurs,	GM_NULL	is returned.
 *
 * POSSIBLE ERRORS
 *	GM_NULLPOINTER
 *	GM_ARGVAL
 *
 *
 *
 * AUTHOR
 *  Jared Levy
 *   Copyright (C) 1988-1990 Greenleaf Software	Inc.  All rights reserved.
 *
 * ALGORITHM
 *	Solve the following formula for	the unknown variable:
 *  0=pv*alpha + (1+intr*begend)*pmt*[1-(1+intr)^-int(nper)]/intr +
 *	 fv*(1+intr)^-int(nper)
 *  where alpha=1 for no odd period,
 *	  alpha=1+intr*frac(nper) for an odd period with simple	interest
 *	  alpha=(1+intr)^frac(nper) for	an odd period with compound interest
 *
 *  The	equation can be	solved in closed form for any variable except intr,
 *  which requires numerical methods.
 *  The	variables used by this routine have the	following meanings:
 *	fun:   1 if called by ComoundInterest, 2 if by CompoundInterestSimple,
 *		3 if by	CompoundInterestCompound, and 0	if by AdvancePayment
 *	intper:	integer	part of	nper
 *	fractper: fractional part of nper (if fun==2 or	fun==3)
 *	intr:  interest	rate per period
 *	pv:    present value
 *	pmt:   period payment
 *	fv:    future value (balloon payment)
 *	begend:	GM_BEGIN if payments at	beginning of period, GM_END if at end
 *		(or if fun==0, number of advanced payments ad)
 *	temp1, temp2: temporary	DEC's
 */

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

DEC	*_InterestAux(fun, intper, fractper, intr, pv, pmt, fv,	begend)
int	fun, intper, begend;
DEC	*fractper, *intr, *pv, *pmt, *fv;
{
	int	i, places, ad;
	mbool	wfDone;
	DEC	dx1, *x1=&dx1, dx2, *x2=&dx2, dpx1, *px1=&dpx1;
	DEC	dpx2, *px2=&dpx2, dy1, *y1=&dy1, dy2, *y2=&dy2;
	DEC	dy0, *y0=&dy0, dxn, *xn=&dxn, dyn, *yn=&dyn;
	DEC	djunk, *junk=&djunk;
	DEC	dtemp1,	*temp1=&dtemp1,	dtemp2,	*temp2=&dtemp2;

/* Can't have zero payments */
	if (intper==0)	{
		_MacErr(GM_ARGVAL);
		return(GM_NULL);
	}

	ad = begend;

/* add advance payments	to pv */
	if (fun==0)  {
		(void) ConvLongToDecimal(junk, (long) ad);
		(void) _MulDec80Bit(junk, pmt, junk);
		(void) _AddDec80Bit(junk, junk,	pv);
		pv = junk;
	}

/* Check that there's a	positive & a negative value among pv, pmt, fv */
	if (!(_MacIsDecP(pv)||_MacIsDecP(pmt)||_MacIsDecP(fv)))	{
		_MacErr(GM_ARGVAL);
		return(GM_NULL);
	}

	if (!(_MacIsDecN(pv)||_MacIsDecN(pmt)||_MacIsDecN(fv)))	 {
		_MacErr(GM_ARGVAL);
		return(GM_NULL);
	}

/* if pmt==0, solve quickly */
/* 0 = pv + fv * (1 + i/100) ^ n ==>
   i = (exp[ ln	(fv / pv) / n] - 1) * 100 */
	if (_MacIsDecZ(pmt)&&fun<=1)  {
		(void) _DivDec80Bit(temp1, fv, pv);
		_MacDChgs(temp1);
		(void) _LnDec80Bit(temp1, temp1);
		(void) ConvLongToDecimal(temp2,	(long) intper);
		(void) _DivDec80Bit(temp1, temp1, temp2);
		(void) _ExpDec80Bit(temp1, temp1);
		(void) _SubDec80Bit(temp1, temp1, &decOne);
		temp1->dc.id-=2;
		(void) _ScaleDec80Bit(intr, temp1, wIntrPrec);
		(void) _Sq5UnsTo4Uns(intr);
		return(intr);
	}

/* find	bounds for intr	*/
/* initial bounds: 0 and 1 (percent) */
	_MacDZero(x1);
	_MacDCopy(x2, (&decOne));
	zcom(y1, fun, intper, fractper,	x1, pv,	pmt, fv, begend);
	zcom(y2, fun, intper, fractper,	x2, pv,	pmt, fv, begend);
	if (_MacIsDecZ(y1))  {
		_MacDCopy(intr,y1);
		return(intr);
	}
	if (_MacIsDecZ(y2))  {
		_MacDCopy(intr,y2);
		return(intr);
	}
	_MacDCopy(y0,y1);
	_MacDCopy(px1, x2);
	_MacDCopy(px2, x1);

/* expanded bounds until zcom has different signs at bounds to find zero */
	i=0;
	while (!_MacIsDecN(y1)^_MacIsDecN(y2))	{
/* change boundary with	smaller	absolute value */
		if ((CompareDecimal(x1,x2)==1)^_MacIsDecN(y1)) {
/* set x2 = 2 *	x2 */
			_MacDCopy(px2,x2);
			_AddDec80Bit(x2,x2,x2);
			zcom(y2, fun, intper, fractper,x2,pv,pmt, fv, begend);
			if (_MacIsDecZ(y2))  {
				_MacDCopy(intr,x2);
				return(intr);
			}
		}

		else  {
/* if x1==0, set x1 = 10%, otherwise x1	= (old x1 + -100%) / 2 */
			_MacDCopy(px1,x1);
			if (_MacIsDecZ(x1))
				(void) ConvLongToDecimal(x1, 10L);
			else  {
				(void) _AddDec80Bit(x1,x1,&decMinusHundred);
				(void) _MulDec80Bit(x1,x1,&decPoint5);
			}
			zcom(y1, fun, intper, fractper,x1,pv,pmt, fv, begend);
			if (_MacIsDecZ(y1))  {
				_MacDCopy(intr,x1);
				return(intr);
			}
		}

/* maximum of 19 expansions */
		i++;
		if (i>19)  {
			_MacErr(GM_ARGVAL);
			return(GM_NULL);
		}

	}

	places=wIntrPrec;
	if (places<GM_MINID)
		places=GM_MINID;
	if (places>GM_MAXID)
		places=GM_MAXID;

	(void) _ScaleDec80Bit(x1,x1,places);
	(void) _ScaleDec80Bit(x2,x2,places);

	if (CompareDecimal(x1,&decMinusHundred)==0)  {
		_MacDCopy(intr,x2);
		return(intr);
	}

/* Restrict range to smaller region with different signs */
	if (_MacIsDecN(y0)^_MacIsDecN(y2))
		_MacDCopy(x1,px2);
	else
		_MacDCopy(x2,px1);


/* Find	root by	bisection (Numerical Recipes in	C, p. 263) */
	wfDone=FALSE;
	while (!wfDone)	 {
/* calculate new position */
		(void) _AddDec80Bit(xn,	x2, x1);
		       _HalveUnsArr(xn->dc.sl, 5);
/*		(void) _SubDec80Bit(temp1, x2, x1);
		(void) _SubDec80Bit(temp2, y2, y1);
		(void) _DivDec80Bit(temp1, temp1, temp2);
		(void) _MulDec80Bit(temp1, temp1, y1);
		(void) _SubDec80Bit(xn,	x1, temp1);	*/

/* calculate corresponding y */
		zcom(yn, fun, intper, fractper,	xn, pv,	pmt, fv, begend);

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

/* 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);
		}

	}

	(void) _Sq5UnsTo4Uns(xn);
	_MacDCopy(intr,xn);
	return(intr);
}


/* calculate function which must be zeroed to solved for intr */
static void	zcom(y,	fun, intper, fractper, intr, pv, pmt, fv, begend)
int	fun, intper, begend;
DEC	*y, *fractper, *intr, *pv, *pmt, *fv;
{
	int	ad;
	DEC	dopi, *opi=&dopi, dopitmn, *opitmn=&dopitmn;
	DEC	dtemp2,	*temp2=&dtemp2,	dopitmna, *opitmna=&dopitmna;
	DEC	dnintr,	*nintr=&dnintr;

	ad=begend;

	_MacDCopy(nintr, intr);
	if (intr->dc.id<=GM_IMAXID-2)
		nintr->dc.id+=2;
	else
		(void) _MulDec80Bit(nintr,intr,&decPoint01);

	(void) _AddDec80Bit(opi, &decOne, nintr);
	(void) _IntPwrDec80Bit(opitmn, opi, -intper);
	if (fun==0)
		(void) _IntPwrDec80Bit(opitmna,	opi, -(intper-ad));

	(void) _MulDec80Bit(y, fv, opitmn);

	if (!_MacIsDecZ(pv))  {
		if (fun<=1)
			_MacDCopy(temp2, pv);
		if (fun==2)  {
			(void) _MulDec80Bit(temp2, nintr, fractper);
			(void) _AddDec80Bit(temp2, temp2, &decOne);
			(void) _MulDec80Bit(temp2, temp2, pv);
		}
		if (fun==3)  {
			(void) PowerDecimal(temp2, opi,	fractper);
			(void) _MulDec80Bit(temp2, temp2, pv);
		}
		(void) _AddDec80Bit(y, y, temp2);
	}

	if (_MacIsDecZ(nintr))	{
		if (fun!=0)
			(void) ConvLongToDecimal(temp2,	(long) intper);
		else
			(void) ConvLongToDecimal(temp2,	(long) (intper-ad));
		(void) _MulDec80Bit(temp2, temp2, pmt);
	}
	else  {
		if (fun!=0)
			(void) _SubDec80Bit(temp2, &decOne, opitmn);
		else
			(void) _SubDec80Bit(temp2, &decOne, opitmna);
		(void) _DivDec80Bit(temp2, temp2, nintr);
		(void) _MulDec80Bit(temp2, temp2, pmt);
		if (begend==GM_BEGIN&&fun!=0)
			(void) _MulDec80Bit(temp2, temp2, opi);
	}
	(void) _AddDec80Bit(y, y, temp2);


}