% function f=schurpoly(l,x). Computes the Schur function corresponding
% to partition l with variables specified in x.
% 
% Plamen Koev, May 2002

function f=schurpoly(l,xx);
  global x sh Lp Lmax n lma
  x=xx;

  n=length(x);

  Lp=0;
  while (length(l)>Lp)&(l(Lp+1)>0)
     Lp=Lp+1;
  end % number of parts in l
  l=l(1:Lp);
  Lmax=l+1;

  sh(prod(Lmax),1,n)=-1; % allocate space

  if Lp>0
     lma(Lp)=1;
     for i=Lp-1:-1:1
        lma(i)=lma(i+1)*Lmax(i+1);
     end

     initialize(1,1,0);

     if isequal(class(xx),'sym')
        sh=sym(sh);
     end

     f=sch(1,n,1,nmu(l));
  else
     f=1;
  end

  clear x Lmax sh Lp n lma;

function f=initialize(t,k,l);
  global Lmax Lp sh n lma

  if k<=Lp
     m=Lmax(k)-1;
     if (k>1)
        m=min(m,part(l,k-1));
     end
     for i=1:m
        l=l+lma(k);
        sh(l+1,1:t,1:n)=-ones(t,n);
        initialize(t,k+1,l);
     end
  end

% given an integer pn that represents a partition, this function returns
% part i of the partition in O(1) time

function f=part(pn,i)
global lma

if i>length(lma)
   f=0;
else
   if i==1
      f=floor(pn/lma(i));
   else
      f=floor(mod(pn,lma(i-1))/lma(i));
   end
end


% nmu computes the unique integer that represents the partition l

function f=nmu(l)
global Lmax Lp
  f=0;
  for i=1:Lp
     f=Lmax(i)*f;
     if i<=length(l)
        f=f+l(i);
     end
  end


% Given vector x and partition, represented by an integer l, sch(m,j,k,l)
% computes the Schur function s_l(x_m,...,x_j), where the partitions
% mu, such that l-mu is horizontal strip, are generated to have 
% l(i)=mu(i) for i=1,2,...,k. The parameter k is designed to keep 
% track of the recursion. from outside the function should be called
% with sch(m,j,1,l)

function f=sch(m,j,k,l);
  global x Lp Lmax n sh lma

  s=1;

  if (m<=j)&(l>0)
    t=sh(l+1,m,j);
    if (k==1)&(t~=-1)
      s=t;
    elseif part(l,j-m+2)>0
      s=0;
    elseif part(l,j-m+1)>0
      t=part(l,j-m+1);
      s=s*prod(x(m:j))^t *sch(m,j,1,l-t*sum(lma(1:j-m+1)));  
    elseif j==m
        s=s*x(m)^l(1);
    else
      i=k;
      s=sch(m,j-1,1,l); 
      if (k==1)&(part(l,1)>part(l,2))
        i=2;
        s=s+sch(m,j,1,l-lma(1))*x(j);
      end
      while part(l,i)>0
        if part(l,i)>part(l,i+1)
          if part(l,i)>1
             s=s+x(j)*sch(m,j,i,l-lma(i));
          else
             s=s+x(j)*sch(m,j-1,1,l-lma(i));
          end
        end
        i=i+1;
      end
    end
    if k==1
      sh(l+1,m,j)=s;
    end
  end 
  f=s;
