/* DEC	*StandardDeviation(pDst,pSrc,n)
 *
 * ARGUMENT
 *	DEC	*pDst;
 *	DEC	**pSrc;
 *	int n;
 *
 * DESCRIPTION
 *	Sets pDst to the standard deviation of the n DEC elements
 *   of	the array pSrc.
 *   n <= 1 is not allowed.  pDst is unchanged if an error occurs.
 *
 * SIDE	EFFECTS
 *	None.
 *
 * RETURNS
 *	Returns	pointer	to the average if successful, otherwise
 *   returns GM_NULL.
 *
 * POSSIBLE ERROR CODES
 *
 *	GM_NULLPOINTER
 *	GM_ARGVAL
 *	GM_OVERFLOW
 *
 * AUTHOR
 *	Jared Levy
 *   Copyright (C) 1987-1990 Greenleaf Software	Inc.  All rights reserved.
 *
 * MODIFICATIONS
 *
 */

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

DEC	*StandardDeviation(pDst,pSrc,n)
DEC	*pDst;
DEC	**pSrc;
int n;
{
	int i;
	DEC	*p, dsum, *sum=&dsum, ddn, *dn=&ddn, dmean, *mean=&dmean;
	DEC	ddiff, *diff=&ddiff, ddsq, *dsq=&ddsq;

/*  source must	be supplied !! */
	_MacStart(GM_STNDDEV);

	if (!pSrc) {
		_MacErr(GM_NULLPOINTER);
		_MacRet(GM_NULL);
	}
	_MacOutVarD(pDst);

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

	for (i=0; i<n; i++)  {
		_MacInVarD(pSrc[i]);
		}

	p = pSrc[0];
	_MacDCopy(sum, p);

	for (i=1;i<n;i++)  {
		p = pSrc[i];
		(void) _AddDec80Bit(sum, sum, p);
		/* overflow impossible here since sum is 80 bits */
	}

	(void) ConvLongToDecimal(dn, (long) n);
	_DivDec80Bit(mean, sum,	dn);
	(void) _Sq5UnsTo4Uns(mean);
	/* must	succeed	since |average|	<= |largest element| */

	_MacDZero(sum);
	for (i=0;i<n;i++)  {
		p = pSrc[i];
		(void) _SubDec80Bit(diff, p, mean);
		if (_MulDec80Bit(dsq, diff, diff)!=GM_SUCCESS)	{
			_MacErr(GM_OVERFLOW);
			_MacRet(GM_NULL);
		}
		if (_AddDec80Bit(sum, sum, dsq)!=GM_SUCCESS)  {
			_MacErr(GM_OVERFLOW);
			_MacRet(GM_NULL);
		}
	}

	(void)	ConvLongToDecimal(dn, (long) (n-1));
	(void)	_DivDec80Bit(mean, sum,	dn);
		_SqrtDec80Bit(pDst, mean);
	(void)	_Sq5UnsTo4Uns(pDst);
	_MacRet(pDst);
}