!-----------------------------------------------------------------------
!
!  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