/* normal.c */

#include <stdio.h>
#include <math.h>
#include "random.h"
#include "normal.h"

/*
	normal_rand() returns normally distributed random doubles
	uses polar method for normal deviates, "Algorithm P" on p.117 of Knuth v.2
*/

double normal_rand (void)
{
	static int flag = 0;
	static double z, a = 2147483648.0;
	double v1, v2, s;
	
	if (flag) {
		flag = 0;
		return z;
	}
	flag = 1;

	do {
		v1 = urand()/a - 1;
		v2 = urand()/a - 1;
	}
	while ((s = v1*v1 + v2*v2) > 1.0);
	s = sqrt (-2.0 * log(s) / s);
	z = v1 * s;
	return v2 * s;
}