function pex1tailed=Wilcoxrsexact(n,m,W,ranks);
% DOES NOT WORK YET!!!!!!, p values wrong for 10,10.
% Exact p values for Wilcoxon Rank Sum statistic
% function p1tailed=Wilcoxrsexact(n,m,W,ranks);
% Borrows shamelessly from Strausss' combvals.m
% Input: n size of smaller group
%        m size of larger group
%        W Wilcoxon signed rank statistic
%        ranks, actual ranks of n+m items if there are ties present.
% strauss's combvals.m:
% Help file for combvals:
% COMBVALS: Generates the combinations of n integers taken r at a time.  The 
%           number of such combinations is given by function nc=combin().  
%     Usage: c = combvals(n,r)
%           n = number of integers (1:n) to be combined.
%           r = number to be taken at a time (0 < r <= n).
%           -------------------------------------------------------
%           c = [nc x r] matrix of combinations.

% Based on ACM Algorithm 94, J. Kurtzberg, Comm. ACM, June 1962.
% RE Strauss, 12/18/98

% An exact conditional distribution with ties follows Hollander & Wolfe p. 115
if nargin<4
    ranks=1:n+m;
else
    ranks=sort(ranks);
end
ranks=ranks(:);
  N=n+m;
  r = n;

  % Strauss's combval lists combinations in c in lexicographic
  % order, thus the critical values for sum(C) are larger than
  % observed W. We can speed up the process by using
  % Wstar=min(W,n*(m+n+1)-W) and exiting loop when Wstar fails
  % to be less than critical value
  if W>W-n*(m+n+1);
      W=n*(m+n+1)-W;
      lowertail=logical(1);
  else
      lowertail=logical(0);
  end
      
  ncomb = nchoosek(N,r);  % built-in combination function.
  % c = zeros(ncomb,r);     % Don't need to store values.
  Tally=0;
  j = zeros(1,r);
  contloop=logical(1);
  i=0;
  while i<ncomb & contloop
    i=i+1;
    b = 1;
    endflag = 0;
    while(~endflag)
      if (j(b)>=b)
        a = j(b)-b-1;
        for l = 1:b
          j(l) = l+a;
        end;
        endflag = 1;  
      else
        if (b==r)
          for b = 1:r
            j(b) = N-r-1+b;
          end;
          endflag = 1;
        end;
        b = b+1;
      end;
    end;
    % c(i,:) = N-j(r:-1:1);
    c=N-j(r:-1:1);
    if sum(ranks(c))<=W
        Tally=Tally+1;
        kount=0;
    else
        kount=kount+1;
        if kount>(N/2)^5  % This exit flag could probably be improved.
            contloop=logical(0);
        end
    end
  end;
pex1tailed=Tally/ncomb;