!----------------------------------------------------------------------- ! ! test90.f90 ! ! Test driver for the Fortran 90 portable random number library ! random_pl.f90, implementing a number of tests recommended by ! Knuth in The Art of Computer Programming, v. 2: Seminumerical ! Algorithms. ! !----------------------------------------------------------------------- program test90 use random_pl use globals implicit none integer, parameter :: nlong = 10000000 ! number of deviates for long tests integer, parameter :: nshort = 1000 ! number of deviates for short tests integer, parameter :: nrep = nlong/nshort ! times to repeat short tests real(wp), dimension(nlong) :: rarr real(wp), dimension(nshort, nrep) :: sarr equivalence ( rarr, sarr ) real(wp) :: Knp, Knm real(wp), dimension(nrep) :: Knparr, Knmarr character(128) :: fmtp="(' Kn+ = ',F8.5,', 90% ci lo = ',F8.5,', hi = ',F8.5)" character(128) :: fmtm="(' Kn- = ',F8.5,', 90% ci lo = ',F8.5,', hi = ',F8.5)" integer :: j integer, parameter :: serial_d = 100 integer, dimension(serial_d,serial_d) :: serial_bins integer, dimension(serial_d**2) :: serial_bins_lin equivalence ( serial_bins, serial_bins_lin ) integer :: indx1, indx2 integer :: gap_t integer, parameter :: gap_n = nlong integer, dimension(:), allocatable :: gap_count real(wp), dimension(:), allocatable :: gap_p integer :: r, s, tt real(wp) :: tmp, alpha, beta integer, dimension(5) :: poker_arr, poker_count integer :: itmp integer :: poker_d logical :: newval integer, parameter :: coupon_n = nlong integer, dimension(5:20) :: coupon_count integer, dimension(0:4) :: coupon_occurs integer :: q integer, dimension(120) :: perm_count real(wp), dimension(5) :: perm_stor integer :: t, tfac, f integer, dimension(:), allocatable :: run_stor real(wp), dimension(:), allocatable :: run_prob integer :: run_len real(wp) :: run_last, run_current integer, parameter :: run_n = nlong integer :: max_n integer, dimension(:), allocatable :: max_count integer, dimension(6) :: coll_num real, dimension(6) :: coll_pct real, dimension(2**14) :: coll_a logical, dimension(2**20) :: coll_count real(wp), dimension(20) :: coll_vec integer(rpk) :: coll_indx integer :: m, n, j0, j1, collisions real(wp) :: p character(128) :: fmtc="(' coll = ',I8,', 90% ci lo = ',I8,', hi = ',I8)" integer(rpk), dimension(512) :: birthdays, spacings integer, dimension(4) :: birth_count integer :: eq_spacings real(wp) :: corr, corr_mu, corr_sig character(128) ::fmtcr="(' corr = ',F8.5,', 90% ci lo = ',F8.5,', hi = ',F8.5)" write(0,*) write(0,*) write(0,*) '-------------------------------------------------------------------' write(0,*) 'Beginning random number tests' write(0,*) write(0,*) 'Tests are among those recommended by Knuth' write(0,*) write(0,*) 'Each test will report "passed", "marginally suspect", "suspect",' write(0,*) ' or "failed". Keep in mind that even for a good random' write(0,*) ' number generator, there should be the following frequencies:' write(0,*) write(0,*) ' passed -> 80%' write(0,*) ' marginally suspect -> 10%' write(0,*) ' suspect -> 8%' write(0,*) ' failed -> 2%' write(0,*) write(0,*) 'The purpose of this code is to look for consistent failures' write(0,*) ' across distinct runs with different seeds.' write(0,*) '-------------------------------------------------------------------' write(0,*) write(0,*) !!! use default seeds !!! get deviates !call init_rand_pl(173737_rpk, 673_rpk, 62283_rpk) !call init_rand_pl(87654321_rpk, 31415927_rpk, 12345678_rpk) !call init_rand_pl(54321_rpk, 12345_rpk, 23456_rpk) write(0,*) 'Getting random deviates:' write(0,*) write(0,*) rarr = rand_pl_vec(nlong) !call random_number(harvest=rarr) !!! do Kolmogorov-Smirnov tests write(0,*) 'Performing Kolmogorov-Smirnov test on ', nlong,' numbers:' call kstest(rarr, Knp, Knm, 1) if ( Knp.lt.ks_limit(0.01, nlong) .or. Knp.gt.ks_limit(0.99, nlong) ) then write(0,*) 'Kn+ statistic: FAILED' write(0,fmtp) Knp, ks_limit(0.05, nlong), ks_limit(0.95, nlong) else if ( Knp.lt.ks_limit(0.05, nlong) .or. Knp.gt.ks_limit(0.95, nlong)) then write(0,*) 'Kn+ statistic: SUSPECT' write(0,fmtp) Knp, ks_limit(0.05, nlong), ks_limit(0.95, nlong) else if ( Knp.lt.ks_limit(0.1, nlong) .or. Knp.gt.ks_limit(0.9, nlong)) then write(0,*) 'Kn+ statistic: MARGINALLY SUSPECT' write(0,fmtp) Knp, ks_limit(0.05, nlong), ks_limit(0.95, nlong) else write(0,*) 'Kn+ statistic: passed' end if if ( Knm.lt.ks_limit(0.01, nlong) .or. Knm.gt.ks_limit(0.99, nlong) ) then write(0,*) 'Kn- statistic: FAILED' write(0,fmtm) Knm, ks_limit(0.05, nlong), ks_limit(0.95, nlong) else if ( Knm.lt.ks_limit(0.05, nlong) .or. Knm.gt.ks_limit(0.95, nlong)) then write(0,*) 'Kn- statistic: SUSPECT' write(0,fmtm) Knm, ks_limit(0.05, nlong), ks_limit(0.95, nlong) else if ( Knm.lt.ks_limit(0.1, nlong) .or. Knm.gt.ks_limit(0.9, nlong)) then write(0,*) 'Kn- statistic: MARGINALLY SUSPECT' write(0,fmtm) Knm, ks_limit(0.05, nlong), ks_limit(0.95, nlong) else write(0,*) 'Kn- statistic: passed' end if write(0,*) 'Performing Kolmogorov-Smirnov test on ', nlong,& ' numbers in subsets of ', nshort, ':' ! perform subset tests do j = 1, nrep call kstest(sarr(:,j), Knparr(j), Knmarr(j), 1) end do ! perform tests on resulting statistics write(0,*) 'testing Kn+ statistics: ' call kstest(Knparr, Knp, Knm, 2) if ( Knp.lt.ks_limit(0.01, nrep) .or. Knp.gt.ks_limit(0.99, nrep) ) then write(0,*) 'Kn+ statistic: FAILED' write(0,fmtp) Knp, ks_limit(0.05, nrep), ks_limit(0.95, nrep) else if ( Knp.lt.ks_limit(0.05, nrep) .or. Knp.gt.ks_limit(0.95, nrep)) then write(0,*) 'Kn+ statistic: SUSPECT' write(0,fmtp) Knp, ks_limit(0.05, nrep), ks_limit(0.95, nrep) else if ( Knp.lt.ks_limit(0.1, nrep) .or. Knp.gt.ks_limit(0.9, nrep)) then write(0,*) 'Kn+ statistic: MARGINALLY SUSPECT' write(0,fmtp) Knp, ks_limit(0.05, nrep), ks_limit(0.95, nrep) else write(0,*) 'Kn+ statistic: passed' end if if ( Knm.lt.ks_limit(0.01, nrep) .or. Knm.gt.ks_limit(0.99, nrep) ) then write(0,*) 'Kn- statistic: FAILED' write(0,fmtm) Knm, ks_limit(0.05, nrep), ks_limit(0.95, nrep) else if ( Knm.lt.ks_limit(0.05, nrep) .or. Knm.gt.ks_limit(0.95, nrep)) then write(0,*) 'Kn- statistic: SUSPECT' write(0,fmtm) Knm, ks_limit(0.05, nrep), ks_limit(0.95, nrep) else if ( Knm.lt.ks_limit(0.1, nrep) .or. Knm.gt.ks_limit(0.9, nrep)) then write(0,*) 'Kn- statistic: MARGINALLY SUSPECT' write(0,fmtm) Knm, ks_limit(0.05, nrep), ks_limit(0.95, nrep) else write(0,*) 'Kn- statistic: passed' end if write(0,*) 'testing Kn- statistics: ' call kstest(Knmarr, Knp, Knm, 2) if ( Knp.lt.ks_limit(0.01, nrep) .or. Knp.gt.ks_limit(0.99, nrep) ) then write(0,*) 'Kn+ statistic: FAILED' write(0,fmtp) Knp, ks_limit(0.05, nrep), ks_limit(0.95, nrep) else if ( Knp.lt.ks_limit(0.05, nrep) .or. Knp.gt.ks_limit(0.95, nrep)) then write(0,*) 'Kn+ statistic: SUSPECT' write(0,fmtp) Knp, ks_limit(0.05, nrep), ks_limit(0.95, nrep) else if ( Knp.lt.ks_limit(0.1, nrep) .or. Knp.gt.ks_limit(0.9, nrep)) then write(0,*) 'Kn+ statistic: MARGINALLY SUSPECT' write(0,fmtp) Knp, ks_limit(0.05, nrep), ks_limit(0.95, nrep) else write(0,*) 'Kn+ statistic: passed' end if if ( Knm.lt.ks_limit(0.01, nrep) .or. Knm.gt.ks_limit(0.99, nrep) ) then write(0,*) 'Kn- statistic: FAILED' write(0,fmtm) Knm, ks_limit(0.05, nrep), ks_limit(0.95, nrep) else if ( Knm.lt.ks_limit(0.05, nrep) .or. Knm.gt.ks_limit(0.95, nrep)) then write(0,*) 'Kn- statistic: SUSPECT' write(0,fmtm) Knm, ks_limit(0.05, nrep), ks_limit(0.95, nrep) else if ( Knm.lt.ks_limit(0.1, nrep) .or. Knm.gt.ks_limit(0.9, nrep)) then write(0,*) 'Kn- statistic: MARGINALLY SUSPECT' write(0,fmtm) Knm, ks_limit(0.05, nrep), ks_limit(0.95, nrep) else write(0,*) 'Kn- statistic: passed' end if write(0,*) write(0,*) !!! do serial test write(0,*) 'Performing serial test with ', serial_d, '^2 bins:' serial_bins = 0 do j = 1, nlong-1, 2 indx1 = ceiling(rarr(j )*serial_d) indx2 = ceiling(rarr(j+1)*serial_d) serial_bins(indx1, indx2) = serial_bins(indx1, indx2) + 1 end do call chisqtest(serial_bins_lin, (/(1._wp/serial_d**2, j=1, serial_d**2)/)) write(0,*) write(0,*) !!! do gap test ! compute probabilities for chi-squared test gap_t = log(10._wp/gap_n)/log(0.5_wp) allocate( gap_count(0:gap_t), gap_p(0:gap_t) ) do j = 0, gap_t-1 gap_p(j) = 0.5_wp**(j+1) end do gap_p(gap_t) = 0.5_wp**gap_t do tt = 1, 2 if ( tt .eq. 1 ) then write(0,*) 'Performing gap test (runs above the mean), ', gap_t, ' bins:' alpha = 0 beta = 0.5 else write(0,*) 'Performing gap test (runs below the mean), ', gap_t, ' bins:' alpha = 0.5 beta = 1 end if s = 0 gap_count = 0 do while ( s .lt. gap_n ) r = 0 do while ( .true. ) tmp = rand_pl() if ( tmp .ge. alpha .and. tmp .lt. beta ) then gap_count(min(gap_t,r)) = gap_count(min(gap_t,r)) + 1 exit end if r = r + 1 end do s = s + 1 end do call chisqtest(gap_count, gap_p) end do deallocate( gap_count, gap_p ) write(0,*) write(0,*) !!! do poker test poker_d = (nlong/50)**0.25 write(0,*) 'Performing classical unordered poker test, ', poker_d,' bins:' poker_count = 0 do j = 1, nlong, 5 poker_arr = ceiling(rarr(j:(j+4)) * poker_d) itmp = 1 do s = 2, 5 newval = .true. do r = 1, s-1 if ( poker_arr(r) .eq. poker_arr(s) ) newval = .false. end do if ( newval ) itmp = itmp + 1 end do poker_count(itmp) = poker_count(itmp) + 1 end do call chisqtest(poker_count, & (/ 1.0_wp/poker_d**4, & (1.0_wp*(poker_d-1))/poker_d**4 * 15, & (1.0_wp*(poker_d-1))*(poker_d-2)/poker_d**4 * 25, & (1.0_wp*(poker_d-1))*(poker_d-2)*(poker_d-3)/poker_d**4 * 10, & (1.0_wp*(poker_d-1))*(poker_d-2)*(poker_d-3)*(poker_d-4)/poker_d**4 & /) ) write(0,*) write(0,*) !!! do coupon collector's test write(0,*) 'Performing coupon collector test, 5 coupons, up to 20 tries' s = 0 coupon_count = 0 do while ( s .lt. coupon_n ) q = 0 r = 0 coupon_occurs = 0 1003 continue r = r + 1 itmp = rand_pl() * 5 if ( coupon_occurs(itmp) .ne. 0 ) go to 1003 coupon_occurs(itmp) = 1 q = q + 1 if ( q .lt. 5 ) go to 1003 coupon_count(min(r,20)) = coupon_count(min(r,20)) + 1 s = s + 1 end do if ( 0.018*coupon_n .lt. 10 ) then write(0,*) 'This test may fail due to insufficient number' end if call chisqtest(coupon_count, & (/ 0.0384_wp, & ! r=5 0.0768_wp, & ! r=6 0.09984_wp, & ! r=7 0.10752_wp, & ! r=8 0.10450944_wp, & ! r=9 0.09547776_wp, & ! r=10 0.083816448_wp, & ! r=11 0.07163904_wp, & ! r=12 0.060112994304_wp, & ! r=13 0.049791565824_wp, & ! r=14 0.0408620040192_wp, & ! r=15 0.0333100744704_wp, & ! r=16 0.0270216306622_wp, & ! r=17 0.0218419625460_wp, & ! r=18 0.0176085709986_wp, & ! r=19 0.0714485091740_wp & ! r>=20 /) ) write(0,*) write(0,*) !!! do permutation test do t = 3, 5 write(0,*) 'Performing permutation test, ', t, ' elements at a time:' tfac = 1 do j = 2, t tfac = tfac * j end do perm_count = 0 perm_stor = 0 do j = 1, nlong-t, t perm_stor(1:t) = rarr(j:(j+t-1)) r = t f = 0 2002 continue s = maxloc(perm_stor(1:r), dim=1) f = r*f + s - 1 tmp = perm_stor(s) perm_stor(s) = perm_stor(r) perm_stor(r) = tmp r = r - 1 if ( r .gt. 1 ) go to 2002 perm_count(f+1) = perm_count(f+1) + 1 end do call chisqtest(perm_count(1:tfac), (/ (1._wp/tfac, j=1, tfac) /) ) end do write(0,*) write(0,*) !!! do run test t = 1 do while ( nlong/(sqrt(6.28_wp*t)*(1._wp*t)**t*exp(-1._wp*t)) .gt. 10 ) t = t + 1 end do t = t - 1 write(0,*) 'Performing simple run test (ascending), length ', t, ' max:' allocate( run_stor( t ), run_prob(t) ) run_prob = 0 run_prob(1) = 1 do j = 2, t tfac = 1 do s = 2, j tfac = tfac * s end do run_prob(j) = run_prob(j) + 1._wp/tfac run_prob(j-1) = run_prob(j-1) - 1._wp/tfac end do run_prob(t) = 1./tfac j = 1 run_len = 1 run_stor = 0 s = 0 run_last = rand_pl() do while ( s .lt. run_n ) run_current = rand_pl() if ( run_current .gt. run_last ) then run_len = run_len + 1 run_last = run_current else run_stor(min(run_len,t)) = run_stor(min(run_len,t)) + 1 run_len = 1 run_last = rand_pl() run_last = rand_pl() s = s + 1 end if end do call chisqtest(run_stor, run_prob) write(0,*) 'Performing simple run test (descending), length ', t, ' max:' j = 1 run_len = 1 run_stor = 0 s = 0 run_last = rand_pl() do while ( s .le. run_n ) run_current = rand_pl() if ( run_current .lt. run_last ) then run_len = run_len + 1 run_last = run_current else run_stor(min(run_len,t)) = run_stor(min(run_len,t)) + 1 run_len = 1 run_last = rand_pl() run_last = rand_pl() s = s + 1 end if end do call chisqtest(run_stor, run_prob) deallocate( run_stor, run_prob ) write(0,*) write(0,*) !!! do max-of-t tests max_n = nlong/50 allocate( max_count(max_n) ) write(0,*) 'Performing max-of-t test for t = 2, ', max_n, ' bins:' max_count = 0 do j = 1, nlong/2 indx1 = ceiling(maxval(rarr((2*j-1):(2*j)))**2*max_n) max_count(indx1) = max_count(indx1) + 1 end do call chisqtest(max_count, (/ (1._wp/max_n, j=1, max_n) /) ) write(0,*) 'Performing max-of-t test for t = 5, ', max_n, ' bins:' max_count = 0 do j = 1, nlong/5 indx1 = ceiling(maxval(rarr((5*j-4):(5*j)))**5*max_n) max_count(indx1) = max_count(indx1) + 1 end do call chisqtest(max_count, (/ (1._wp/max_n, j=1, max_n) /) ) deallocate( max_count ) write(0,*) write(0,*) !!! do collision test, only on a relatively small sample write(0,*) 'Collision test, 20 dimensions, 2**20 bins, 2**14 vectors:' coll_count = .false. collisions = 0 n = 2**14 m = 2**20 do j = 1, n coll_vec = rand_pl_vec(20) coll_indx = 0 do s = 1, 20 coll_indx = ishft(coll_indx, 1) coll_indx = coll_indx + int(coll_vec(s)*2) end do if ( coll_count(coll_indx) ) then collisions = collisions + 1 else coll_count(coll_indx) = .true. end if end do coll_a = 0 coll_a(2) = 1 j0 = 1 j1 = 1 do s = 1, n-1 j1 = j1 + 1 do j = j1, j0, -1 coll_a(j+1) = real(j,kind=wp)/m*coll_a(j+1) + & ((1+1._wp/m)-real(j,kind=wp)/m)*coll_a(j) if ( coll_a(j) .lt. 1.e-20 ) then coll_a(j) = 0 if ( j .eq. j1 ) j1 = j1 - 1 if ( j .eq. j0 ) j0 = j0 + 1 end if end do end do p = 0 t = 1 j = j0 - 1 coll_pct = (/ 0.01, 0.05, 0.1, 0.9, 0.95, 0.99 /) do t = 1, size(coll_pct) do while (p .le. coll_pct(t)) j = j + 1 p = p + coll_a(j) end do coll_num(7-t) = n - j - 1 end do if ( collisions.lt.coll_num(1) .or. collisions.gt.coll_num(6) ) then write(0,*) 'collision statistic: FAILED' write(0,fmtc) collisions, coll_num(2), coll_num(5) else if ( collisions.lt.coll_num(2) .or. collisions.gt.coll_num(5)) then write(0,*) 'collision statistic: SUSPECT' write(0,fmtc) collisions, coll_num(2), coll_num(5) else if ( collisions.lt.coll_num(3) .or. collisions.gt.coll_num(4)) then write(0,*) 'collision statistic: MARGINALLY SUSPECT' write(0,fmtc) collisions, coll_num(2), coll_num(5) else write(0,*) 'collision statistic: passed' end if write(0,*) write(0,*) !!! do birthday spacings test write(0,*) 'Performing birthday spacings test, m=2**25, n=512, 3 max eq spcngs:' m = 2**25 n = 512 birth_count = 0 do j = 1, nlong-n, n birthdays = rarr(j:(j+n-1)) * m call isort(birthdays) spacings(1:(n-1)) = birthdays(2:n) - birthdays(1:(n-1)) spacings(n) = birthdays(1) + m - birthdays(n) call isort(spacings) eq_spacings = 0 do s = 2, n if ( spacings(s) .eq. spacings(s-1) ) eq_spacings = eq_spacings + 1 end do birth_count(min(3,eq_spacings)+1) = birth_count(min(3,eq_spacings)+1) + 1 end do call chisqtest(birth_count, & (/ 0.368801577_wp, 0.369035243_wp, 0.183471182_wp, 0.078691997_wp /) ) write(0,*) write(0,*) !!! do serial correlation test write(0,*) 'Performing serial correlation test:' corr = (nlong*dot_product(rarr,cshift(rarr,1)) - sum(rarr)**2) / & (nlong*dot_product(rarr,rarr) - sum(rarr)**2) corr_mu = -1._wp/(nlong-1) corr_sig = (1._wp*nlong)/(n-1._wp)/sqrt(n-2._wp) if ( corr.lt.corr_mu-2.34*corr_sig .or. corr.gt.corr_mu+2.34*corr_sig ) then write(0,*) 'correlation statistic: FAILED' write(0,fmtcr) corr, corr_mu-1.64*corr_sig, corr_mu+1.64*corr_sig else if ( corr.lt.corr_mu-1.64*corr_sig .or. corr.gt.corr_mu+1.64*corr_sig )then write(0,*) 'correlation statistic: SUSPECT' write(0,fmtcr) corr, corr_mu-1.64*corr_sig, corr_mu+1.64*corr_sig else if ( corr.lt.corr_mu-1.28*corr_sig .or. corr.gt.corr_mu+1.28*corr_sig )then write(0,*) 'correlation statistic: MARGINALLY SUSPECT' write(0,fmtcr) corr, corr_mu-1.64*corr_sig, corr_mu+1.64*corr_sig else write(0,*) 'correlation statistic: passed' end if write(0,*) write(0,*) contains !----------------------------------------------------------------------- ! ! subroutine kstest(arr, Knp, Knm, functionflag) ! ! Implements k-s test on random deviates, given input array arr ! output is the two statistics Knp and Knm. The integer functionflag ! means to perform the test on uniform random deviates on (0, 1) ! if functionflag = 1 (i.e., F(x) = x), and if functionflag = 2, test ! for the distribution function F(x) = 1 - exp(-2*x*x), useful for ! doing a K-S test on the Kn+ and Kn- statistics. ! !----------------------------------------------------------------------- subroutine kstest(arr, Knp, Knm, functionflag) implicit none real(wp), dimension(:), intent(in) :: arr real(wp), intent(out) :: Knp, Knm integer, intent(in) :: functionflag real(wp), dimension(size(arr)) :: cpy, a, b integer, dimension(size(arr)) :: c integer :: j, k, n, m real(wp) :: y, rp, rm n = size(arr) m = n a = 1 b = 0 c = 0 cpy = arr do j = 1, n if ( functionflag .eq. 1 ) y = cpy(j) if ( functionflag .eq. 2 ) y = 1._wp - exp(-2*cpy(j)*cpy(j)) k = ceiling(m * y) a(k) = min(a(k), y) b(k) = max(b(k), y) c(k) = c(k) + 1 end do j = 0 rp = 0 rm = 0 do k = 1, m if ( c(k) .gt. 0 ) then rm = max(rm, a(k) - real(j,kind=wp)/n) j = j + c(k) rp = max(rp, real(j,kind=wp)/n - b(k)) end if end do Knp = sqrt(real(n,kind=wp)) * rp Knm = sqrt(real(n,kind=wp)) * rm end subroutine kstest !----------------------------------------------------------------------- ! ! real function ks_limit(p, n) ! ! Calculates percentage points for the distributions of Kn+ and Kn- ! for the K-S test. Asymptotic formula is used here, with an ! accuracy of O(n**(-3/2)). Here p is the percentage (expressed as a ! decimal, not percentage) and n is the number of observations in ! the test. ! !----------------------------------------------------------------------- function ks_limit(p, n) implicit none real, intent(in) :: p integer, intent(in) :: n real :: ks_limit ks_limit = sqrt(-0.5*log(1.-p)) - sqrt(1./n)/6. & -(4./9.*p*p-2./3.)*p*p/n end function ks_limit !----------------------------------------------------------------------- ! ! subroutine chisqtest(arr, parr) ! ! Implements chi-squared test on random deviates in a 1-D integer ! array, where the corresponding probabilites are given in the real(wp) ! array parr. ! ! Also checks statistical consistency, using approximate formula ! for chi-squared integrals given by Knuth. ! !----------------------------------------------------------------------- subroutine chisqtest(arr, parr) implicit none integer, dimension(:), intent(in) :: arr real(wp), dimension(:), intent(in) :: parr integer :: k, n real(wp) :: chisq, nu, lim99, lim95, lim90, lim10, lim5, lim1, xp character(128) :: & fmt="(' chisq = ',F12.3,', 90% ci lo = ',F12.3,', hi = ',F12.3)" k = size(arr) n = sum(arr) chisq = sum(real(arr,kind=wp)**2/parr)/n - n nu = k - 1 xp = 2.3263 lim99 = nu + sqrt(2*nu)*xp + 2._wp/3._wp*(xp*xp - 1) xp = 1.6449 lim95 = nu + sqrt(2*nu)*xp + 2._wp/3._wp*(xp*xp - 1) xp = 1.2816 lim90 = nu + sqrt(2*nu)*xp + 2._wp/3._wp*(xp*xp - 1) xp = -1.2816 lim10 = nu + sqrt(2*nu)*xp + 2._wp/3._wp*(xp*xp - 1) xp = -1.6449 lim5 = nu + sqrt(2*nu)*xp + 2._wp/3._wp*(xp*xp - 1) xp = -2.3263 lim1 = nu + sqrt(2*nu)*xp + 2._wp/3._wp*(xp*xp - 1) if ( chisq.lt.lim1 .or. chisq.gt.lim99 ) then write(0,*) 'chisq statistic: FAILED' write(0,fmt) chisq, lim5, lim95 else if ( chisq.lt.lim5 .or. chisq.gt.lim95 ) then write(0,*) 'chisq statistic: SUSPECT' write(0,fmt) chisq, lim5, lim95 else if ( chisq.lt.lim10 .or. chisq.gt.lim90 ) then write(0,*) 'chisq statistic: MARGINALLY SUSPECT' write(0,fmt) chisq, lim5, lim95 else write(0,*) 'chisq statistic: passed' end if end subroutine chisqtest !----------------------------------------------------------------------- ! ! subroutine isort(p) ! ! Sorts an array of integers p. ! This is an implementation of the merge sort routine. ! ! Uses subroutine irecursivesort ! !----------------------------------------------------------------------- subroutine isort(p) implicit none integer, dimension(:), intent(inout) :: p integer :: n n = size(p) call irecursivesort(p) end subroutine isort recursive subroutine irecursivesort(p) implicit none integer, dimension(:), intent(inout) :: p integer, dimension(1:size(p)) :: pc integer n, n2, ptr, ptr1, ptr2 ! a list of length 1 doesn't need sorting n = size(p) if ( n .lt. 2 ) return ! otherwise, split p into two lists and sort them n2 = n/2 call irecursivesort(p(1:n2)) call irecursivesort(p((n2+1):n)) ! combine the two sorted lists ptr = 1 ptr1 = 1 ptr2 = 1 + n2 do while ( ptr .le. n ) if ( ptr1 .gt. n2 ) then pc(ptr) = p(ptr2) ptr2 = ptr2 + 1 else if ( ptr2 .gt. n ) then pc(ptr) = p(ptr1) ptr1 = ptr1 + 1 else if ( p(ptr1) .le. p(ptr2) ) then pc(ptr) = p(ptr1) ptr1 = ptr1 + 1 else pc(ptr) = p(ptr2) ptr2 = ptr2 + 1 end if end if ptr = ptr + 1 end do p = pc end subroutine irecursivesort end program test90