campo-sirio/gfm/linest.c

/* DEC	*LinearEstimate(pCorr, pLinA, pLinB, pSrcX, pSrcY, wSize);
 *
 * ARGUMENT
 *	DEC	*pCorr,	*pLinA,	*pLinB;
 *	DEC	**pSrcX, **pSrcY;
 *	int	wSize;
 *
 * DESCRIPTION
 *	Determines the least squares line y = A	+ B * x	and the
 *   linear correlation	coefficient of a given set of (X, Y) values.
 *   The wSize (X, Y) values, where wSize > 1, are stored in arrays
 *   pSrcX and pSrcY.  The correlation coefficient is put in pCorr.
 *   The value A is put	in the global pGMStatA,	along with pLinA,
 *   while B is	put in the global pGMStatB and pLinB.  Storing A and
 *   B in globals allow	their later usage in PredictX and PredictY.
 *   The correlation coefficient is returned.
 *
 * SIDE	EFFECTS
 *	None.
 *
 * RETURNS
 *	Returns	pointer	to the correlation coefficient if successful,
 *   otherwise returns GM_NULL.
 *
 * POSSIBLE ERROR CODES
 *
 *	GM_NULLPOINTER
 *	GM_ARGVAL
 *	GM_OVERFLOW
 *	GM_UNDERFLOW
 *
 * AUTHOR
 *	Jared Levy
 *   Copyright (C) 1987-1990 Greenleaf Software	Inc.  All rights reserved.
 *
 * MODIFICATIONS
 *
 */

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

DEC	*LinearEstimate(pCorr, pLinA, pLinB, pSrcX, pSrcY, wSize)
DEC	*pCorr,	*pLinA,	*pLinB;
DEC	**pSrcX, **pSrcY;
int	wSize;
{
	int i;
	DEC	dSumX,	*pSumX=&dSumX,	 dSumY,	 *pSumY=&dSumY;
	DEC	dSumXX,	*pSumXX=&dSumXX, dSumYY, *pSumYY=&dSumYY;
	DEC	dSumXY,	*pSumXY=&dSumXY, dSize,	 *pSize=&dSize;
	DEC	dProd,	*pProd=&dProd,	 dProd2, *pProd2=&dProd2;
	DEC	dVarX,	*pVarX=&dVarX,	 dVarY,	 *pVarY=&dVarY;
	DEC	dNum,	*pNum=&dNum,	 dDen,	 *pDen=&dDen;
	_MacStart(GM_LINEST);

	if (!pSrcX||!pSrcY||!pCorr||!pLinA||!pLinB)  {
		_MacErr(GM_NULLPOINTER);
		_MacRet(GM_NULL);
	}

	if (wSize<=1)  {
		_MacErr(GM_ARGVAL);
		_MacRet(GM_NULL);
	}

	for (i=0; i<wSize; i++)	 {
		_MacInVarD(pSrcX[i]);
		_MacInVarD(pSrcY[i]);
	}

/* |pSrcX[i]| && |pSrcY[i]| must be < decMaxLinEst (1e6) to avoid overflow */
/* values less than 1e-6 may also lead to incorrect results */
	for (i=0; i<wSize; i++)
		if ((CompareDecimal(pSrcX[i],&decMaxLinEst)==1)	||
			(CompareDecimal(pSrcY[i],&decMaxLinEst)==1) ||
			(CompareDecimal(pSrcX[i],&decMinLinEst)==-1) ||
			(CompareDecimal(pSrcY[i],&decMinLinEst)==-1))  {
			_MacErr(GM_ARGVAL);
			_MacRet(GM_NULL);
		}

/* calculate the sum of	X, Y, X^2, Y^2,	and XY */
	_MacDZero(pSumX);
	_MacDZero(pSumY);
	_MacDZero(pSumXX);
	_MacDZero(pSumYY);
	_MacDZero(pSumXY);

	for (i=0;i<wSize;i++)  {
		(void) _AddDec80Bit(pSumX, pSumX, pSrcX[i]);
		(void) _AddDec80Bit(pSumY, pSumY, pSrcY[i]);

		(void) _MulDec80Bit(pProd, pSrcX[i], pSrcX[i]);
		(void) _AddDec80Bit(pSumXX, pSumXX, pProd);

		(void) _MulDec80Bit(pProd, pSrcY[i], pSrcY[i]);
		(void) _AddDec80Bit(pSumYY, pSumYY, pProd);

		(void) _MulDec80Bit(pProd, pSrcX[i], pSrcY[i]);
		(void) _AddDec80Bit(pSumXY, pSumXY, pProd);
		}

/*  Compute n*sum(X^2) - sum(X)*sum(X)	*/
	(void)	ConvLongToDecimal(pSize, (long)	wSize);
	(void) _MulDec80Bit(pProd, pSumXX, pSize);
	(void) _MulDec80Bit(pProd2, pSumX, pSumX);
	(void) _SubDec80Bit(pVarX, pProd, pProd2);
	if (_MacIsDecZ(pVarX))	{
		_MacErr(GM_UNDERFLOW);
		_MacRet(GM_NULL);
	}

/*  Compute A  */
	(void) _MulDec80Bit(pProd, pSumXX, pSumY);
	(void) _MulDec80Bit(pProd2, pSumX, pSumXY);
	(void) _SubDec80Bit(pNum, pProd, pProd2);
	if (!DivideDecimal(pLinA,pNum,pVarX))
		_MacRet(GM_NULL);
	_MacDCopy(pGMStatA, pLinA);

/*  Compute B  */
	(void) _MulDec80Bit(pProd, pSumXY, pSize);
	(void) _MulDec80Bit(pProd2, pSumX, pSumY);
	(void) _SubDec80Bit(pNum, pProd, pProd2);
	if (!DivideDecimal(pLinB,pNum,pVarX))
		_MacRet(GM_NULL);
	_MacDCopy(pGMStatB, pLinB);

/*  Compute correlation	coefficient  */
	(void) _MulDec80Bit(pProd, pSumYY, pSize);
	(void) _MulDec80Bit(pProd2, pSumY, pSumY);
	(void) _SubDec80Bit(pVarY, pProd, pProd2);
	(void) _MulDec80Bit(pProd, pVarX, pVarY);
	       _SqrtDec80Bit(pDen, pProd);
	if (_MacIsDecZ(pDen))  {
		_MacErr(GM_UNDERFLOW);
		_MacRet(GM_NULL);
	}
	if (!DivideDecimal(pCorr, pNum,	pDen))
		_MacRet(GM_NULL);

	_MacRet(pCorr);
}