/* normaldist.c */

#include <stdio.h>
#include <math.h>
#include "normaldist.h"

/*
	normal_dist (x) returns cumulative normal distribution
	adapted from erfcc(), p.221, Numerical Recipes in C, 2nd Ed.
*/

double normal_dist (double x)
{
	double t,y,z;
	const double sqrt2 = 1.41421356;
	const double
		a0 = -1.26551223,
		a1 = 1.00002368,
		a2 = 0.37409196,
		a3 = 0.09678418,
		a4 = -0.18628806,
		a5 = 0.27886807,
		a6 = -1.13520398,
		a7 = 1.48851587,
		a8 = -0.82215223,
		a9 = 0.17087277;

	y = fabs (x/sqrt2);
	t = 1.0 / (1.0 + 0.5*y);
	z = a0+t*(a1+t*(a2+t*(a3+t*(a4+t*(a5+t*(a6+t*(a7+t*(a8+t*a9))))))));
	z = 0.5 * t * exp (-y*y + z);
	return x >= 0.0 ? 1.0 - z : z;
}