* A C-program for MT19937: Integer     version
*  genrand() generates one pseudorandom unsigned integer (32bit)
* which is uniformly distributed among 0 to 2^32-1  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).
*   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()      -> integer 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.
*
************************************************************************
* Reverse Generator by Katsumi Hagita.  Jun. 10, 2000.
*
*   revrand() means integer functions of grnd() in reverse
*
************************************************************************
* this main() outputs first 1000 generated numbers
*             and the 1000 numbers in reverse of the order
      program main
*
      implicit integer(a-z)
*
* Period parameters
      parameter(N     =  624)
      parameter(NO    =  100000)
*
      dimension mtn(0:NO)
      dimension mtr(0:NO)

      dimension mt(0:N-1)
*                     the array for the state vector
      dimension i(0:7)
*
      common /block/mti,mt
      save   /block/
*
      call sgrnd(4357)
*                         any nonzero integer can be used as a seed
      do 1000 j=0,no
CC        i(mod(j,8))=grnd()
        mtn(j)=grnd()
        i(mod(j,8))=mtn(j)
        if(mod(j,8).eq.7) then
CCCC          write(*,'(8(i12,'' ''))') (i(k),k=0,7)
        else if(j.eq.no) then
CCCC          write(*,'(8(i12,'' ''))') (i(k),k=0,mod(no-1,8))
        endif
 1000 continue
*

      do 1100 j=0,no
CC        i(mod(j,8))=revgrnd()
        mtr(no-j)=revgrnd()
        i(mod(j,8))=mtr(no-j)
        if(mod(j,8).eq.7) then
CCCC          write(*,'(8(i12,'' ''))') (i(k),k=0,7)
        else if(j.eq.no) then
CCCC          write(*,'(8(i12,'' ''))') (i(k),k=0,mod(no-1,8))
        endif
 1100 continue
      do 1200 j=0,no
         if( mtn(j)-mtr(j).ne.0) then
         write(6,*) j,mtn(j),mtr(j), mtn(j)-mtr(j)
         endif
 1200 continue

      stop
      end
************************************************************************
      subroutine sgrnd(seed)
*
      implicit integer(a-z)
*
* Period parameters
      parameter(N     =  624)
*
      dimension mt(0:N-1)
*                     the array for the state vector
      common /block/mti,mt
      save   /block/
*
*      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]
*
      mt(0)= iand(seed,-1)
      do 1000 mti=1,N-1
        mt(mti) = iand(69069 * mt(mti-1),-1)
 1000 continue
*
      return
      end
************************************************************************
      integer function grnd()
*
      implicit integer(a-z)
*
* Period parameters
      parameter(N     =  624)
      parameter(N1    =  N+1)
      parameter(M     =  397)
      parameter(MATA  = -1727483681)
*                                    constant vector a
      parameter(UMASK = -2147483648)
*                                    most significant w-r bits
      parameter(LMASK =  2147483647)
*                                    least significant r bits
* Tempering parameters
      parameter(TMASKB= -1658038656)
      parameter(TMASKC= -272236544)
*
      dimension mt(0:N-1)
*                     the array for the state vector
      common /block/mti,mt
      save   /block/
      data   mti/N1/
*                     mti==N+1 means mt[N] is not initialized
*
      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.N) then
*                       generate N words at one time
        if(mti.eq.N+1) then
*                            if sgrnd() has not been called,
          call sgrnd(4357)
*                              a default initial seed is used
        endif
*
        do 1000 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)))
 1000   continue
        do 1100 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)))
 1100   continue
        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
      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))
*
      grnd=y
*
      return
      end


**********************************
**********************************

      integer function revgrnd()
*
      implicit integer(a-z)
*
* Period parameters
      parameter(N     =  624)
      parameter(N1    =  N+1)
      parameter(M     =  397)
      parameter(MATA  = -1727483681)
*                                    constant vector a
      parameter(UMASK = -2147483648)
*                                    most significant w-r bits
      parameter(LMASK =  2147483647)
*                                    least significant r bits
      parameter(YMASK =  -2)
* Tempering parameters
      parameter(TMASKB= -1658038656)
      parameter(TMASKC= -272236544)
*
      dimension mt(0:N-1)
*                     the array for the state vector
      common /block/mti,mt
      save   /block/
***      data   mti/N1/
*                     mti==N+1 means mt[N] is not initialized
*
      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)
*

      mti=mti-1
      y=mt(mti)

      if(mti.eq.0) then

          z=ieor(mt(N-1),mt(M-1))
          x=ishft(z,-31)
          y=ieor(z,mag01(x))
          p=ior(iand(ishft(y,1),YMASK),x)
          mt0L=iand(LMASK,p)
          mt(0)=ior(iand(mt(0),UMASK),mt0L)
          mt0=mt(0)

       do 2000 kk=N-1,N-M+1,-1
          z=ieor(mt(kk-1),mt(kk-1+M-N))
          x=ishft(z,-31)
          y=ieor(z,mag01(x))
          q=ior(iand(ishft(y,1),YMASK),x)
          mt(kk)=ior(iand(UMASK,p),iand(LMASK,q))
          p=q
2000   continue

       do 2010 kk=N-M,1,-1
          z=ieor(mt(kk-1),mt(kk-1+M))
          x=ishft(z,-31)
          y=ieor(z,mag01(x))
          q=ior(iand(ishft(y,1),YMASK),x)
          mt(kk)=ior(iand(UMASK,p),iand(LMASK,q))
          p=q
2010   continue

          mt(0)=ior(iand(UMASK,p),mt0L)

        y=mt0

        mti = N
      endif
*

      y=ieor(y,TSHFTU(y))
      y=ieor(y,iand(TSHFTS(y),TMASKB))
      y=ieor(y,iand(TSHFTT(y),TMASKC))
      y=ieor(y,TSHFTL(y))
*
      revgrnd=y

      return
      end

