Random subset

[Beliavsky]

Stdlib has procedures to generate random numbers between 1 and N. When sampling, one may also want random subsets -- k distinct integers between 1 and N. I suggest that stdlib include this. I wrote my own code to do this, by generating successive random integers between 1 and N, 1 and N-1, 1 and N-2 etc. and choosing the associated integer between 1 and N. For N = 100000 and k = 40000 John Burkardt's ksub_random2 code was about 200 times faster. Here is a driver code for his subset.f90 that runs in 0.05 second on my PC.

program xksub_random
  implicit none
  integer, parameter :: k = 40000, n = 100000, niter = 1
  integer            :: i,j,a(k),seed,ncount(n)
  write ( *, "(a)" ) ""
  write ( *, "(a)" ) "KSUB_RANDOM2_TEST"
  write ( *, "(a)" ) "  KSUB_RANDOM2 generates a random K subset of an N set."
  write ( *, "(a,i8)" ) "  Set size is N =    ", n
  write ( *, "(a,i8)" ) "  Subset size is K = ", k
  write ( *, "(a)" ) ""
  seed = 123456789
  do i = 1, niter
    ncount = 0   
    call ksub_random2(n,k,seed,a)
    do j=1,k      
      ncount(a(j)) = ncount(a(j)) + 1
    end do
    write (*,"(4(i0,1x),f0.4)" ) minval(a),maxval(a),sum(ncount),maxval(ncount),sum(dble(a))/k
!    print*,"a=",a
  end do
end program xksub_random

with output

KSUB_RANDOM2_TEST
  KSUB_RANDOM2 generates a random K subset of an N set.
  Set size is N =      100000
  Subset size is K =    40000

1 99999 40000 1 49931.3902