#ifndef CA_UWATERLOO_ALUMNI_DWHARDER_MARSAGLIA_POLAR
#define CA_UWATERLOO_ALUMNI_DWHARDER_MARSAGLIA_POLAR

#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-6.

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

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

		private:
			static int n;
			static double second;

			static double rand();
	};

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

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

		double u1, u2, s;

		do {
			u1 = 2.0*drand48() - 1.0;
			u2 = 2.0*drand48() - 1.0;
			s = u1*u1 + u2*u2;
		} while ( s >= 1 );

		double c = std::sqrt( -2.0*log(s)/s );

		second = u2*c;

		return u1*c;
	}
}

#endif