/*
 * spd
 * Mar 12  2003
 *
 * plot values for standard normal distribution (0,1)
 * (f, probability density function) or
 * plot values for the cumulative distribution
 * function F.
 */

#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 );
double F ( double x );

int main( int argc, char *argv[] )
{
	int		i, n;
	double	s2, m, x, sum, med;

	if (argv[1] == NULL)
	for( x = -5; x<=5; x+=0.1)
	{
		printf ("%g %g\n",x ,f(x, 0, 1) );
	}
	else
	for( x = -5; x<=5; x+=0.1)
	{
		printf ("%g %g\n",x ,F(x) );
	}

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

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