#define P 624
#define Q 397
#define MATRIX_A X'9908b0df'
#define UPPER_MASK X'80000000'
#define LOWER_MASK X'7fffffff'
#define MULTIPLIER 1664525
#define INITSEED 19660809
#define MASK X'ffffffff'
#define nofseq 10000
#define N0 100
#define ns0 50
module mt_var
    implicit none
    integer*8 mt(1:P)
    integer :: mti = P + 2
    integer*8 mag01(1:2)
end module

program main
!   Generation of uniform random numbers by Mersenne Twister. 
!   Rndom walk test in 2 dimensional.
    use mt_var
    implicit none
    integer*8 r,genrand32

    integer :: nodisp = 0
    character(32) :: argv

    integer rmax
    integer rmin
    integer outn
    integer irseq(1:nofseq)
    real rseq(1:nofseq)
    integer r2
    
    real chi2       !Chi-square statistics
    integer df      !Degrees of freedom
    integer k       !Dimension of random walk test
    integer N       !No. of iterations
    integer ns      !No. of steps for 1 iteration
    integer freq(1:2,1:2)   !No. of random numbers in each quadrant.
    
    integer i
    integer j
    
    real x(N0,ns0,2)
    
    outn=nofseq

    rmax=0
    rmin=1000000

    mag01(1) = 0
    mag01(2) = MATRIX_A

    call init_genrand(INITSEED)
    !Generate random numbers.
    do i=1,outn
        r = genrand32()
        r2=ishft(r, -1)
        if(rmax<r2) then
            rmax=r2
        end if
        if(rmin>r2) then
            rmin=r2
        end if           
        irseq(i)=r2
        rseq(i)=r2
    end do
    rmax=rmax+10000
    !Converted to [0,1) uniform random numbers.
    do i=1,outn
        rseq(i)=rseq(i)/rmax
    end do

    k=2
    N=N0
    ns=outn/(k*N)
    k=1
    do i=1,N
        do j=1,ns
           x(i,j,1)=rseq(k) 
           x(i,j,2)=rseq(k+1) 
           k=k+2
        end do
    end do

    !2-dimensional random walk test
    call ranwalk(N0,ns,x,chi2,df,freq)
    print *, 'chi-square value,',chi2
    print *, 'degrees of freedom,',df
    print *, 'No. of frequencies in each quadrant'
    do i=1,2
        print *, (freq(i,j),j=1,2)
    end do

end
    
    
subroutine ranwalk(N,ns,x,chi2,df,freq)
    !2-dimensional random walk test
    implicit none
    integer N               !No. of iterations
    integer ns              !No. of steps for 1 iteration
    real x(1:N,1:ns,1:2) !Uniform random numbers [0,1].
    real chi2               !Chi-square statistics
    integer df              !Degrees of freedom
    integer freq(1:2,1:2)   !No. of random numbers in each quadrant.
    integer i
    integer j
    integer xpos
    integer ypos
    do i=1,2
        do j=1,2
            freq(i,j)=0
        end do
    end do
    do i=1,N
        xpos=0
        ypos=0
        do j=1,ns
            if (x(i,j,1).lt.0.5) then
                xpos=xpos-1
            else
                xpos=xpos+1
            end if
            if (x(i,j,2).lt.0.5) then
                ypos=ypos-1
            else
                ypos=ypos+1
            end if
        end do
        if (xpos.ge.0.and.ypos.ge.0) then
            freq(2,2)=freq(2,2)+1
        else if (xpos.lt.0.and.ypos.ge.0) then
            freq(1,2)=freq(1,2)+1
        else if (xpos.ge.0.and.ypos.lt.0) then
            freq(2,1)=freq(2,1)+1
        else
            freq(1,1)=freq(1,1)+1
        end if
    end do
               
    chi2=0.0
    do i=1,2
        do j=1,2
            chi2=chi2+(freq(i,j)-N/4.0)**2/(N/4.0)
        end do
    end do
    df=2**2-1
end

subroutine init_genrand(s)
    use mt_var
    implicit none
    integer s

    mt(1) = iand(s, MASK)
    do mti=2,P
        mt(mti) = MULTIPLIER * mt(mti - 1) + 1
        mt(mti) = iand(mt(mti), MASK)
    end do
end

function genrand32()
    use mt_var
    implicit none
    integer*8 genrand32
    integer*8 y
    integer :: v1 = 5489
    integer*8 :: v2 = X'9d2c5680'
    integer*8 :: v3 = X'efc60000'
    integer kk

    if(mti >= P+1) then
        if(mti == P+2) call init_genrand(v1)
        do kk=1,P-Q
            y = ior(iand(mt(kk), UPPER_MASK), iand(mt(kk+1), LOWER_MASK))
            mt(kk) = ieor(ieor(mt(kk+Q), ishft(y,-1)), mag01(iand(y, 1)+1))
        end do
        do kk=P-Q+1,P-1
            y = ior(iand(mt(kk), UPPER_MASK), iand(mt(kk+1), LOWER_MASK))
            mt(kk) = ieor(ieor(mt(kk+Q-P), ishft(y,-1)), mag01(iand(y, 1)+1))
        end do
        y = ior(iand(mt(P), UPPER_MASK), iand(mt(1), LOWER_MASK))
        mt(P) = ieor(ieor(mt(Q), ishft(y,-1)), mag01(iand(y, 1)+1))
        mti = 1
    end if
    y = mt(mti)
    mti = mti + 1
    y = ieor(y, ishft(y, -11))
    y = ieor(y, iand(ishft(y, 7), v2))
    y = ieor(y, iand(ishft(y, 15), v3))
    y = ieor(y, ishft(y, -18))
    genrand32 = y
end