#ifndef CA_UWATERLOO_ALUMNI_DWHARDER_BOX_MULLER
#define CA_UWATERLOO_ALUMNI_DWHARDER_BOX_MULLER

// Author:  Douglas Wilhelm Harder
// Copyright (c) 2012 by Douglas Wilhelm Harder.  All rights reserved.
// Reference:  Donald E. Knuth, "The Art of Computer Programming, Vol. 2, Seminumerical Algorithms", 3rd Ed., Addison Wesley, Toronto, 1998, p.122.

// On successive calls to randn_bm(), this function returns
// a random standard normal value (a mean of zero and standard deviation of one).

#include <stdlib.h>
#include <math.h>

static double Box_muller_rand() {
	double x = drand48();
	return (x == 0.0) ? 1.0 : x;
}

double randn_bm() {
	static unsigned int n = 1;
	static double second;

	++n;
	if ( n & 1 ) {
		return second;
	}

	double u1 = Box_muller_rand();
	double u2 = Box_muller_rand();

	double c = sqrt( -2.0*log(u1) );
	// atan(1) == pi/4
	// Multiplication by 8 adds 3 to the exponent
	double pi2 = 8.0*atan( 1 );
	second = c*sin( pi2*u2 );

	return c*cos( pi2*u2 );
}

#endif