mt19937.f90 Source File

!> Mersenne-Twister pseudorandom number generator (PRNG), specifically
!> [MT19937](https://en.wikipedia.org/wiki/Mersenne_Twister)
!>
!> This module has a type, [[mt19937_type]], which can be used to have multiple,
!> independent PRNGs within the same program. There is a module-level instance
!> that can be used via the public procedures [[sgrnd]] and [[grnd]], which can
!> be used to have a consistent, reproducible PRNG across a whole program.
!>
!> The original documentation comment follows:
!>
!>  A C-program for MT19937: Real number version
!>    genrand() generates one pseudorandom real number (double)
!>  which is uniformly distributed on [0,1]-interval, for each
!>  call. sgenrand(seed) set initial values to the working area
!>  of 624 words. Before genrand(), sgenrand(seed) must be
!>  called once. (seed is any 32-bit integer except for 0).
!>  Integer generator is obtained by modifying two lines.
!>    Coded by Takuji Nishimura, considering the suggestions by
!>  Topher Cooper and Marc Rieffel in July-Aug. 1997.
!>
!>  This library is free software; you can redistribute it and/or
!>  modify it under the terms of the GNU Library General Public
!>  License as published by the Free Software Foundation; either
!>  version 2 of the License, or (at your option) any later
!>  version.
!>  This library 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 Library General Public License for more details.
!>  You should have received a copy of the GNU Library General
!>  Public License along with this library; if not, write to the
!>  Free Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA
!>  02111-1307  USA
!>
!>  Copyright (C) 1997 Makoto Matsumoto and Takuji Nishimura.
!>  When you use this, send an email to: matumoto@math.keio.ac.jp
!>  with an appropriate reference to your work.
!>
!> ***********************************************************************
!>  Fortran translation by Hiroshi Takano.  Jan. 13, 1999.
!>
!>  This program uses the following non-standard intrinsics.
!>    ishft(i,n): If n>0, shifts bits in i by n positions to left.
!>                If n<0, shifts bits in i by n positions to right.
!>    iand (i,j): Performs logical AND on corresponding bits of i and j.
!>    ior  (i,j): Performs inclusive OR on corresponding bits of i and j.
!>    ieor (i,j): Performs exclusive OR on corresponding bits of i and j.
!>
!> ************************************************************************
!>  Fortran 95 translation and modularized to use in gyrokinetics project
!>  by R. Numata, June 2, 2010.
!>   * removed statement functions (TSHFTU,TSHFTS,TSHFTT,TSHFTL)
!>   * uses [0,1) interval
!>   * the above bit manipulation functions became standard.
!>   * definition of UMASK is corrected (the original value cannot be
!>     represented as a kind=4 integer.
!> ***********************************************************************
module mt19937

  implicit none

  private
  public :: sgrnd, grnd, mt19937_type

  integer, parameter :: default_seed=4357

  ! Period parameters
  integer, parameter :: N_mt19937 = 624 !< state size
  integer, parameter :: M_mt19937 = 397 !< state size

  !> Mersenne-Twister object. Allows multiple, independent instances
  !> to be used within a program
  type mt19937_type
    integer :: mt(0:N_mt19937-1) !< The array for the state vector
    integer :: mti = N_mt19937+1 !< mti==N+1 means mt[N] is not initialized
  contains
    procedure :: set_seed
    procedure :: generate
  end type mt19937_type

  !> Module-level instance
  type(mt19937_type) :: module_mt19937

contains

  !> Set initial seed of module-level [[mt19937_type]] instance
  subroutine sgrnd(seed)
    implicit none
    integer, intent(in) :: seed
    call module_mt19937%set_seed(seed)
  end subroutine sgrnd

  !> Set initial seed of a [[mt19937_type]] instance
  !>
  !> Set 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]
  subroutine set_seed(this, seed)
    class(mt19937_type), intent(inout) :: this
    integer, intent(in) :: seed
    this%mt(0) = iand(seed, -1)
    associate(mti => this%mti)
      do mti = 1, N_mt19937 - 1
        this%mt(mti) = iand(69069 * this%mt(mti - 1), -1)
      end do
    end associate
    ! Explicitly make sure this is set, as loop index may be undefined
    this%mti = N_mt19937
  end subroutine set_seed

  !> Generate a random number from the module-level [[mt19937_type]] instance
  function grnd ()
    implicit none
    real :: grnd
    grnd = module_mt19937%generate()
  end function grnd

  !> Generate a random number from a [[mt19937_type]]
  function generate(this)
    class(mt19937_type), intent(inout) :: this
    real :: generate

    ! Period parameters
    integer, parameter :: LMASK = 2147483647 ! least significant r bits
    integer, parameter :: UMASK = -LMASK - 1 ! most significant w-r bits
    integer, parameter :: MATA = -1727483681 !< constant vector a
    integer, parameter :: mag01(0:1) = [0, MATA] ! mag01(x) = x * MATA for x=0,1

    ! Tempering parameters
    integer, parameter :: TMASKB = -1658038656
    integer, parameter :: TMASKC = -272236544

    real, parameter :: pow = 2.**32       ! range of possible `integer` values
    real, parameter :: div = 1./pow       ! divided by 2**32 [0,1)-real-interval
    real :: random_temp
    integer :: y, kk

    if (this%mti >= N_mt19937) then ! generate N words at one time
       if (this%mti == N_mt19937+1) then ! if sgrnd() has not been called,
          call this%set_seed(default_seed) ! a default initial seed is used
       endif

      do kk = 0, N_mt19937 - M_mt19937 - 1
        y = ior(iand(this%mt(kk), UMASK), iand(this%mt(kk + 1), LMASK))
        this%mt(kk) = ieor(ieor(this%mt(kk + M_mt19937), ishft(y, -1)), mag01(iand(y, 1)))
      end do

      do kk = N_mt19937 - M_mt19937, N_mt19937 - 2
        y = ior(iand(this%mt(kk), UMASK), iand(this%mt(kk + 1), LMASK))
        this%mt(kk) = ieor(ieor(this%mt(kk + (M_mt19937 - N_mt19937)), ishft(y, -1)), mag01(iand(y, 1)))
      end do

      y = ior(iand(this%mt(N_mt19937 - 1), UMASK), iand(this%mt(0), LMASK))
      this%mt(N_mt19937 - 1) = ieor(ieor(this%mt(M_mt19937 - 1), ishft(y, -1)), mag01(iand(y, 1)))
      this%mti = 0
    end if

    y = this%mt(this%mti)
    this%mti = this%mti + 1
    y = ieor(y, ishft(y, -11))
    y = ieor(y, iand(ishft(y, 7), TMASKB))
    y = ieor(y, iand(ishft(y, 15), TMASKC))
    y = ieor(y, ishft(y, -18))

    random_temp = real(y)
    if (random_temp < 0.) random_temp = random_temp + pow
    generate = random_temp*div
  end function generate

end module mt19937