module modrandom
  implicit none
  ! Period parameters.
  integer, parameter :: n=624, m=397
  ! Constant vector A.
  integer, parameter :: mata=-1727483681
  ! Most significant w-r bits, note the way to set its value.
  integer, parameter :: umask=-2147483647-1
  ! Least significant r bits.
  integer, parameter :: lmask=2147483647
  ! Tempering parameters.
  integer, parameter :: tmaskb=-1658038656
  integer, parameter :: tmaskc=-272236544
  ! mag01(x)=x*mata for x=0,1.
  integer, parameter :: mag01(0:1)=(/0,mata/)
  ! State vector and index, mti=n+1 means st vector is
  ! not initialized.
  integer, save :: mt(0:n-1)=0, mti=n+1
end module modrandom
subroutine initrand(seed)
  use modrandom
  implicit none
  integer, intent(in) :: seed
  mt(0)=iand(seed,-1)
  do mti=1,n-1
     ! See Knuth's TAOCP vol2 for multiplier.
     mt(mti)=1812433253*ieor(mt(mti-1),ishft(mt(mti-1),-30))+mti
     mt(mti)=iand(mt(mti),-1)
  end do
  mti=mti+1
end subroutine initrand

subroutine genrand(rnd)
  use modrandom
  implicit none
  real(8), intent(out) :: rnd
  integer :: y, kk
  ! Generate N words at once.
  if (mti.gt.n) then
     ! If seed not set, use the default seed 5489.
     if (mti.eq.n+1) call initrand(5489)
     do kk=0,n-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)))
     end do
     do kk=n-m,n-2
        y=ior(iand(mt(kk),umask),iand(mt(kk+1),lmask))
        mt(kk)=ieor(ieor(mt(kk+m-n),ishft(y,-1)),mag01(iand(y,1)))
     end do
     y=ior(iand(mt(n-1),umask),iand(mt(0),lmask))
     mt(n-1)=ieor(ieor(mt(m-1),ishft(y,-1)),mag01(iand(y,1)))
     mti=0
  end if
  ! Tempering
  y=mt(mti)
  mti=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))
  ! Convert to double in [0, 1] interval
  if (y .lt. 0) then
     rnd=(dble(y)+2.d0**32)/(2.d0**32-1.d0)
  else
     rnd=dble(y)/(2.d0**32-1.d0)
  end if
end subroutine genrand

program test
i=6
j=10
x = ran2(i)
y = ran2(j)

write(*,*) x, y

end program test


FUNCTION ran2(idum)

INTEGER, INTENT(inOUT)                     :: idum

REAL :: ran2
INTEGER, PARAMETER :: im1=2147483563
INTEGER, PARAMETER :: im2=2147483399
REAL, PARAMETER :: am=1./im1
INTEGER, PARAMETER :: imm1=im1-1
INTEGER, PARAMETER :: ia1=40014
INTEGER, PARAMETER :: ia2=40692
INTEGER, PARAMETER :: iq1=53668
INTEGER, PARAMETER :: iq2=52774
INTEGER, PARAMETER :: ir1=12211
INTEGER, PARAMETER :: ir2=3791
INTEGER, PARAMETER :: ntab=32
INTEGER, PARAMETER :: ndiv=1+imm1/ntab
REAL, PARAMETER :: eps=1.2E-7
REAL, PARAMETER :: rnmx=1.-eps
INTEGER :: idum2,j,k,iv(ntab),iy
SAVE iv,iy,idum2
DATA idum2/123456789/, iv/ntab*0/, iy/0/

IF (idum <= 0) THEN
  idum=MAX(-idum,1)
  idum2=idum
  DO  j=ntab+8,1,-1
    k=idum/iq1
    idum=ia1*(idum-k*iq1)-k*ir1
    IF (idum < 0) idum=idum+im1
    IF (j <= ntab) iv(j)=idum
  END DO
  iy=iv(1)
END IF
k=idum/iq1
idum=ia1*(idum-k*iq1)-k*ir1
IF (idum < 0) idum=idum+im1
k=idum2/iq2
idum2=ia2*(idum2-k*iq2)-k*ir2
IF (idum2 < 0) idum2=idum2+im2
j=1+iy/ndiv
iy=iv(j)-idum2
iv(j)=idum
IF(iy < 1)iy=iy+imm1
print *, am*iy, rnmx
ran2=MIN(am*iy,rnmx)

RETURN
END FUNCTION ran2