#include "rand.h"

namespace uty {

namespace {
	long s1 = -1;
	long s2 = -1;
	double norm = 2.3283163396834613e-10;
	double m1   = 4294949027.0;
	double m2   = 4294934327.0;
	double a12  =    1154721.0;
	double a14  =    1739991.0;
	double a15n =    1108499.0;
	double a21  =    1776413.0;
	double a23  =     865203.0;
	double a25n =    1641052.0;
	double s10  = -1;
	double s11  = -1;
	double s12  = -1;
	double s13  = -1;
	double s14  = -1;
	double s20  = -1;
	double s21  = -1;
	double s22  = -1;
	double s23  = -1;
	double s24  = -1;
}

//******************************************************************************

#define MODMULT(a,b,c,m,s) q=s/a; s=b*(s-a*q)-c*q; if (s<0) s+=m;

double ecuyer0()
{
	long q,z;
	if (s1 < 0)
	{
		s1 = long(2147483563L * SimpleUniform());
		s2 = long(2147483399L * SimpleUniform());
	}
	MODMULT(53668,40014,12211,2147483563L,s1);
	MODMULT(52774,40692,3791,2147483399L,s2);
	z = s1 - s2;
	if (z < 1)
		z += 2147483562;
	return z * 4.65661305956e-10;
}

//******************************************************************************

// CMRG with 2 components of order 5
//  period = 1067929815034781460332887887493625356722083142367607302091170291084074196106172588975581863502918
//         = 2^318.9999139445751006757668477712126100260
// ~ 2^319 = 1067993517960455041197510853084776057301352261178326384973520803911109862890320275011481043468288
// ~ 10^96.02854271154581696754941894339674477547

double ecuyer2()
{
    long   k;
    double d,p1, p2;
	if (s10 < 0)
	{
		s10 = floor(m1 * SimpleUniform());
		s11 = floor(m1 * SimpleUniform());
		s12 = floor(m1 * SimpleUniform());
		s13 = floor(m1 * SimpleUniform());
		s14 = floor(m1 * SimpleUniform());
		s20 = floor(m2 * SimpleUniform());
		s21 = floor(m2 * SimpleUniform());
		s22 = floor(m2 * SimpleUniform());
		s23 = floor(m2 * SimpleUniform());
		s24 = floor(m2 * SimpleUniform());
	}
    /* Component 1 */
    p1 = a12 * s13 - a15n * s10;
    if (p1 > 0.0)
		p1 -= a14 * m1;
    p1 += a14 * s11;
	k = long(p1 / m1);
	p1 -= k * m1;
    if (p1 < 0.0)
		p1 += m1;
    s10 = s11;
	s11 = s12;
	s12 = s13;
	s13 = s14;
	s14 = p1;
    /* Component 2 */
    p2 = a21 * s24 - a25n * s20;
    if (p2 > 0.0)
		p2 -= a23 * m2;
    p2 += a23 * s22;
	k  = long(p2 / m2);
	p2 -= k * m2;
    if (p2 < 0.0)
		p2 += m2;
    s20 = s21;
	s21 = s22;
	s22 = s23;
	s23 = s24;
	s24 = p2;
    /* Combination */
    if (p1 <= p2)
		d = ((p1 - p2 + m1) * norm);
    else
		d = ((p1 - p2) * norm);
	return d;
}

//******************************************************************************

double SimpleUniform()
{
	return (double)rand() / RAND_MAX;
}

}	// namespace