mersenne_twister.f

Go to the documentation of this file.

 
 
      module mersenne_twister
        private
!  Public declarations
        public random_stat
        public random_seed
        public random_setseed
        public random_number
        public random_number_f
        public random_gauss
        public random_gauss_f
        public random_index
        public random_index_f
!  Parameters
        integer,parameter:: n=624
        integer,parameter:: m=397
        integer,parameter:: mata=-1727483681 ! constant vector a
        integer,parameter:: umask=-2147483648 ! most significant w-r bits
        integer,parameter:: lmask =2147483647 ! least significant r bits
        integer,parameter:: tmaskb=-1658038656 ! tempering parameter
        integer,parameter:: tmaskc=-272236544 ! tempering parameter
        integer,parameter:: mag01(0:1)=(/0,mata/)
        integer,parameter:: iseed=4357
        integer,parameter:: nrest=n+6
!  Defined types
        type random_stat
          private
          integer:: mti=n+1
          integer:: mt(0:n-1)
          integer:: iset
          real:: gset
        end type
!  Saved data
        type(random_stat),save:: sstat
!  Overloaded interfaces
        interface random_setseed
          module procedure random_setseed_s
          module procedure random_setseed_t
        end interface
        interface random_number
          module procedure random_number_i
          module procedure random_number_s
          module procedure random_number_t
        end interface
        interface random_gauss
          module procedure random_gauss_i
          module procedure random_gauss_s
          module procedure random_gauss_t
        end interface
        interface random_index
          module procedure random_index_i
          module procedure random_index_s
          module procedure random_index_t
        end interface
!  All the subprograms
      contains
 
        subroutine random_seed(size,put,get,stat)
          implicit none
          integer,intent(out),optional:: size
          integer,intent(in),optional:: put(nrest)
          integer,intent(out),optional:: get(nrest)
          type(random_stat),intent(inout),optional:: stat
          if(present(size)) then     ! return size of seed array
!     if(present(put).or.present(get))&
!      call errmsg('RANDOM_SEED: more than one option set - some ignored')
            size=nrest
          elseif(present(put)) then  ! restore from seed array
!     if(present(get))&
!      call errmsg('RANDOM_SEED: more than one option set - some ignored')
            if(present(stat)) then
              stat%mti=put(1)
              stat%mt=put(2:n+1)
              stat%iset=put(n+2)
              stat%gset=transfer(put(n+3:nrest),stat%gset)
              if(stat%mti.lt.0.or.stat%mti.gt.n.or.any(stat%mt.eq.0).or.
     &         stat%iset.lt.0.or.stat%iset.gt.1) then
                call random_setseed_t(iseed,stat)
!         call errmsg('RANDOM_SEED: invalid seeds put - default seeds used')
              endif
            else
              sstat%mti=put(1)
              sstat%mt=put(2:n+1)
              sstat%iset=put(n+2)
              sstat%gset=transfer(put(n+3:nrest),sstat%gset)
              if(sstat%mti.lt.0.or.sstat%mti.gt.n.or.any(sstat%mt.eq.0)
     &         .or.sstat%iset.lt.0.or.sstat%iset.gt.1) then
                call random_setseed_t(iseed,sstat)
