Generate a random number from a mt19937_type
Type | Intent | Optional | Attributes | Name | ||
---|---|---|---|---|---|---|
class(mt19937_type), | intent(inout) | :: | this |
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