BlockedWilcoxon=function(datlist,alternative="greater",Nsim=10000,PDF=F){
#========================================================================
# This function needs as input a list of lists, say b lists.
# Each of the b lists should consist of two vectors x and y
# in that order, where y represents the treatment scores
# and x represents the control scores for the block 
# represented by this sublist.
# It is assumed that the random assigment of treatments 
# from block to block are independent.
# alternative="greater" means that we expect that the y-scores
# will tend to be larger than the x-scores under the alternative.
# Other values for alternative are "less" and "two.sided", with
# corresponding meanings.
# This function evaluates the exact null distribution and 
# exact p-value when the number of combined randomization
# combinations over all blocks does not exceed Nsim. Otherwise it 
# estimates the null distribution by Nsim simulations. 
# The time to run the simulations is propertional to Nsim. 
# For Nsim=1000000 it took 260 seconds
# using a slow 900MHz processor and 500MB RAM. 
#========================================================================
if(alternative!="greater" & alternative!="less" & alternative!="two.sided"){
return("alternative must be `greater' or `less' or `two.sided'")
}
if(PDF==T) pdf(file="BlockedWilcoxon.pdf",width=11)
b=length(datlist)
mean.W=0
var.W=0
W.obs=0
Nvec=NULL
nvec=NULL
MM=1
for(i in 1:b){
   dat.i=datlist[[i]]
   x=dat.i[[1]]
   y=dat.i[[2]]
   nvec=c(nvec,length(y))
   Nvec=c(Nvec,length(x)+length(y))
   MM=MM*choose(Nvec[i],nvec[i])
}
if(MM<=Nsim){
  Wvec=rep(0,MM)
  for(i in 1:b){
     dat.i=datlist[[i]]
     x=dat.i[[1]]
     y=dat.i[[2]]
     m=length(x)
     n=length(y)
     N.i=m+n
     R=rank(c(y,x))
     mean.W=mean.W+n*mean(R)/(N.i+1)
     var.W=var.W+var(R)*n*m/(N.i*(N.i+1)^2)
     W.obs=W.obs+sum(R[1:n])/(N.i+1)
     dist.new=combn(R,n,FUN=sum)/(N.i+1)
     if(i>1){
         Wvec=outer(Wvec,dist.new,"+")}else{
         Wvec=dist.new
     }
  }
}else{
  # computing W.obs and rank sets
  R.list=list()
  for(i in 1:b){
        dat.i=datlist[[i]]
        x=dat.i[[1]]
        y=dat.i[[2]]
        m=length(x)
        n=length(y)
        N.i=m+n
        R=rank(c(y,x))
        mean.W=mean.W+n*mean(R)/(N.i+1)
        var.W=var.W+var(R)*n*m/(N.i*(N.i+1)^2)
        R.list[[i]]=R
        W.obs=W.obs+sum(R[1:n])/(N.i+1)
  }
  # simulation loop
  Wvec=rep(0,Nsim)
  for(j in 1:Nsim){
       WW=0
       for(i in 1:b){
           WW=WW+sum(sample(R.list[[i]],nvec[i]))/(Nvec[i]+1)
       }
       Wvec[j]=WW  
  }
}
Wvec=round(Wvec,9)
sig.W=sqrt(var.W)
m=min(Wvec,mean.W-5*sig.W)
M=max(Wvec,mean.W+5*sig.W)
W.u=unique(sort(Wvec))
d.u=diff(W.u)
k=length(d.u)
br=c(W.u[1]-d.u[1]/2,W.u[1:k]+d.u/2,W.u[1+k]+d.u[k]/2)
if(MM<=Nsim){
out.hist=hist(Wvec,breaks=br,probability=T,xlim=c(m,M),col=c("grey","grey90"),
xlab="",main="Block Combined Wilcoxon Test")}else{
out.hist=hist(Wvec,nclass=151,probability=T,xlim=c(m,M),col=c("grey","grey90"),
xlab="",main="Block Combined Wilcoxon Test")
}
mtext(substitute(W[s]==sum(W[s]^"(i)"/(N[i]+1), i==1,bx),list(bx=b)),1,3.5)
if(alternative=="greater" | alternative=="less"){
  abline(v=W.obs,col="red")}else{
  offset=abs(W.obs-mean.W)
  abline(v=mean.W+c(-offset,offset),col="red")
}
if(alternative=="less"){
  pval=mean(Wvec<=W.obs)
  pval.norm=pnorm((W.obs-mean.W)/sig.W)
  if(MM<=Nsim){
    text(W.obs-.05*sig.W,max(out.hist$density),"exact",adj=1)}else{
    text(W.obs-.05*sig.W,max(out.hist$density),"estimated",adj=1)
    text(W.obs-.05*sig.W,.70*max(out.hist$density),substitute(N[sim]==xNsim,list(xNsim=Nsim)),adj=1)
  }
  text(W.obs-.05*sig.W,.95*max(out.hist$density),paste("p-value =",round(pval,4)),adj=1)
  text(W.obs-.05*sig.W,.85*max(out.hist$density),paste("normal approximation"),adj=1)
  text(W.obs-.05*sig.W,.80*max(out.hist$density),paste("p-value =",round(pval.norm,4)),adj=1)
}
if(alternative=="greater"){
  pval=mean(Wvec>=W.obs)
  pval.norm=1-pnorm((W.obs-mean.W)/sig.W)
  if(MM<=Nsim){
    text(W.obs+.05*sig.W,max(out.hist$density),"exact",adj=0)}else{
    text(W.obs+.05*sig.W,max(out.hist$density),"estimated",adj=0)
    text(W.obs+.05*sig.W,.70*max(out.hist$density),substitute(N[sim]==xNsim,list(xNsim=Nsim)),adj=0)
  }
  text(W.obs+.05*sig.W,.95*max(out.hist$density),paste("p-value =",round(pval,4)),adj=0)
  text(W.obs+.05*sig.W,.85*max(out.hist$density),paste("normal approximation"),adj=0)
  text(W.obs+.05*sig.W,.80*max(out.hist$density),paste("p-value =",round(pval.norm,4)),adj=0)
}
if(alternative=="two.sided"){

  pval=mean(abs(Wvec-mean.W)>=abs(W.obs-mean.W))
  pval.norm=2*pnorm(-abs(W.obs-mean.W)/sig.W)
  if(MM<=Nsim){
    text(W.obs+.05*sig.W,max(out.hist$density),"exact",adj=0)}else{
    text(W.obs+.05*sig.W,max(out.hist$density),"estimated",adj=0)
    text(W.obs+.05*sig.W,.70*max(out.hist$density),substitute(N[sim]==xNsim,list(xNsim=Nsim)),adj=0)
  }
  text(W.obs+.05*sig.W,.95*max(out.hist$density),paste("2-sided p-value =",round(pval,4)),adj=0)
  text(W.obs+.05*sig.W,.85*max(out.hist$density),paste("normal approximation"),adj=0)
  text(W.obs+.05*sig.W,.80*max(out.hist$density),paste("2-sided p-value =",round(pval.norm,4)),adj=0)
}
x=seq(m,M,length.out=201)
y=dnorm(x,mean.W,sig.W)
lines(x,y,col="blue")
if(PDF==T) dev.off()
}