package edu.cornell.lassp.houle.RngPack;
import java.util.*;
import java.io.Serializable;

//
// RngPack 1.1a by Paul Houle
// http://www.honeylocust.com/RngPack/
//

/**
 *
 * <A HREF="http://www.math.keio.ac.jp/matumoto/emt.htm">Mersenne Twister</A> --
 * advanced psuedorandom generator with a period of 2<SUP>19937</SUP>-1
 * <P>
 *
 * <P>
 * <A HREF="/RngPack/src/edu/cornell/lassp/houle/RngPack/RanMT.java">
 * Source code </A> is available. 
 * <P>
 * This class is derived from Sean Luke's 
 * <A HREF="http://www.cs.umd.edu/users/seanl/gp/">implementation</A>
 *
 *
 *
 *
 * <h3>License</h3>
 *
 * Copyright (c) 2003 by Paul Houle. <br>
 * Derived from a work copyright (c) 2003 by Sean Luke. <br>
 * Portions copyright (c) 1993 by Michael Lecuyer. <br>
 * All rights reserved. <br>
 *
 * <p>Redistribution and use in source and binary forms, with or without 
 * modification, are permitted provided that the following conditions are met:
 * <ul>
 * <li> Redistributions of source code must retain the above copyright notice, 
 * this list of conditions and the following disclaimer.
 * <li> Redistributions in binary form must reproduce the above copyright notice, 
 * this list of conditions and the following disclaimer in the documentation 
 * and/or other materials provided with the distribution.
 * <li> Neither the name of the copyright owners, their employers, nor the 
 * names of its contributors may be used to endorse or promote products 
 * derived from this software without specific prior written permission.
 * </ul>
 * <p>THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" 
 * AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE 
 * IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE 
 * DISCLAIMED.  IN NO EVENT SHALL THE COPYRIGHT OWNERS OR CONTRIBUTORS BE 
 * LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR 
 * CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF 
 * SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS 
 * INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN 
 * CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) 
 * ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE 
 * POSSIBILITY OF SUCH DAMAGE.
 *
 *
 *
 * @author <A HREF="http://www.honeylocust.com/"> Paul Houle </A> (E-mail: <A HREF="mailto:paul@honeylocust.com">paul@honeylocust.com</A>)
 * @version 1.1a
 */


public class RanMT extends RandomSeedable implements Serializable {

    private static final int N = 624;
    private static final int M = 397;
    private static final int MATRIX_A = 0x9908b0df;   //    private static final * constant vector a
    private static final int UPPER_MASK = 0x80000000; // most significant w-r bits
    private static final int LOWER_MASK = 0x7fffffff; // least significant r bits
    
    
    // Tempering parameters
    private static final int TEMPERING_MASK_B = 0x9d2c5680;
    private static final int TEMPERING_MASK_C = 0xefc60000;
    
    private int mt[]; // the array for the state vector
    private int mti; // mti==N+1 means mt[N] is not initialized
    private int mag01[];

    public RanMT() {
	this.setSeed(4357);
    }
    
    /*
     *
     * Note that this uses only the LSB 32 bits from the long
     *
     */

    public RanMT(long seed) {
	this.setSeed(seed);
    }

    public RanMT(Date d) {
	this.setSeed(d.getTime());
    }

    /**
     *
     * If a 32 bit seed isn't enough for you,  you can pass an array of
     * 624 integers.  Any array of integers is fine so long as they 
     * aren't all zero.
     *
     */

    public RanMT(int array[]) {
	this.setSeed(array);
    }

    private void setSeed(final long seed)
    {
	mt = new int[N];
	
	mag01 = new int[2];
	mag01[0] = 0x0;
	mag01[1] = MATRIX_A;
	
        mt[0]= (int)(seed & 0xfffffff);
        for (mti=1; mti<N; mti++) 
            {
		mt[mti] = 
		    (1812433253 * (mt[mti-1] ^ (mt[mti-1] >>> 30)) + mti); 
		/* See Knuth TAOCP Vol2. 3rd Ed. P.106 for multiplier. */
		/* In the previous versions, MSBs of the seed affect   */
		/* only MSBs of the array mt[].                        */
		/* 2002/01/09 modified by Makoto Matsumoto             */
		mt[mti] &= 0xffffffff;
		/* for >32 bit machines */
            }
    }

