
(* :Title: RandomMT.m -- a package of Mersenne Twister random number generator*)

(* :Context: ProgrammingInMathematica`RandomMT` *)

(* :Author: Takashi Yoshino *)

(* :Summary:
   The skeleton package is a new type random generator called Mersenne Twister.
 *)

(* :Copyright: © <2006> by <Takashi Yoshino> *)

(* :Package Version: 0.8 *)

(* :Mathematica Version: 3.0 *)

(* :History:
   0.8 for Private use.
*)

(* :Keywords: Mersenne Twister, Random number *)

(* :Sources:
   
*)

(* :Warnings:
   <description of global effects, incompatibilities>
*)

(* :Limitations:
   <special cases not handled, known problems>
*)

(* :Discussion:
   <description of algorithm, information for experts>
*)

(* :Requirements:
*)

(* :Examples:
   <sample input that demonstrates the features of this package>
*)


(* set up the package context, including public imports *)

BeginPackage["RandomMT`"]

(* usage messages for the exported functions and the context itself *)

initGenrand::usage ="initGenrand initialize the seed of the MT scheme"

RandomMT::usage = "RandomMT generates a [0,1) Real Random number"
(* error messages for the exported objects *)

RandomMT::badarg = "No argument is needed."

GaussianRandomMT::usage ="GaussianRandomMT generates a random number whose PDF is N(0,1)"

RandomMT2::usage = "RandomMT2 generates a (0,1] Real Random number"
RandomMT2::badarg = "No argument is needed."

Begin["`Private`"]    (* begin the private context (implementation part) *)

(* protected = Unprotect[ Sin, Cos ] *)

(* definition of auxiliary functions and local (static) variables *)

n = 624;
m = 397;
matrixA = 16^^9908b0df;
upperMask = 16^^80000000;
lowerMask = 16^^7fffffff;
Array[mt, n];
mti = n + 1;
BMFlag=True;
BMRest=0;

(* definition of the exported functions *)
initGenrand[s_Integer] := Module[{},
    mt[0] = BitAnd[s, 16^^ffffffff];
    Table[
      mt[mti] = 
        1812433253*BitXor[mt[mti - 1], IntegerPart[mt[mti - 1]/2^30]] + mti;
      mt[mti] = BitAnd[mt[mti], 16^^ffffffff];,
      {mti, 1, (n - 1)}];
    mti = n; 
    ]

genrandInt32[] := Module[{y, mag01, kk},
    mag01[0] = 0;
    mag01[1] = matrixA;
    If[mti >= n,
      If[mti == (n + 1), initGenrand[5489]];
      Table[
        y = BitOr[BitAnd[mt[kk], upperMask], BitAnd[mt[kk + 1], lowerMask]];
        mt[kk] = 
          BitXor[mt[kk + m], IntegerPart[y/2], mag01[BitAnd[y, 1]]];, {kk, 
          0, (n - m - 1)}
        ];
      Table[
        y = BitOr[BitAnd[mt[kk], upperMask], BitAnd[mt[kk + 1], lowerMask]];
        mt[kk] = 
          BitXor[mt[kk + (m - n)], IntegerPart[y/2], 
            mag01[BitAnd[y, 1]]];, {kk, n - m, n - 2}
        ];
      y = BitOr[BitAnd[mt[n - 1], upperMask], BitAnd[mt[0], lowerMask]];
      mt[n - 1] = BitXor[mt[m - 1], IntegerPart[y/2], mag01[BitAnd[y, 1]]];
      mti = 0;
      ];
    y = mt[mti]; 
    mti++;
    y = BitXor[y, IntegerPart[y/(2^11)]]; 
    y = BitXor[y, BitAnd[y*2^7, 16^^9d2c5680]];
    y = BitXor[y, BitAnd[y*2^15, 16^^efc60000]];
    y = BitXor[y, IntegerPart[y/(2^18)]]
    ]

RandomMT[] := genrandReal2[]
RandomMT2[] := genrandReal4[]

GaussianRandomMT[] := Module[{ran1, ran2},
    If[BMFlag, 
        ran1=RandomMT[]; ran2=RandomMT[];
        BMRest=Sqrt[-2*Log[ran1]]*Cos[2*Pi*ran2];(* nokori wo simau *)
        BMFlag=False; 
        Sqrt[-2*Log[ran1]]*Sin[2*Pi*ran2],
        BMFlag=True;
        BMRest
    ]
]

genrandInt31[] := IntegerPart[genrandInt32[]/2]

genrandReal1[] := genrandInt32[]/4294967295.0

genrandReal2[] := genrandInt32[]/4294967296.0

genrandReal3[] := (genrandInt32[] + 0.5)/4294967296.0

genrandReal4[] := (genrandInt32[] + 1)/4294967296.0

genrandRes53[] := (IntegerPart[genrandInt32[]/(2^5)]*67108864. + IntegerPart[genrandInt32[]/(2^6)])/9007199254740992.0


(* The following comments are needed in future. *)
(* definitions for system functions *)
(* Sin/: Sin[x_]^2 := 1 - Cos[x]^2 *)

(* Protect[ Evaluate[protected] ] *)    (* restore protection of system symbols *)

End[ ]         (* end the private context *)

Protect[initGenrand, genrandInt32, genrandInt31, genrandReal1, genrandReal2, genrandReal3,genrandRes53] 
(* protect exported symbols *)

EndPackage[ ]  (* end the package context *)
