mersennetwister.F90

!***********************************************************************
!*                   GNU Lesser General Public License
!*
!* This file is part of the GFDL Flexible Modeling System (FMS).
!*
!* FMS is free software: you can redistribute it and/or modify it under
!* the terms of the GNU Lesser General Public License as published by
!* the Free Software Foundation, either version 3 of the License, or (at
!* your option) any later version.
!*
!* FMS is distributed in the hope that it will be useful, but WITHOUT
!* ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
!* FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
!* for more details.
!*
!* You should have received a copy of the GNU Lesser General Public
!* License along with FMS.  If not, see <http://www.gnu.org/licenses/>.
!***********************************************************************
 
!> @defgroup mersennetwister_mod MersenneTwister_mod
!> @ingroup random_numbers
!> @brief Fortran-95 implementation of the Mersenne Twister 19937 algorithm
!> @author Robert Pincus
!!
!! 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
 
!> @addtogroup mersennetwister_mod
!> @{
module mersennetwister_mod
! -------------------------------------------------------------
  use platform_mod, only: r8_kind, i8_kind
 
  implicit none
  private
 
  ! Algorithm parameters
  ! -------
  ! Period parameters
  integer, parameter :: blocksize = 624,                 &
                        m         = 397,                 &
                        matrix_a  = -1727483681,         & !< constant vector a         (0x9908b0dfUL)
                        umask     = -2147483648_i8_kind, & !< 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
  !> @ingroup mersennetwister_mod
  type randomnumbersequence
    integer                            :: currentElement ! = blockSize
    integer, dimension(0:blockSize -1) :: state ! = 0
  end type randomnumbersequence
 
  !> @ingroup mersennetwister_mod
  interface new_randomnumbersequence
    module procedure initialize_scalar, initialize_vector
  end interface new_randomnumbersequence
 
  public :: randomnumbersequence
  public :: new_randomnumbersequence, finalize_randomnumbersequence, &
            getrandomint, getrandompositiveint, getrandomreal
!> @addtogroup mersennetwister_mod
!> @{
! -------------------------------------------------------------
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
  ! --------------------
   !> 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
  function getrandomint(twister)
    type(randomnumbersequence), intent(inout) :: twister
    integer                      :: getrandomint
 
    if(twister%currentElement >= blocksize) call nextstate(twister)
 
    getrandomint = temper(twister%state(twister%currentElement))
    twister%currentElement = twister%currentElement + 1
 
  end function getrandomint
  ! --------------------
  !> Generate a random integer on the interval [0,0x7fffffff]
  !! or [0,2**31].
  !! Equivalent to genrand_int31 in the C code.
  function getrandompositiveint(twister)
    type(randomnumbersequence), intent(inout) :: twister
    integer                      :: getrandompositiveint
 
    ! Local integers
    integer :: localint
 
    localint = getrandomint(twister)
    getrandompositiveint = ishft(localint, -1)
 
  end function getrandompositiveint
  ! --------------------
  !> 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
  function getrandomreal(twister)
    type(randomnumbersequence), intent(inout) :: twister
    real(r8_kind)                             :: getrandomreal
 
    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
!> @}
! close documentation grouping