!         call errmsg('RANDOM_SEED: invalid seeds put - default seeds used')
              endif
            endif
          elseif(present(get)) then  ! save to seed array
            if(present(stat)) then
              if(stat%mti.eq.n+1) call random_setseed_t(iseed,stat)
              get(1)=stat%mti
              get(2:n+1)=stat%mt
              get(n+2)=stat%iset
              get(n+3:nrest)=transfer(stat%gset,get,nrest-(n+3)+1)
            else
              if(sstat%mti.eq.n+1) call random_setseed_t(iseed,sstat)
              get(1)=sstat%mti
              get(2:n+1)=sstat%mt
              get(n+2)=sstat%iset
              get(n+3:nrest)=transfer(sstat%gset,get,nrest-(n+3)+1)
            endif
          else                       ! reset default seed
            if(present(stat)) then
              call random_setseed_t(iseed,stat)
            else
              call random_setseed_t(iseed,sstat)
            endif
          endif
        end subroutine
 
        subroutine random_setseed_s(inseed)
          implicit none
          integer,intent(in):: inseed
          call random_setseed_t(inseed,sstat)
        end subroutine
 
        subroutine random_setseed_t(inseed,stat)
          implicit none
          integer,intent(in):: inseed
          type(random_stat),intent(out):: stat
          integer ii,mti
          ii=inseed
          if(ii.eq.0) ii=iseed
          stat%mti=n
          stat%mt(0)=iand(ii,-1)
          do mti=1,n-1
            stat%mt(mti)=iand(69069*stat%mt(mti-1),-1)
          enddo
          stat%iset=0
          stat%gset=0.
        end subroutine
 
        function random_number_f() result(harvest)
          implicit none
          real:: harvest
          real h(1)
          if(sstat%mti.eq.n+1) call random_setseed_t(iseed,sstat)
          call random_number_t(h,sstat)
          harvest=h(1)
        end function
 
        subroutine random_number_i(harvest,inseed)
          implicit none
          real,intent(out):: harvest(:)
          integer,intent(in):: inseed
          type(random_stat) stat
          call random_setseed_t(inseed,stat)
          call random_number_t(harvest,stat)
        end subroutine
 
        subroutine random_number_s(harvest)
          implicit none
          real,intent(out):: harvest(:)
          if(sstat%mti.eq.n+1) call random_setseed_t(iseed,sstat)
          call random_number_t(harvest,sstat)
        end subroutine
 
        subroutine random_number_t(harvest,stat)
          implicit none
          real,intent(out):: harvest(:)
          type(random_stat),intent(inout):: stat
          integer j,kk,y
          integer tshftu,tshfts,tshftt,tshftl
          tshftu(y)=ishft(y,-11)
          tshfts(y)=ishft(y,7)
          tshftt(y)=ishft(y,15)
          tshftl(y)=ishft(y,-18)
          do j=1,size(harvest)
            if(stat%mti.ge.n) then
              do kk=0,n-m-1
                y=ior(iand(stat%mt(kk),umask),iand(stat%mt(kk+1),lmask))
                stat%mt(kk)=ieor(ieor(stat%mt(kk+m),ishft(y,-1)),
     &   mag01(iand(y,1)))
              enddo
              do kk=n-m,n-2
                y=ior(iand(stat%mt(kk),umask),iand(stat%mt(kk+1),lmask))
                stat%mt(kk)=ieor(ieor(stat%mt(kk+(m-n)),ishft(y,-1)),
     &   mag01(iand(y,1)))
              enddo
              y=ior(iand(stat%mt(n-1),umask),iand(stat%mt(0),lmask))
              stat%mt(n-1)=ieor(ieor(stat%mt(m-1),ishft(y,-1)),
     &   mag01(iand(y,1)))
              stat%mti=0
            endif
            y=stat%mt(stat%mti)
            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
              harvest(j)=(real(y)+2.0**32)/(2.0**32-1.0)
            else
              harvest(j)=real(y)/(2.0**32-1.0)
            endif
            stat%mti=stat%mti+1
          enddo
        end subroutine
 
        function random_gauss_f() result(harvest)
          implicit none
          real:: harvest
          real h(1)
          if(sstat%mti.eq.n+1) call random_setseed_t(iseed,sstat)
          call random_gauss_t(h,sstat)
          harvest=h(1)
        end function
 
        subroutine random_gauss_i(harvest,inseed)
          implicit none
          real,intent(out):: harvest(:)
          integer,intent(in):: inseed
          type(random_stat) stat
          call random_setseed_t(inseed,stat)
          call random_gauss_t(harvest,stat)
        end subroutine
 
        subroutine random_gauss_s(harvest)
          implicit none
          real,intent(out):: harvest(:)
          if(sstat%mti.eq.n+1) call random_setseed_t(iseed,sstat)
          call random_gauss_t(harvest,sstat)
        end subroutine
 
        subroutine random_gauss_t(harvest,stat)
          implicit none
          real,intent(out):: harvest(:)
          type(random_stat),intent(inout):: stat
          integer mx,my,mz,j
          real r2(2),r,g1,g2
          mz=size(harvest)
          if(mz.le.0) return
          mx=0
          if(stat%iset.eq.1) then
            mx=1
            harvest(1)=stat%gset
            stat%iset=0
          endif
          my=(mz-mx)/2*2+mx
          do
            call random_number_t(harvest(mx+1:my),stat)
            do j=mx,my-2,2
              call rgauss(harvest(j+1),harvest(j+2),r,g1,g2)
              if(r.lt.1.) then
                harvest(mx+1)=g1
                harvest(mx+2)=g2
                mx=mx+2
              endif
            enddo
            if(mx.eq.my) exit
          enddo
          if(my.lt.mz) then
            do
              call random_number_t(r2,stat)
              call rgauss(r2(1),r2(2),r,g1,g2)
              if(r.lt.1.) exit
            enddo
            harvest(mz)=g1
            stat%gset=g2
            stat%iset=1
          endif
        contains
 
          subroutine rgauss(r1,r2,r,g1,g2)
            real,intent(in):: r1,r2
            real,intent(out):: r,g1,g2
            real v1,v2,fac
            v1=2.*r1-1.
            v2=2.*r2-1.
            r=v1**2+v2**2
            if(r.lt.1.) then
              fac=sqrt(-2.*log(r)/r)
              g1=v1*fac
              g2=v2*fac
            endif
          end subroutine
        end subroutine
 
        function random_index_f(imax) result(iharvest)
          implicit none
          integer,intent(in):: imax
          integer:: iharvest
          integer ih(1)
          if(sstat%mti.eq.n+1) call random_setseed_t(iseed,sstat)
          call random_index_t(imax,ih,sstat)
          iharvest=ih(1)
        end function
 
        subroutine random_index_i(imax,iharvest,inseed)
          implicit none
          integer,intent(in):: imax
          integer,intent(out):: iharvest(:)
          integer,intent(in):: inseed
          type(random_stat) stat
          call random_setseed_t(inseed,stat)
          call random_index_t(imax,iharvest,stat)
        end subroutine
 
        subroutine random_index_s(imax,iharvest)
          implicit none
          integer,intent(in):: imax
          integer,intent(out):: iharvest(:)
          if(sstat%mti.eq.n+1) call random_setseed_t(iseed,sstat)
          call random_index_t(imax,iharvest,sstat)
        end subroutine
 
        subroutine random_index_t(imax,iharvest,stat)
          implicit none
          integer,intent(in):: imax
          integer,intent(out):: iharvest(:)
          type(random_stat),intent(inout):: stat
          integer,parameter:: mh=n
          integer i1,i2,mz
          real h(mh)
          mz=size(iharvest)
          do i1=1,mz,mh
            i2=min((i1-1)+mh,mz)
            call random_number_t(h(:i2-(i1-1)),stat)
            iharvest(i1:i2)=max(ceiling(h(:i2-(i1-1))*imax),1)
          enddo
        end subroutine
      end module