#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(rmaxr2) 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