#ifndef CA_UWATERLOO_ALUMNI_DWHARDER_BOX_MULLER
#define CA_UWATERLOO_ALUMNI_DWHARDER_BOX_MULLER

#include <cassert>
#include <cmath>

// Author:  Douglas Wilhelm Harder
// Copyright (c) 2011 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 Normal_distributions::Box_muller::randn(), this function returns
// a random standard normal value (a mean of zero and standard deviation of one).

namespace Normal_distributions {
	class Box_muller {
		public:
			static double randn();

		private:
			static int n;
			static double second;

			static double rand();
	};

	int Box_muller::n = 1;
	double Box_muller::second = 0.0;

	double Box_muller::randn() {
		++n;
		if ( n & 1 ) {
			return second;
		}

		double u1 = rand();
		double u2 = rand();

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

		second = c*std::sin( pi2*u2 );

		return c*std::cos( pi2*u2 );
	}

	double Box_muller::rand() {
		double x = drand48();
		return (x == 0.0) ? 1.0 : x;
	}
}

#endif