using System;

namespace Orthogonal.Common
{
	#region START_ECUYER

	/// <summary>
	/// A class of highly effective and reliable random number generators invented by Pierre L'Ecuyer.
	/// </summary>
	public sealed class RandEcuyer
	{
		public RandEcuyer()
		{
			iv = new long[32];
			AM = 1.0 / 2147483563.0;
			NDIV = 1.0 + 2147483562.0 / 32.0;
			SeedRandom();
		}

		/// <summary>
		/// Uses the standard .NET Random class to reseed all generators with random values.
		/// </summary>
		public void SeedRandom()
		{
			Random rand = new Random();
			e1m1 = (int)(rand.NextDouble() * E1Max1);
			e1m2 = (int)(rand.NextDouble() * E1Max2);

			e2m1 = (int)(rand.NextDouble() * IM1);
			e2m2 = (int)(rand.NextDouble() * IM2);

			s10 = (long)(rand.NextDouble() * E3Max1);
			s11 = (long)(rand.NextDouble() * E3Max1);
			s12 = (long)(rand.NextDouble() * E3Max1);
			s20 = (long)(rand.NextDouble() * E3Max2);
			s21 = (long)(rand.NextDouble() * E3Max2);
			s22 = (long)(rand.NextDouble() * E3Max2);
		}

		#region Ecuyer1

		private int e1m1;
		private int e1m2;
		private const int E1Max1 = 2147483563;
		private const int E1Max2 = 2147483399;

		/// <summary>
		/// Pierre L'Ecuyer's famous 1988 random number generator that combines two
		/// MLGCs to create a period of around 10^18. This is the original algorithm
		/// without any enhancements or shuffles.
		/// </summary>
		/// <returns>A random number in the range 0.0 to 1.0 exclusive of endpoints.</returns>
		public double Ecuyer1()
		{
			MODMULT(53668, 40014, 12211, E1Max1, ref e1m1);
			MODMULT(52774, 40692, 3791, E1Max2, ref e1m2);
			int z = e1m1 - e1m2;
			if (z < 1)
				z += (E1Max1 - 1); //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;
		}

		#endregion

		#region Ecuyer2

		private const long IM1 = 2147483563;
		private const long IM2 = 2147483399;
		private long[] iv;
		private long iy;
		private long e2m1;
		private long e2m2;
		private double AM;
		private double NDIV;
		private long s10;
		private long s11;
		private long s12;
		private long s20;
		private long s21;
		private long s22;

		/// <summary>
		/// Pierre L'ecuyer's 1988 random number generator that combines two MLCGs, with the
		/// additional safeguard of a Bayes-Durham shuffle to reduce correlation.
		/// </summary>
		/// <returns>A random number in the range 0.0 to 1.0 exclusive of endpoints.</returns>
		public double Ecuyer2()
		{
			const int NTAB = 32;
			const long IQ1 = 53668;
			const long IQ2 = 52774;
			const long IR1 = 12211;
			const long IR2 = 3791;
			const long IA1 = 40014;
			const long IA2 = 40692;
			const long IMM1 = 2147483562;
			const double RNMX = 1.0 - 1.2e-7;
			int j;
			long k;
			double temp;
			if (e2m1 <= 0)
			{
				if (-(e2m1) < 1) e2m1 = 1;
				else e2m1 = -(e2m1);
				e2m2 = (e2m1);
				for (j = NTAB + 7; j >= 0; j--)
				{
					k = (e2m1) / IQ1;
					e2m1 = IA1 * (e2m1 - k * IQ1) - k * IR1;
					if (e2m1 < 0) e2m1 += IM1;
					if (j < NTAB) iv[j] = e2m1;
				}
				iy = iv[0];
			}
			k = (e2m1) / IQ1;
			e2m1 = IA1 * (e2m1 - k * IQ1) - k * IR1;
			if (e2m1 < 0) e2m1 += IM1;
			k = e2m2 / IQ2;
			e2m2 = IA2 * (e2m2 - k * IQ2) - k * IR2;
			if (e2m2 < 0) e2m2 += IM2;
			j = Convert.ToInt32(iy / NDIV) % NTAB;
			iy = iv[j] - e2m2;
			iv[j] = e2m1;
			if (iy < 1) iy += IMM1;
			if ((temp = AM * iy) > RNMX) return RNMX;
			else return temp;
		}

		#endregion

		#region Encuyer3

		private const long E3Max1 = 9223372036854769163L;
		private const long E3Max2 = 9223372036854754679L;

		/// <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 Ecuyer3()
		{
			const double norm = 1.0842021724855052e-19;
			const long a13n = 3182104042L;
			const long a23n = 6199136374L;
			const long r23 = 985240079L;
			const long q23 = 1487847900L;
			const long r21 = 143639429L;
			const long q21 = 293855150L;
			const long a21 = 31387477935L;
			const long r13 = 394451401L;
			const long q13 = 2898513661L;
			const long r12 = 251304723L;
			const long q12 = 5256471877L;
			const long a12 = 1754669720L;
			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 += E3Max1;
			if (p12 < 0)
				p12 += E3Max1 - p13;
			else
				p12 -= p13;
			if (p12 < 0)
				p12 += E3Max1;
			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 += E3Max2;
			if (p21 < 0)
				p21 += E3Max2 - p23;
			else
				p21 -= p23;
			if (p21 < 0)
				p21 += E3Max2;
			s20 = s21;
			s21 = s22;
			s22 = p21;
			/* Combination */
			if (p12 > p21)
				return norm * (p12 - p21);
			else
				return norm * (p12 - p21 + E3Max1);
		}

		#endregion
	}

	#endregion END_ECUYER
}