!!!!!!!!! OBTAINED AT http://www.math.sci.hiroshima-u.ac.jp/~m-mat/MT/VERSIONS/FORTRAN/fortran.html
!!!!!!!!! VERSION IS A restorable MT19937 code in fortran 90 by Richard Woloshyn: http://www.math.sci.hiroshima-u.ac.jp/~m-mat/MT/VERSIONS/FORTRAN/mtfort90.f
!____________________________________________________________________________
! 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.
!
! genrand() -> double precision function grnd()
! sgenrand(seed) -> subroutine sgrnd(seed)
! integer seed
!
! 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 version rewritten as an F90 module and mt state saving and getting
! subroutines added by Richard Woloshyn. (rwww@triumf.ca). June 30, 1999
module mtmod
! Default seed
integer, parameter :: defaultsd = 4357
! Period parameters
integer, parameter :: Nmt = 624, N1 = Nmt + 1
! the array for the state vector
integer, save, dimension(0:Nmt-1) :: mt
integer, save :: mti = N1
! Overload procedures for saving and getting mt state
interface mtsave
module procedure mtsavef
module procedure mtsaveu
end interface
interface mtget
module procedure mtgetf
module procedure mtgetu
end interface
contains
!Initialization subroutine
subroutine sgrnd(seed)
implicit none
!
! setting 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]
!
integer, intent(in) :: seed
mt(0) = iand(seed,-1)
do mti=1,Nmt-1
mt(mti) = iand(69069 * mt(mti-1),-1)
enddo
!
return
end subroutine sgrnd
!Random number generator
real(8) function grnd()
implicit integer(a-z)
! Period parameters
integer, parameter :: M = 397, MATA = -1727483681
! constant vector a
integer, parameter :: LMASK = 2147483647
! least significant r bits
integer, parameter :: UMASK = -LMASK - 1
! most significant w-r bits
! Tempering parameters
integer, parameter :: TMASKB= -1658038656, TMASKC= -272236544
dimension mag01(0:1)
data mag01/0, MATA/
save mag01
! mag01(x) = x * MATA for x=0,1
TSHFTU(y)=ishft(y,-11)
TSHFTS(y)=ishft(y,7)
TSHFTT(y)=ishft(y,15)
TSHFTL(y)=ishft(y,-18)
if(mti.ge.Nmt) then
! generate N words at one time
if(mti.eq.Nmt+1) then
! if sgrnd() has not been called,
call sgrnd( defaultsd )
! a default initial seed is used
endif
do kk=0,Nmt-M-1
y=ior(iand(mt(kk),UMASK),iand(mt(kk+1),LMASK))
mt(kk)=ieor(ieor(mt(kk+M),ishft(y,-1)),mag01(iand(y,1)))
enddo
do kk=Nmt-M,Nmt-2
y=ior(iand(mt(kk),UMASK),iand(mt(kk+1),LMASK))
mt(kk)=ieor(ieor(mt(kk+(M-Nmt)),ishft(y,-1)),mag01(iand(y,1)))
enddo
y=ior(iand(mt(Nmt-1),UMASK),iand(mt(0),LMASK))
mt(Nmt-1)=ieor(ieor(mt(M-1),ishft(y,-1)),mag01(iand(y,1)))
mti = 0
endif
y=mt(mti)
mti = mti + 1
y=ieor(y,TSHFTU(y))
y=ieor(y,iand(TSHFTS(y),TMASKB))
y=ieor(y,iand(TSHFTT(y),TMASKC))
y=ieor(y,TSHFTL(y))
if(y .lt. 0) then
grnd=(dble(y)+2.0d0**32)/(2.0d0**32-1.0d0)
else
grnd=dble(y)/(2.0d0**32-1.0d0)
endif
return
end function grnd
!State saving subroutines.
! Usage: call mtsave( file_name, format_character )
! or call mtsave( unit_number, format_character )
! where format_character = 'u' or 'U' will save in unformatted form, otherwise
! state information will be written in formatted form.
subroutine mtsavef( fname, forma )
!NOTE: This subroutine APPENDS to the end of the file "fname".
character(*), intent(in) :: fname
character, intent(in) :: forma
select case (forma)
case('u','U')
open(unit=10,file=trim(fname),status='UNKNOWN',form='UNFORMATTED', &
position='APPEND')
write(10)mti
write(10)mt
case default
open(unit=10,file=trim(fname),status='UNKNOWN',form='FORMATTED', &
position='APPEND')
write(10,*)mti
write(10,*)mt
end select
close(10)
return
end subroutine mtsavef
subroutine mtsaveu( unum, forma )
integer, intent(in) :: unum
character, intent(in) :: forma
select case (forma)
case('u','U')
write(unum)mti
write(unum)mt
case default
write(unum,*)mti
write(unum,*)mt
end select
return
end subroutine mtsaveu
!State getting subroutines.
! Usage: call mtget( file_name, format_character )
! or call mtget( unit_number, format_character )
! where format_character = 'u' or 'U' will read in unformatted form, otherwise
! state information will be read in formatted form.
subroutine mtgetf( fname, forma )
character(*), intent(in) :: fname
character, intent(in) :: forma
select case (forma)
case('u','U')
open(unit=10,file=trim(fname),status='OLD',form='UNFORMATTED')
read(10)mti
read(10)mt
case default
open(unit=10,file=trim(fname),status='OLD',form='FORMATTED')
read(10,*)mti
read(10,*)mt
end select
close(10)
return
end subroutine mtgetf
subroutine mtgetu( unum, forma )
integer, intent(in) :: unum
character, intent(in) :: forma
select case (forma)
case('u','U')
read(unum)mti
read(unum)mt
case default
read(unum,*)mti
read(unum,*)mt
end select
return
end subroutine mtgetu
end module mtmod
! program main
! this main() outputs first 1000 generated numbers
!
! use mtmod
! implicit none
!
! integer, parameter :: no=1000
! real(8), dimension(0:7) :: r
! integer j,k
!! real(8) grnd
!!
!! call sgrnd(4357)
!! any nonzero integer can be used as a seed
! do j=0,no-1
! r(mod(j,8))=grnd()
! if(mod(j,8).eq.7) then
! write(*,'(8(f9.6,'' ''))') (r(k),k=0,7)
! else if(j.eq.no-1) then
! write(*,'(8(f9.6,'' ''))') (r(k),k=0,mod(no-1,8))
! endif
! enddo
!
! end program main