/*
 * spd
 * May 2006
 *
 * Get a normal[] array of normally distributed random values.
 * Plot them and verify mean and variance (sigma2)
 *
 * Sample usage:

./normal_aprox 0 1000 16000 1  | plot
#######################################################################
m: -0.208904    s: 1004.8
T: 16000 total samples for 16000 used
 * 
 */

#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <errno.h>
#include <sys/time.h>
#include <math.h>

#define	sqr(a)	((a)*(a))
double NormalDist(double m, double s2);

#define DIST_INTERVALS 400
#define SCREEN_W 70
int		dist[DIST_INTERVALS];
#define RANGE 1000.0

int main( int argc, char *argv[] )
{
	int		total=0, i, w=-1, wn;
	double	s2, m, x, sum, med;
	struct timeval  tp;
	int plot;
	double *normal;

	int normal_samples;

	m = (double) atof( argv[1] );
	s2 = (double) atof( argv[2] );
	normal_samples = (int) atoi( argv[3] );
	plot = (int) atoi( argv[4] );

	if(NULL==(normal=malloc(normal_samples*sizeof(double))))
	{
		perror("malloc");
		exit(1);
	}

#if defined _POSIX_VERSION || defined _HPUX_SOURCE
	gettimeofday(&tp, NULL);
#else
	gettimeofday(&tp);
#endif

	/* drand48 return values uniformly distributed [0.0,~1.0) */

	srand48((long)tp.tv_usec);
	for( i = 0; i<normal_samples; total++)
	{
		normal[i++] = NormalDist(m,s2);
	}

	/* normal[] should be a normal (m,s2) distribution */


	if (plot)
	{
		for( i=0; i<normal_samples ; i++ )
		{
			double y;

			y= normal[i] + (double)RANGE * 0.5;
			y *= (double)DIST_INTERVALS/(double)RANGE;
			if (( y >= 0.0 ) && (y<DIST_INTERVALS))
			dist[(int)y]++;
		}
		for( i=0; i<DIST_INTERVALS ; i++ )
			printf("%d %d\n",i, dist[i]);
	}

	for( sum = 0, i=0; i<normal_samples ; i++ )
		sum += normal[i];
	med = sum / (double)normal_samples;

	fprintf( stderr, "m: %g\t", med);

	for( sum = 0.0, i=0; i<normal_samples ; i++ )
		sum += sqr (normal[i] - med);

	fprintf(stderr, "s: %g\n", sum/(double)normal_samples);
	fprintf(stderr, "T: %d total samples for %d used\n", 
		total,
		normal_samples
		);
	
	return 0;
}

/******************************************************************************
 * Devuelve un valor aleatorio con distribucion normal N(m,s2)                *
 * Turon
 *****************************************************************************/
double NormalDist(double m, double s2)
{
   double aux=0.0;
   int i;

   for (i = 1; i < 13; i++)
      aux = aux + drand48();
   return ((aux-6.0) * sqrt(s2) + m);
}