Program
Program (Link to Latest)[mt_stream_f90.tar.gz] [2011/03/31]

• [mt_stream_f90-1.11.tar.gz][2011/03/31] (Latest)
  A minor bug is corrected. An extra ``space" in the line 628 of program mt_stream.F90 is removed. (Reported by Tridib Sadhu)
• [mt_stream_f90-1.10.tar.gz][2010/07/29]
  Michael Briggs kindly corrected the comments on the jump length in the source codes and improved the interoperability to the NTL/C++ code using the ISO C Binding feature of Fortran 2003. Although the ISO C Binding feature is not the standard of Fortran 90/95, it is already supported by many compilers.
  He also provided test programs to check the consistency between this Fortran module and the original (non-multi-stream) Mersenne Twister C version mt19937ar.c of Takuji Nishimura and Makoto Matsumoto (2002/1/26) [http://www.math.sci.hiroshima-u.ac.jp/~m-mat/MT/MT2002/CODES/mt19937ar.c]. I put them in the `contrib' directory.
  Michael Briggs pointed out a bug in the GF(2)[x] Fortran module of [mt_stream_f90-1.00.tar.gz], and this is fixed.
• [mt_stream_f90-1.00.tar.gz] [2010/03/08] A slow library to compute GF(2)[x], written in Fortran, is included.
  
• [mt_stream_f90-0.95.tar.gz] [2010/02/18] README.html, README.jp.html are included. Comment typos in the source file are corrected.
  
• [mt_stream_f90-0.9.tar.gz] [2010/02/16]

This program module splits single long period random number series from Mersenne Twister (MT) into multiple (almost) independent streams. This enables us the use of parallel Mersenne Twister in a large scale parallel simulation with MPI or OpenMP. This module provides several routines to manipulate multiple stream state and to generate random numbers.

The main part of the program is written in Fortran90/95. A part of the program uses some external libraries written in C/C++ language.
ISO C binding feature of Fortran 2003 is used to link with NTL/GF2X C++ libraries (version 1.10). The list of the compilers which support the feature, to my knowledge, is;

• Intel Fortran version 11.0 and higher.
• gfortran 4.3 and the version higher.

The parameter included in this program is taken from the original code (MT19937) by Makoto Matsumoto and Takuji Nishimura [http://www.math.sci.hiroshima-u.ac.jp/~m-mat/MT/MT2002/emt19937ar.html]. The stream length picked out from the long period is fixed to 2^256.

The method to divide the main stream into several streams is based on the manipulation on the polynomial representation of transition matrix with Cayley-Hamilton theorem for the characteristic polynomial of MT. I only use Horner's method to jump long ahead from the initial state to the next stream initial state. For the details on this method see the following original paper:
[Ref.: H. Haramoto, M. Matsumoto, T. Nishimura, F. Panneton, and P. L'Ecuyer, ``Efficient Jump Ahead for F_2-Linear Random Number Generators'', GERAD Report G-2006-62. INFORMS Journal on Computing, 20, 3 (2008), 385-390.]

In order to obtain the polynomial coefficients to make the MT state to jump by 2^256 steps, we employ the following external libraries in default setting:

• Fast arithmetic in GF(2)[x], [http://wwwmaths.anu.edu.au/~brent/software.html]


• NTL : A Library for doing Number Theory, [http://www.shoup.net/ntl/index.html]

The MT state initialization routine used in this program is the same as that of the original MT19931 with improved initialization by Matusmoto-Nishimura (see also their original web page). I have checked that my generator produces the same first 1000 integer numbers as those of the original MT19931.

To check whether the state jump ahead routine works correctly or not, I have tested the following points. By setting the stream length to 2^15, make a branch form the master stream with the jump ahead routine. I have checked that the branch stream produces the same numbers as those from the master stream, where the master stream simply drops 2^15 random numbers to catch up the branch stream. I also checked this simple test for 1024 multiple streams separated by 2^15 steps each.

Note that I have not completely checked/understood the autocorrelation among multiple disjoint streams picked out from the long period stream.

This MT stream module library can yield the following types of random numbers:

• 32bit integer: Fortran integer(4).
• 53bit double precision: Fortran real(8), in [0,1],(0,1],[0,1).
• 52bit double precision: Fortran real(8), in (0,1).

This program is NO WARRANTY. License follows the New BSD License.

Compile: This program has been developed on Linux with gcc/g++/gfortran and Intel compilers.

1. If C/C++/NTL/GF2X are available, compile and Install the following external libraries:
   
   
   • Fast arithmetic in GF(2)[x], [http://wwwmaths.anu.edu.au/~brent/software.html]
   • NTL : A Library for doing Number Theory, [http://www.shoup.net/ntl/index.html]
   
   If C/C++/NTL/GF2X are not available skip step 1 go next.
   
   
2. Untar program tarball (Latest) [mt_stream_f90.tar.gz].
   Check and Edit the compiler settings in the Makefile.
   Set the include and library paths to the above external libraries in the Makefile.
   If C/C++/NTL/NG2X are not available, edit Makefile to define USE_NTL = no.


3. Type make. You will see some auto consistency checks during compilation.


4. mt_stream.o and mt_stream.mod are the main product of this MT Stream library.
   If you choose to use NTL/GF2X libraries, the following object file is generated.
   • jump_ahead_coeff/get_coeff.o
   If you do not use NTL/GF2X libraries, the following object files are generated.
   • f_jump_ahead_coeff/gf2xe.o
   • f_jump_ahead_coeff/f_get_coeff.o
   You can write programs using the routines from the MT Stream library. To compile your program you need to put 'use mt_stream' statement and to link to mt_stream.o and the above jump ahead object files. Compiler will require mt_stream.mod in the include search path.

The subroutine/function list contained in this MT Stream module.

• Use this module.
  
  use mt_stream


• MT State type
  
  type(mt_state) :: mts


• Set MT19937 parameters to this module.
  
  call set_mt19937


• Initialize MT state type (allocate state vector in the type element). The state is still empty by this initialization.
  
  call new(mts)
  type(mt_state) :: mts
  


• Set seed into MT state type and initialize state.
  
  call init(mts,iseed)
  type(mt_state) :: mts
  integer :: iseed          ! initialization by scalar number
  or
  integer :: iseed(:)        ! initialization by scalar array numbers
  
  The first MT state initialized has the stream id equals to zero.


• Create a new stream (mts_new) from the stream (mts). The (mts_new) state is far from (id*2^256) steps from the (mts) state. The (mts_new) has stream id equals to id.
  
  call create_stream(mts,mts_new,id)
  type(mt_state) :: mts           ! MT stream (id = 0 is required)(input/output)
  type(mt_state) :: mts_new      ! MT stream with id (output). (id*2^256) steps ahead from mts
  integer :: id                    ! stream id for mt_new (id > 0 is required)
  
  mts should be initialized and have state. When the mts state is not aligned at proper step, the mts drops several random numbers to reach aligned step. Then mts_new is created from the aligned mts state. The mts_new points at (id*2^256) steps beyond the mts aligned point. After calling create_stream mts state points aligned point by dropping several random numbers.


• Obtain one 32bit signed integer(integer(4)) random number.
  
  integer :: k
  type(mt_state) :: mts
  k = genrand_int32(mts)
  
  "mts" should be initialized and seeded or created with create_stream. The bit pattern of "k" is the same as that of original C code by Matsumoto-Nishimura. If you want to compare "k" and original one as decimal numbers, convert integer(4) to 64bit integer(8) and add 2_8^32 when k<0. After calling genrand_int32, the state of "mts" proceeds by 32bit.


• Save MT state into a file.
  
  call save(mts,unit)
  type(mt_state) :: mts
  integer :: unit          ! Fortran file unit number
  
  Save the mt state into a file described by the Fortran file unit number(unit). The file should be opened with form='unformatted' before calling this routine. The state mts does not changed after calling save.


• Restore MT state from a file.
  
  call read(mts,unit)
  type(mt_state) :: mts
  integer :: unit          ! Fortran file unit number
  
  Read the mts status information from a file described by the Fortran file unit number(unit) and the status information is restored into mts. The file should be opened with form='unformatted' before calling this routine. The file should contain the proper information which stored by the save routine.


• Free and clear the memory in MT state type.
  
  call delete(mts)
  type(mt_state) :: mts
  


• Obtain one 53/52 bit resolution double precision real uniform random number (Fortran real(8)) from mts.
  
  real(8) :: r
  type(mt_state) :: mts
  r = genrand_double1(mts)          ! r in [0,1]
  r = genrand_double2(mts)          ! r in [0,1)
  r = genrand_double3(mts)          ! r in (0,1)
  r = genrand_double4(mts)          ! r in (0,1]
  
  The random number generator from (0,1) has 52bit resolution. The random number generators from [0,1],[0,1),(0,1] have 53bit resolution.
  After calling these routines, mts state proceeds by 64bits.
  I have not coded the versions with 32bit resolution.

Acknowledgment
I would like to greatly thank to The Mersenne Twister web page by Makoto Matsumoto, at which I could find the original papers and related documents on MT/RNG and very useful links to samples and related programs. Without this page I could not develop this program and understand the philosophy on MT and better PRNG.
I also express my gratitude to the authors of various implementations in C/C++/Fortran etc. on his web page.
I thank Michael Briggs for his kind comments and contributions on the code.
I thank Tridib Sadhu for his bug report.

Comments, Improvements, Contact:
Ken-Ichi Ishikawa
ishikawa[at]theo.phys.sci.hiroshima-u.ac.jp