MersenneTwister.F90

Go to the documentation of this file.

! Fortran-95 implementation of the Mersenne Twister 19937, following 
!   the C implementation described below (code mt19937ar-cok.c, dated 2002/2/10), 
!   adapted cosmetically by making the names more general.  
! Users must declare one or more variables of type randomNumberSequence in the calling 
!   procedure which are then initialized using a required seed. If the 
!   variable is not initialized the random numbers will all be 0. 
! For example: 
! program testRandoms 
!   use RandomNumbers
!   type(randomNumberSequence) :: randomNumbers
!   integer                    :: i
!   
!   randomNumbers = new_RandomNumberSequence(seed = 100)
!   do i = 1, 10
!     print ('(f12.10, 2x)'), getRandomReal(randomNumbers)
!   end do
! end program testRandoms
! 
! Fortran-95 implementation by 
!   Robert Pincus
!   NOAA-CIRES Climate Diagnostics Center
!   Boulder, CO 80305 
!   email: Robert.Pincus@colorado.edu
!
! This documentation in the original C program reads:
! -------------------------------------------------------------
!    A C-program for MT19937, with initialization improved 2002/2/10.
!    Coded by Takuji Nishimura and Makoto Matsumoto.
!    This is a faster version by taking Shawn Cokus's optimization,
!    Matthe Bellew's simplification, Isaku Wada's real version.
! 
!    Before using, initialize the state by using init_genrand(seed) 
!    or init_by_array(init_key, key_length).
! 
!    Copyright (C) 1997 - 2002, Makoto Matsumoto and Takuji Nishimura,
!    All rights reserved.                          
! 
!    Redistribution and use in source and binary forms, with or without
!    modification, are permitted provided that the following conditions
!    are met:
! 
!      1. Redistributions of source code must retain the above copyright
!         notice, this list of conditions and the following disclaimer.
! 
!      2. 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.
! 
!      3. The names of its contributors may not be used to endorse or promote 
!         products derived from this software without specific prior written 
!         permission.
! 
!    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 OWNER 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.
! 
! 
!    Any feedback is very welcome.
!    http://www.math.keio.ac.jp/matumoto/emt.html
!    email: matumoto@math.keio.ac.jp
! -------------------------------------------------------------

module mersennetwister_mod
! -------------------------------------------------------------
  implicit none
  private
  
  ! Algorithm parameters
  ! -------
  ! Period parameters
  integer, parameter :: blocksize = 624,         &
                        m         = 397,         &
                        matrix_a  = -1727483681, & ! constant vector a         (0x9908b0dfUL)
                        umask     = -2147483648_8, & ! most significant w-r bits (0x80000000UL)
                        lmask     =  2147483647    ! least significant r bits  (0x7fffffffUL)
  ! Tempering parameters
  integer, parameter :: tmaskb= -1658038656, & ! (0x9d2c5680UL)
                        tmaskc= -272236544     ! (0xefc60000UL)
  ! -------

  ! The type containing the state variable  
  type randomnumbersequence
    integer                            :: currentelement ! = blockSize
    integer, dimension(0:blockSize -1) :: state ! = 0
  end type randomnumbersequence

  interface new_randomnumbersequence
    module procedure initialize_scalar, initialize_vector
  end interface new_randomnumbersequence 

  public :: randomnumbersequence
  public :: new_randomnumbersequence, finalize_randomnumbersequence, &
            getrandomint, getrandompositiveint, getrandomreal
