/*
 * spd
 * Mar 2003
 *
 * Get a normal[] array of normally distributed random values.
 * Plot them and verify mean and variance (sigma2)
 *
 * We will chose random values from -RANGE/2 to RANGE/2
 * For every value we will discard it according its gaussian probability
 *
 * Sample usage:

./normal 0 1000 16000 1  | plot
#######################################################################
 m: 0.297407     s: 1006.32
 T: 16004241 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 f ( double x, 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++)
	{
		x = (drand48() -0.5) * RANGE;
		/* x is uniformly distributed [-RANGE/2, RANGE/2) */

		if ( drand48() < f(x, m, s2) )
		{
			normal[i++] = x;
			wn=(int)((float)i/normal_samples*SCREEN_W);
			if (wn > w) {fprintf(stderr,"#"); w=wn;}
/*			printf("%d %g\n",i, x);*/
		}
	}
	fprintf(stderr,"\n");

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

/**************************************************************************
 *
 * Normal Distribution
 *
 *                           (x - m)^2
 *              1          - -----------
 * f(x) = -------------- e     2 s2
 *        s (2 M_PI)^0.5
 */
double f ( double x, double m, double s2 )
{
	x -= m;
	x *= x;
	x /= (2.0 * s2);
	x = exp(-x) / sqrt( 2.0 * M_PI * s2 );
	return x;
}

/**************************************************************************
 *
 * Normal Cumulative Distribution
 * Hull's book p.243 
 *
 *    / x
 *    |       f(t)     dt
 *   / -inf       0,1
 */
double F ( double x )
{
	double t,Fx;


	if ( x == 0.0 )
		return 0.5;

	t = 1.0 / (1 + 0.2316419 * fabs(x));
	Fx= f(x, 0.0, 1.0) *
		(
		 (
		  (
		   (1.330274429 * t - 1.821255978) * t
		   + 1.78147937
		  ) *t
		  - 0.356563782
		 ) *t
		 + 0.31938153
		) *t ;
		
	if ( x >= 0.0 )
		return( 1 - Fx );
	else
		return Fx;
}

/**************************************************************************
 *
 * Uniform distribution
 */
double u ( double x, double min, double max )
{
	if (( x <= min ) || ( x >= max )) return 0;
	else
		return ( 1 / ( max - min ) );
}