using System;

namespace Orthogonal.Common
{
	/// <summary>
	/// A class of highly effective and reliable random number generators invented by Pierre L'Ecuyer.
	/// </summary>
	public sealed class RandEcuyer
	{
		private int ecs1 = -1;
		private int ecs2 = -1;

		private long s10 = -1;
		private long s11 = -1;
		private long s12 = -1;
		private long s20 = -1;
		private long s21 = -1;
		private long s22 = -1;

		private const double norm = 1.0842021724855052e-19;
		private const int m21 = 2147483563;
		private const int m22 = 2147483399;
		private const long m1   = 9223372036854769163L;
		private const long m2   = 9223372036854754679L;
		private const long a12  = 1754669720L;
		private const long q12  = 5256471877L;
		private const long r12  = 251304723L;
		private const long a13n = 3182104042L;
		private const long q13  = 2898513661L;
		private const long r13  = 394451401L;
		private const long a21  = 31387477935L;
		private const long q21  = 293855150L;
		private const long r21  = 143639429L;
		private const long a23n = 6199136374L;
		private const long q23  = 1487847900L;
		private const long r23  = 985240079L;

		public RandEcuyer()
		{
			Random rand = new Random();
			// Seed the generators using the standard Rand class.
			ecs1 = (int)(rand.NextDouble() * m21);
			ecs2 = (int)(rand.NextDouble() * m22);
			s10 = (long)(rand.NextDouble() * m1);
			s11 = (long)(rand.NextDouble() * m1);
			s12 = (long)(rand.NextDouble() * m1);
			s20 = (long)(rand.NextDouble() * m2);
			s11 = (long)(rand.NextDouble() * m2);
			s22 = (long)(rand.NextDouble() * m2);
		}

		/// <summary>
		/// Pierre L'Ecuyer's famous 1988 random number generator that combines two
		/// MLGCs to create a period of around 10^18. Numerical Recipies in C improves
		/// this raw algorithm by using a Bayes-Durham shuffle to remove any low order
		/// bit correlations.
		/// </summary>
		/// <returns>A random number in the range 0-1.</returns>
		public double Uniform2()
		{
			MODMULT(53668, 40014, 12211, m21, ref ecs1);
			MODMULT(52774, 40692,  3791, m22, ref ecs2);
			int z = ecs1 - ecs2;
			if (z < 1)
				z += 2147483562;
			return z * 4.65661305956e-10;
		}

		/// <summary>
		/// A helper method that replaces the original C++ macro.
		/// </summary>
		private void MODMULT(int a, int b, int c, int m, ref int s)
		{
			int q = s / a;
			s = b*(s-a*q)-c*q;
			if (s<0)
				s += m;
		}

		/// <summary>
		/// Pierre L'Ecuyer's Combined Multiple Recursive Generator with 2 components of order 3 (see
		/// page 14, <I>Good Parameters and Implementations for Combined Multiple Recursive Random
		/// Number Generators</I>, May 4, 1998). The period length is (m1^3-1)(m2^3-1)/2, which is
		/// around 10^113. The implementation assumes that all integers from -m1 to m1 can
		/// be represented exactly, which is the case with C#'s 64-bit longs. This generator not
		/// only has an astronomically long period, it has excellent structural properties according
		/// to the spectral test up to dimension 30.
		/// </summary>
		/// <returns>A random number in the range 0-1.</returns>
		public double Uniform()
		{
			long h, p12, p13, p21, p23;
			/* Component 1 */
			h = s10 / q13;
			p13 = a13n * (s10 - h * q13) - h * r13;
			h = s11 / q12;
			p12 = a12  * (s11 - h * q12) - h * r12;
			if (p13 < 0)
				p13 += m1;
			if (p12 < 0)
				p12 += m1 - p13;
			else
				p12 -= p13;
			if (p12 < 0)
				p12 += m1;
			s10 = s11;  s11 = s12;  s12 = p12;
			/* Component 2 */
			h = s20 / q23;
			p23 = a23n * (s20 - h * q23) - h * r23;
			h = s22 / q21;
			p21 = a21  * (s22 - h * q21) - h * r21;
			if (p23 < 0)
				p23 += m2;
			if (p21 < 0)
				p21 += m2 - p23;
			else
				p21 -= p23;
			if (p21 < 0)
				p21 += m2;
			s20 = s21;
			s21 = s22;
			s22 = p21;
			/* Combination */
			if (p12 > p21)
				return norm * (p12 - p21);
			else
				return norm * (p12 - p21 + m1);
		}
	}
}