! -------------------------------------------------------------
contains
  ! -------------------------------------------------------------
  ! Private functions
  ! ---------------------------
  function mixbits(u, v)
    integer, intent( in) :: u, v
    integer              :: mixbits
    
    mixbits = ior(iand(u, umask), iand(v, lmask))
  end function mixbits
  ! ---------------------------
  function twist(u, v)
    integer, intent( in) :: u, v
    integer              :: twist

    ! Local variable
    integer, parameter, dimension(0:1) :: t_matrix = (/ 0, matrix_a /)
    
    twist = ieor(ishft(mixbits(u, v), -1), t_matrix(iand(v, 1)))
    twist = ieor(ishft(mixbits(u, v), -1), t_matrix(iand(v, 1)))
  end function twist
  ! ---------------------------
  subroutine nextstate(twister)
    type(randomNumberSequence), intent(inout) :: twister
    
    ! Local variables
    integer :: k
    
    do k = 0, blocksize - m - 1
      twister%state(k) = ieor(twister%state(k + m), &
                              twist(twister%state(k), twister%state(k + 1)))
    end do 
    do k = blocksize - m, blocksize - 2
      twister%state(k) = ieor(twister%state(k + m - blocksize), &
                              twist(twister%state(k), twister%state(k + 1)))
    end do 
    twister%state(blocksize - 1) = ieor(twister%state(m - 1), &
                                        twist(twister%state(blocksize - 1), twister%state(0)))
    twister%currentElement = 0

  end subroutine nextstate
  ! ---------------------------
  elemental function temper(y)
    integer, intent(in) :: y
    integer             :: temper
    
    integer :: x
    
    ! Tempering
    x      = ieor(y, ishft(y, -11))
    x      = ieor(x, iand(ishft(x,  7), tmaskb))
    x      = ieor(x, iand(ishft(x, 15), tmaskc))
    temper = ieor(x, ishft(x, -18))
  end function temper
  ! -------------------------------------------------------------
  ! Public (but hidden) functions
  ! --------------------
  function initialize_scalar(seed) result(twister)
    integer,       intent(in   ) :: seed
    type(randomnumbersequence)                :: twister 
    
    integer :: i
    ! See Knuth TAOCP Vol2. 3rd Ed. P.106 for multiplier. In the previous versions, 
    !   MSBs of the seed affect only MSBs of the array state[].                       
    !   2002/01/09 modified by Makoto Matsumoto            
    
    twister%state(0) = iand(seed, -1)
    do i = 1,  blocksize - 1 ! ubound(twister%state)
       twister%state(i) = 1812433253 * ieor(twister%state(i-1), &
                                            ishft(twister%state(i-1), -30)) + i
       twister%state(i) = iand(twister%state(i), -1) ! for >32 bit machines
    end do
    twister%currentElement = blocksize
  end function initialize_scalar
  ! -------------------------------------------------------------
  function initialize_vector(seed) result(twister)
    integer, dimension(0:), intent(in) :: seed
    type(randomnumbersequence)                      :: twister 
    
    integer :: i, j, k, nfirstloop, nwraps
    
    nwraps  = 0
    twister = initialize_scalar(19650218)
    
    nfirstloop = max(blocksize, size(seed))
    do k = 1, nfirstloop
       i = mod(k + nwraps, blocksize)
       j = mod(k - 1,      size(seed))
       if(i == 0) then
         twister%state(i) = twister%state(blocksize - 1)
         twister%state(1) = ieor(twister%state(1),                                 &
                                 ieor(twister%state(1-1),                          & 
                                      ishft(twister%state(1-1), -30)) * 1664525) + & 
                            seed(j) + j ! Non-linear
         twister%state(i) = iand(twister%state(i), -1) ! for >32 bit machines
         nwraps = nwraps + 1
       else
         twister%state(i) = ieor(twister%state(i),                                 &
                                 ieor(twister%state(i-1),                          & 
                                      ishft(twister%state(i-1), -30)) * 1664525) + & 
                            seed(j) + j ! Non-linear
         twister%state(i) = iand(twister%state(i), -1) ! for >32 bit machines
      end if
    end do
    
    !
    ! Walk through the state array, beginning where we left off in the block above
    ! 
    do i = mod(nfirstloop, blocksize) + nwraps + 1, blocksize - 1
      twister%state(i) = ieor(twister%state(i),                                 &
                              ieor(twister%state(i-1),                          & 
                                   ishft(twister%state(i-1), -30)) * 1566083941) - i ! Non-linear
      twister%state(i) = iand(twister%state(i), -1) ! for >32 bit machines
    end do
    
    twister%state(0) = twister%state(blocksize - 1) 
    
    do i = 1, mod(nfirstloop, blocksize) + nwraps
      twister%state(i) = ieor(twister%state(i),                                 &
                              ieor(twister%state(i-1),                          & 
                                   ishft(twister%state(i-1), -30)) * 1566083941) - i ! Non-linear
      twister%state(i) = iand(twister%state(i), -1) ! for >32 bit machines
    end do
    
    twister%state(0) = umask 
    twister%currentElement = blocksize
    
  end function initialize_vector
  ! -------------------------------------------------------------
  ! Public functions
  ! --------------------
  function getrandomint(twister)
    type(randomnumbersequence), intent(inout) :: twister
    integer                      :: getrandomint
    ! Generate a random integer on the interval [0,0xffffffff]
    !   Equivalent to genrand_int32 in the C code. 
    !   Fortran doesn't have a type that's unsigned like C does, 
    !   so this is integers in the range -2**31 - 2**31
    ! All functions for getting random numbers call this one, 
    !   then manipulate the result
    
    if(twister%currentElement >= blocksize) call nextstate(twister)
      
    getrandomint = temper(twister%state(twister%currentElement))
    twister%currentElement = twister%currentElement + 1
  
  end function getrandomint
  ! --------------------
  function getrandompositiveint(twister)
    type(randomnumbersequence), intent(inout) :: twister
    integer                      :: getrandompositiveint
    ! Generate a random integer on the interval [0,0x7fffffff]
    !   or [0,2**31]
    !   Equivalent to genrand_int31 in the C code. 
    
    ! Local integers
    integer :: localint

    localint = getrandomint(twister)
    getrandompositiveint = ishft(localint, -1)
  
  end function getrandompositiveint
  ! --------------------
  function getrandomreal(twister)
    type(randomnumbersequence), intent(inout) :: twister
    double precision             :: getrandomreal
    ! Generate a random number on [0,1]
    !   Equivalent to genrand_real1 in the C code
    !   The result is stored as double precision but has 32 bit resolution
    
    integer :: localint
    
    localint = getrandomint(twister)
    if(localint < 0) then
      getrandomreal = dble(localint + 2.0d0**32)/(2.0d0**32 - 1.0d0)
    else
      getrandomreal = dble(localint            )/(2.0d0**32 - 1.0d0)
    end if
  end function getrandomreal
  ! --------------------
  subroutine finalize_randomnumbersequence(twister)
    type(randomnumbersequence), intent(inout) :: twister
    
      twister%currentElement = blocksize
      twister%state(:) = 0
  end subroutine finalize_randomnumbersequence
  ! --------------------  
end module mersennetwister_mod