import java.util.Random;

public final class UtilsRand
{

    static Random standardRand = null;

    // Global state variables for the L'Ecuyer randomizer

    static long s10=-1, s11=-1, s12=-1, s20=-1, s21=-1, s22=-1;

    //----------------------------------------------------------------------------------------------
    /**
     * 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</>, May 4, 1998). The period length is (m1^3-1)(m2^3-1)/2 =
     * <CODE>
     *    307828173409329087991058016384928047770447385542991980602
     *    564803055628462831272662068106119198862352993963568683574
     * </CODE>
     * Which is around 10^113. The implementation assumes that all integers from -m1 to m1 can
     * be represented exactly, which is the case with Java'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.
     */

    public static double uniform()
    {
        final double norm = 1.0842021724855052e-19;
        final long m1   = 9223372036854769163L;
        final long m2   = 9223372036854754679L;
        final long a12  = 1754669720L;
        final long q12  = 5256471877L;
        final long r12  = 251304723L;
        final long a13n = 3182104042L;
        final long q13  = 2898513661L;
        final long r13  = 394451401L;
        final long a21  = 31387477935L;
        final long q21  = 293855150L;
        final long r21  = 143639429L;
        final long a23n = 6199136374L;
        final long q23  = 1487847900L;
        final long r23  = 985240079L;

	    long h, p12, p13, p21, p23;

        if (s10 < 0)
        {
        	// First time through, seed the 3x2 state variables of the generator
        	s10 = (long)(getStandardRand().nextDouble() * m1);
            s11 = (long)(getStandardRand().nextDouble() * m1);
            s12 = (long)(getStandardRand().nextDouble() * m1);
            s20 = (long)(getStandardRand().nextDouble() * m2);
            s21 = (long)(getStandardRand().nextDouble() * m2);
            s22 = (long)(getStandardRand().nextDouble() * m2);
        }
        /* 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);
        }
    //----------------------------------------------------------------------------------------------
    /**
     * Get an instance of the standard randomizer.
     */
    private static Random getStandardRand()
    {
        if (standardRand == null)
            standardRand = new Random();
        return standardRand;
    }

}