    /**
     * An alternative, more complete, method of seeding the
     * pseudo random number generator.  array must be an 
     * array of 624 ints, and they can be any value as long as
     * they're not *all* zero.
     */

    private void setSeed(final int[] array) {
        int i, j, k;
        setSeed(19650218);
        i=1; j=0;
        k = (N>array.length ? N : array.length);
        for (; k!=0; k--) {
            mt[i] = (mt[i] ^ ((mt[i-1] ^ (mt[i-1] >>> 30)) * 1664525)) + array[j] + j; /* non linear */
            mt[i] &= 0xffffffff; /* for WORDSIZE > 32 machines */
            i++;
            j++;
            if (i>=N) { mt[0] = mt[N-1]; i=1; }
            if (j>=array.length) j=0;
	}
        for (k=N-1; k!=0; k--) {
            mt[i] = (mt[i] ^ ((mt[i-1] ^ (mt[i-1] >>> 30)) * 1566083941)) - i; /* non linear */
            mt[i] &= 0xffffffff; /* for WORDSIZE > 32 machines */
            i++;
            if (i>=N) {
                mt[0] = mt[N-1]; i=1; 
	    }
	}
        mt[0] = 0x80000000; /* MSB is 1; assuring non-zero initial array */ 
    }




    public final double raw() {
	int y;
	int z;
	
	if (mti >= N) {   // generate N words at one time
	    
	    int kk;
	    
	    for (kk = 0; kk < N - M; kk++) {
		y = (mt[kk] & UPPER_MASK) | (mt[kk+1] & LOWER_MASK);
		mt[kk] = mt[kk+M] ^ (y >>> 1) ^ mag01[y & 0x1];
	    }
	    for (; kk < N-1; kk++) {
		y = (mt[kk] & UPPER_MASK) | (mt[kk+1] & LOWER_MASK);
		mt[kk] = mt[kk+(M-N)] ^ (y >>> 1) ^ mag01[y & 0x1];
	    }
	    y = (mt[N-1] & UPPER_MASK) | (mt[0] & LOWER_MASK);
	    mt[N-1] = mt[M-1] ^ (y >>> 1) ^ mag01[y & 0x1];
	    
	    mti = 0;
	}
  
	y = mt[mti++];
	y ^= y >>> 11;                          // TEMPERING_SHIFT_U(y)
	y ^= (y << 7) & TEMPERING_MASK_B;       // TEMPERING_SHIFT_S(y)
	y ^= (y << 15) & TEMPERING_MASK_C;      // TEMPERING_SHIFT_T(y)
	y ^= (y >>> 18);                        // TEMPERING_SHIFT_L(y)
	
	if (mti >= N) {  // generate N words at one time
	    int kk;
	    
	    for (kk = 0; kk < N - M; kk++) {
		z = (mt[kk] & UPPER_MASK) | (mt[kk+1] & LOWER_MASK);
		mt[kk] = mt[kk+M] ^ (z >>> 1) ^ mag01[z & 0x1];
	    }
	    for (; kk < N-1; kk++) {
		z = (mt[kk] & UPPER_MASK) | (mt[kk+1] & LOWER_MASK);
		mt[kk] = mt[kk+(M-N)] ^ (z >>> 1) ^ mag01[z & 0x1];
	    }
	    z = (mt[N-1] & UPPER_MASK) | (mt[0] & LOWER_MASK);
	    mt[N-1] = mt[M-1] ^ (z >>> 1) ^ mag01[z & 0x1];
	    
	    mti = 0;
	}
	
	z = mt[mti++];
	z ^= z >>> 11;                          // TEMPERING_SHIFT_U(z)
	z ^= (z << 7) & TEMPERING_MASK_B;       // TEMPERING_SHIFT_S(z)
	z ^= (z << 15) & TEMPERING_MASK_C;      // TEMPERING_SHIFT_T(z)
	z ^= (z >>> 18);                        // TEMPERING_SHIFT_L(z)
	
	/* derived from nextDouble documentation in jdk 1.2 docs, see top */
	return ((((long)(y >>> 6)) << 27) + (z >>> 5)) / (double)(1L << 53);
    }
}

