#ifndef CA_UWATERLOO_ALUMNI_DWHARDER_RATIO
#define CA_UWATERLOO_ALUMNI_DWHARDER_RATIO

#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.130-2.

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

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

	double Ratio::randn() {
		while ( true ) {
			double u1 = drand48();
			double u2 = drand48();
			double s = std::sqrt( 8.0/std::exp(1.0) ) * (u2 - 0.5)/u1;

			double ss = s*s;

			if ( ( ss <= 5.0 - 4.0*std::exp(0.25)*u1 ) || (
				( ss < 4.0*std::exp(-1.35)/u1 + 1.4 ) && ( ss <= -4.0/std::log( u1 ) )
			) ) {
				return s;
			}
		}
	}
}

#endif