function(y=c(14 , 18 , 2 , 4 , -5 , 14 , -3 , -1 , 1 , 6 , 3 , 3),
         x=c(8 , 26 , -7 , -1 , 2 , 9 , 0 , -4 , 13 , 3 , 3 , 4),
         alternative="greater",flag=T,PDF=F){
# flag = T means that the absolute differences are ranked including the zero
#          differences, but then the zeros are deleted.
# flag = F means that the zero differences are deleted and then the remaining
#        absolute differences are ranked.
#
if(PDF==T) pdf("SignedRankExact.pdf",width=11)
D=y-x
absD=abs(D)
D.sign=sign(D)
rabsD=rank(absD)
rabsD=rabsD[D.sign!=0]
D.sign=D.sign[D.sign!=0]
if(flag==F){
  rabsD=rank(rabsD) 
  # This reranks the ranks of absolute differences after
  # the ranks corresponding to zero differences are removed.
  # This amounts to the same as removing zero differences first
  # and then ranking the remaining differences.
}
r.sign=rabsD*D.sign
Vs.star.obs=sum(r.sign[r.sign>0])
Vs.star=0
mVs.star=sum(rabsD)/2 # mean of Vs.star
sigVs.star=sqrt(sum(rabsD^2)/4) # standard dev. of Vs.star
N=length(rabsD)
for( k in 1:floor((N-1)/2)){
Vs.star=c(Vs.star,combn(rabsD,k,FUN=sum))
}
Vs.star=c(Vs.star,2*mVs.star-Vs.star) #this uses the symmetry of the Vs.star
                                      # distribution around the mean mVs.star
if(N%%2==0) Vs.star=c(Vs.star,combn(rabsD,N/2,FUN=sum))
NN=length(Vs.star)
if(alternative=="greater"){
p.val=mean(Vs.star>=Vs.star.obs)}else{
p.val=mean(Vs.star<=Vs.star.obs)
}
xm=mVs.star-3.5*sigVs.star
xM=mVs.star+3.5*sigVs.star
xp=sort(unique(Vs.star))
dx=diff(xp)
k=length(xp)
xxp=c(xp[1]-dx[1]/2,xp[1:(k-1)]+dx/2,xp[k]+dx[k-1]/2)
hist.out=hist(Vs.star,breaks=xxp,probability=T,
  col=c("grey","grey30"),xlim=c(xm,xM),main="",xlab=expression(V[s]^'*'))

abline(h=0)
abline(v=Vs.star.obs,lwd=2,col="red")
if(alternative=="greater"){
text(Vs.star.obs+.05*sigVs.star,.95*max(hist.out$density),
   substitute(P[H](V[s]^'*'>=vx)==px,list(vx=Vs.star.obs,px=round(p.val,4))),adj=0)}else{
text(Vs.star.obs-.05*sigVs.star,.95*max(hist.out$density),
   substitute(P[H](V[s]^'*'<=vx)==px,list(vx=Vs.star.obs,px=round(p.val,4))),adj=1)
}
x=seq(xm,xM,length.out=250)
y=dnorm(x,mVs.star,sigVs.star)
lines(x,y,col="blue")
Z=(Vs.star.obs-mVs.star)/sigVs.star
if(alternative=="greater"){
p.val.normal=1-pnorm(Z)}else{
p.val.normal=pnorm(Z)
}
out=c(Vs.star.obs,mVs.star,sigVs.star,Z,p.val,p.val.normal)
names(out)=c("Vs.star.obs","meanVs.star","sigVs.star","Z","p.val","p.val.normal")
if(PDF==T) dev.off()
out
}