#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 100 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. ! outn is number of random numbers. ! Uniformity test using Kolmogorov-Smirnov test use mt_var implicit none integer i integer*8 r,genrand32 integer :: nodisp = 0 character(32) :: argv integer :: fno1 =1 integer rmax integer rmin integer outn integer irseq(1:nofseq) real rseq(1:nofseq) real rseq2(0:nofseq) integer r2 real Knalpha real D !Kolmogorov-Smirnov test statistics outn=nofseq open(fno1, file='ranout.csv', status='replace') rmax=0 rmin=1000000 mag01(1) = 0 mag01(2) = MATRIX_A call init_genrand(INITSEED) 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 do i=1,outn rseq(i)=rseq(i)/rmax end do !Kolmogorov-Smirnov test. call kstest(outn,rseq,D,Knalpha) print *, 'D value,',D print *, 'Knalpha,',Knalpha close(fno1) end subroutine kstest(n,x,D,Knalpha) ! Kolmogorov-Smirnov test implicit none integer n !No. of uniform random numbers. real x(n) !Uniform random numbers in (0,1]. real D !Kolmogorov-Smirnov test statistics. integer i integer j integer icount real Knalpha !Critical value. real alpha integer k k=n alpha=0.05 !Significant level. D=0 do i=1,k icount=0 do j=1,n if (x(j)<=real(i)/k) then icount=icount+1 end if end do if (D= 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