#
# set data with "nt" treatment obs followed by "nt" control obs
#
n <- 6
yobs <- c(3,5,0,4,0,1)
nt <- n/2
#
# compute observed test statistic
#
sumy <- sum(yobs)
tobs <- sum(yobs[1:nt])/nt - sum(yobs[(nt+1):n])/nt
#
# set up vector to hold randomization distribution of test stat
#
out <- rep(0,choose(n,nt))
cnt <- 0
#
# loop to identify unique randomizations
#
for (i1 in (1:(nt+1))) {
   for (i2 in ((i1+1):(nt+2))) {
      for (i3 in ((i2+1):n)) {
#
# compute test statistic for randomization and store
#
         cnt <- cnt + 1
         trtsum <- sum(yobs[c(i1,i2,i3)])
         out[cnt] <- trtsum/nt - (sumy-trtsum)/nt
         }}}
#
# compute p-value 
#
pval <- sum(abs(out) >= abs(tobs))/choose(n,nt)

#
# simulation version using the R sample command to choose the sample
#
outsim <- rep(0,1000)
for (i in (1:1000)) {
   ysamp <- sample(yobs,nt)
   trtsum <- sum(ysamp)
   outsim[i] <- trtsum/nt - (sumy-trtsum)/nt
   }
pvalsim <- sum(abs(outsim) >= abs(tobs))/1000