// Filename:    mersenne.v
// Author:      David Ljung Madison <DaveSource.com>
// See License: http://MarginalHacks.com/License
// Version:     1.01
// Description: Verilog implementation of the Mersenne Twister,
//              A great random number generator.
//              With improvements by Paul Mateman (nvidia)
//                (distributions + external seed generator)
//
// Other implementations at: http://www.math.keio.ac.jp/~matumoto/emt.html
//

module mersenne ( );
`define MERSENNE_N        624
`define MERSENNE_M        397
`define MERSENNE_MATRIX_A    32'h9908b0df    /* constant vector a */ 
`define MERSENNE_UPPER_MASK    32'h80000000    /* most significant w-r bits */
`define MERSENNE_LOWER_MASK    32'h7fffffff    /* least significant r bits */

/* Tempering parameters */
`define TEMPERING_MASK_B    32'h9d2c5680
`define TEMPERING_MASK_C    32'hefc60000
`define TEMPERING_SHIFT_U    11
`define TEMPERING_SHIFT_S    7
`define TEMPERING_SHIFT_T    15
`define TEMPERING_SHIFT_L    18

  integer    mt    [0:`MERSENNE_N];   //the array for the state vector
  integer    mti;                     // index in array 
  real    rand_real;
  reg  toggle=1'b1;                   // used to generate rdist_normal in pairs



//============================================================================
//
// task seed(seed_value)
//
//      initialise the generator with a know seed
//
//============================================================================
    task seed;
        input    [31:0]    seed_value;
        begin
            // Setting initial seeds to mt[N] using the generator 
            // Line 25 of Table 1 in [KNUTH 1981, The Art of Computer
            // Programming Vol. 2 (2nd Ed.), pp102]
            mt[0] = seed_value;
            for (mti=1; mti<`MERSENNE_N; mti=mti+1)
                mt[mti] = (69069 * mt[mti-1]);
            $display("Mersenne twister seed=%1d\n", seed_value);                          
        end
    endtask




//============================================================================
//
// task random_seed
//
//      initialise the generator with a random seed
//
//============================================================================
    task random_seed;
        integer rseed_value;
	    integer file_handle;
        integer IOresult;   // dummy
        begin
            $system("bash -c 'echo $RANDOM' > _random_.tmp" );  // hack: first dump result in file
            file_handle = $fopen("_random_.tmp","r");           // open the file
            IOresult =$fscanf(file_handle,"%d",rseed_value);   // read random value (ignore IOresult)
            $fclose(file_handle);
            seed(rseed_value);
        end
    endtask


//============================================================================
//
// task rand(x)
//
//      return a 32 bit random number in x
//
//============================================================================
    task rand;
        output [31:0]    out;
        integer    y;
        integer    kk;
        begin

            if (|mti===1'bx || mti>=`MERSENNE_N) begin
                if (|mti===1'bx) begin
                    $display("You have to initialise a mersenne with a seed. Using default seed=4357");
                    seed(4357);        // Default seed
                end

                for (kk=0; kk<`MERSENNE_N-`MERSENNE_M; kk=kk+1) begin
                    y = (mt[kk]&`MERSENNE_UPPER_MASK)|(mt[kk+1]&`MERSENNE_LOWER_MASK);
                    mt[kk] = mt[kk+`MERSENNE_M] ^ (y >> 1) ^ (y&1 ? `MERSENNE_MATRIX_A : 0);
                end

                for (kk=kk; kk<`MERSENNE_N-1; kk=kk+1) begin
                    y = (mt[kk]&`MERSENNE_UPPER_MASK)|(mt[kk+1]&`MERSENNE_LOWER_MASK);
                    mt[kk] = mt[kk+(`MERSENNE_M-`MERSENNE_N)] ^ (y >> 1) ^ (y&1 ? `MERSENNE_MATRIX_A : 0);
                end
                y = (mt[`MERSENNE_N-1]&`MERSENNE_UPPER_MASK)|(mt[0]&`MERSENNE_LOWER_MASK);
                mt[`MERSENNE_N-1] = mt[`MERSENNE_M-1] ^ (y >> 1) ^ (y&1 ? `MERSENNE_MATRIX_A : 0);

                mti = 0;

            end

            y = mt[mti];
            mti = mti + 1;
            y = y ^ (y >> `TEMPERING_SHIFT_U);
            y = y ^ ((y << `TEMPERING_SHIFT_S) & `TEMPERING_MASK_B);
            y = y ^ ((y << `TEMPERING_SHIFT_T) & `TEMPERING_MASK_C);
            y = y ^ (y >> `TEMPERING_SHIFT_L);

            out = y;
        end
    endtask
    





//============================================================================
//
// task rdist_uniform
//
//      returns nothing, but computes a uniform distributed random number from 0..1, 
//      note that only 32 bits are valid
//
//      access result via
//          mersenne_object_name.rand_real
//      
//
//============================================================================
    task rdist_uniform;
     
       integer rand_int;
    
       begin
            rand(rand_int);
            rand_real = $itor(rand_int) / 32'hffffffff + 0.5;
       end

    endtask //rdist_uniform
 
 
 
 
    
//============================================================================
//
// task rdist_uniform
//
//      returns nothing, but computes a normal distributed random number with mu=0, sigma=1
//      note that only 32 bits are valid
//
//      access result via
//          mersenne_object_name.rand_real
//
//  Box Muller transform
//  see http://en.wikipedia.org/wiki/Box%E2%80%93Muller_transform
//
//============================================================================
    task rdist_normal;
     
       integer rand_int;
       real u1;
       real u2;
       real R;
       real theta;
    
       begin
            if (toggle) begin
                // draw two uniform distributed samples
                rand(rand_int);
                u1 = $itor(rand_int) / 32'hffffffff + 0.5;
                
                rand(rand_int);
                u2 = $itor(rand_int) / 32'hffffffff + 0.5;
                
                // Box-Muller transform
                R = sqrt(-2.0*ln(u1));
								// If `M_TWO_PI is not defined for you, it's (2*pi).  :)
                theta = `M_TWO_PI * u2;
                rand_real = R * cos(theta);
            end
            else begin
                rand_real = R * sin(theta);
            end
            toggle = ~toggle;
       end

    endtask //rdist_normal

endmodule //mersenne