#
# Symmetric Functions package: SF version 2.4, vanilla edition.
# This version/edition has been tested on
#   Maple V R3, R4, R5, Maple 6, Maple 7, Maple 9, and Maple 9.5.
#
# This is *not* a Maple worksheet.
#
# After loading this file during a Maple session, each function in the
# package can be accessed using the calling sequence
#
#   SF[<functionname>](<arguments>).
#
# In order to use package functions in the abbreviated form
#
#   <functionname>(<arguments>),
#
# run the command 'withSF()' after loading this file. If there is a
# conflict between the names of one of these functions and another name
# in the same session, a warning is printed.
#
# In order to use the abbreviated form for a subset of the procedures in
# the package, run the command
#
#   withSF(<name1>,<name2>,...).
#
# For an introduction, see
#   http://www.math.lsa.umich.edu/~jrs/software/SF.ps
# For documentation on the individual functions, see
#   http://www.math.lsa.umich.edu/~jrs/software/SFhelp.html
#
# Copyright (c) 2005 by John R. Stembridge
#
#########################################################################
#
SF:=table():
e:='e': p:='p': h:='h': s:='s': cl:='cl':
`SF/Bases`:={e,h,p,s[]}:
#
# assign short names, printing warnings if conflicts occur.
#
withSF:=proc() local install,f;
  install:=proc(x)
    if not assigned(SF[x]) then
      ERROR(cat(x,` is not a top level function in SF`))
    elif eval(x)<>eval(SF[x]) then
      if x='scalar' or x='conjugate' then unprotect(x) fi;
      if assigned(x) then printf(`Warning: new definition for  %a\n`,x) fi;
      assign(x,SF[x])
    fi; x
  end;
  if nargs>0 then map(install,[args]) else
    f:=proc() map(op,[args]) end; # hack the names w/o full evaluation!
    map(install,f(indices(SF)))
  fi
end:
#
#
# Define a new basis b that is orthogonal with respect to some inner
#  product ip and triangular w.r.t conjugate lex-ordered monomials.
# Supply an optional partition-eating procedure as args[3] to specify the
#  leading coefficient of each b[mu] (default = 1).
# If the final argument (args[3] or args[4]) is not a procedure, set a
#  flag indicating that this basis is triangular w.r.t dominance order.
#
# Note: `SF/iprod`, `SF/lcoeff` and `SF/nondom` are subversively global.
#
`SF/add_basis`:=proc(b,ip) local a,n,flag,lc;
  if type(b,'indexed') then a:=op(0,b) else a:=b fi;
  if member(a[],`SF/Bases`) or member(a,`SF/Bases`) then
    ERROR(cat(`base `,a,` is already in use`)) fi;
  n:=nargs; flag:=type(args[n],'procedure');
  if not flag then n:=n-1 fi;
  if n>2 then lc:=args[3] else lc:=1 fi;
  assign([evaln(`SF/Bases`)={op(`SF/Bases`), a[]},
    `SF/iprod`[a]=ip, `SF/lcoeff`[a]=lc, `SF/nondom`[a]=flag]);
  assign(cat(`to`,b)=proc() local c,f,d,B,i,v,vars,kill,nr,res,mu,nu;
    B:=substring(procname,3..length(procname));
    f:=SF['toe'](args); res:=0;
    nr:=`SF/getrows`(args);
    d:=SF['stdeg'](f,'e');
    if nr=0 then nr:=d else nr:=min(nr,d) fi;
    kill:=[seq(cat('e',i)=0,i=nr+1..d)];
    vars:=[seq(cat('e',-i),i=-nr..-1)];
    while f<>0 do
      c:=tcoeff(f,vars,'v');
      nu:=`SF/shape`(v,'e',nr); mu:=SF['conjugate'](nu);
      c:=normal(c/`SF/lcoeff`[B](mu));
      res:=c*B[op(mu)]+res;
      f:=f-c*subs(kill,`SF/added`(B,nu));
      f:=`SF/normal`(f,vars,v);
    od; res
  end);
  'Okay';
end:
#
# Assume b is a basis created by 'add_basis', and that sp is either
#  a list of partitions or a single partition.
# With three arguments, replace all occurrences of b[mu] in f0 for
#  all mu in sp with their corresponding e-expansions.
# With two arguments, assume sp is a single partition,
#  and return the e-expansion of b[conjugate(sp)].
#
`SF/added`:=proc(b,sp,f0) local f,mu;
  if nargs>2 then
    f:=subs({seq(b[op(mu)]=
      `SF/added/gs`(b,SF['conjugate'](mu))[2],mu=sp)},f0)
  else
    f:=`SF/added/gs`(b,sp)[2]
  fi;
  subs(map(x->`SF/added/etable`(x)[1],SF['varset'](f,e[])),f)
end:
#
# Given a user-added basis b, compute and remember the e-monomial
#  expansion of b[mu'] (NOTICE: this is conjugate(mu), not mu)
#  via the Gram-Schmidt algorithm. The initial ordered basis is
#  the set of e-monomials, ordered lexicographically as in 'Par'.
# The results are returned in the form [N,f], where N is the squared norm
#  of b[mu'] and f is an expression for b[mu'] as a collected linear
#  combination of terms e[...], where (for example) e[3,3,2,1] represents
#  the e-monomial e3^2*e2*e1. 
# If `SF/nondom`[b] is false, assume that basis b
#  is triangular with respect to dominance order.
# FUTURE: The optimal frequency of collecting/normalizing the coeffs may
#  need further tuning depending on the ground ring and the Maple version.
#
`SF/added/gs`:=proc(b,mu) local g,lc,f,f0,sc,nu,N,c;
  option remember;
  lc:=`SF/lcoeff`[b](SF['conjugate'](mu));
  f:=lc*e[op(mu)]; g:=`SF/added/etable`(mu)[2];
  sc:=[]; nu:=subs(0=NULL,[convert(mu,`+`)]);
  N:=lc^2*SF['scalar'](g,g,0,0,`SF/iprod`[b]);
  while nu<>mu do
    if `SF/nondom`[b] or SF['dominate'](nu,mu) then
      c:=`SF/added/etable`(nu)[2];
      c:=normal(SF['scalar'](g,c,0,0,`SF/iprod`[b]));
      sc:=[op(sc),e[op(nu)]=c];
      f0:=`SF/added/gs`(b,nu);
      c:=normal((lc/f0[1])*subs(sc,f0[2])); f:=f-c*f0[2];
      if length(f)>20000 then  #see FUTURE
        f:=collect(f,map(x->op(1,x),sc),'distributed',normal) fi;
      N:=N-c^2*f0[1];
      if length(N)>1000 then N:=normal(N) fi;  #see FUTURE
    fi;
    nu:=SF['nextPar'](nu)
  od;
  [factor(N),collect(f,map(x->op(1,x),sc),'distributed',normal)];
end:
#
# Remember two things about each e-monomial:
# 1. A substitution that replaces the dummy variable e[mu] used by
#   SF/added/gs with the actual e-monomial, and
# 2. The p-expansion of the e-monomial in the format used by 'scalar'.
#
`SF/added/etable`:=proc(mu) local v,g,tf,vars,j,d;
  option remember;
  v:=convert([seq(cat('e',j),j=mu)],`*`);
  d:=op(1,[op(mu),0]);
  vars:=[seq(cat('p',j),j=1..d)];
  g:=coeffs(SF['top'](v,'e'),vars,'tf');
  [e[op(mu)]=v,[[g],[tf],d]]
end:
#
# char2sf(chi) applies the characteristic map to the class function chi,
# producing a symmetric function as the output. The class funtion must be
# expressed as a linear combination of characteristic functions cl[mu]
# for various partitions mu. By default, the output will be a p-polynomial.
# Use char2sf(chi,b) to require output in base b.
#
`SF/char2sf`:=proc(chi,b) local res,mu;
  res:=SF['varset'](chi,'cl[]');
  res:=convert([seq(coeff(chi,cl[op(mu)])*
    convert(map(x->cat('p',x),mu),`*`)/SF['zee'](mu),mu=res)],`+`);
  if nargs>1 then `SF/apply`(b,res,'p') else res fi;
end: 
#
# conjugate(mu) returns the conjugate of partition mu.
# mu must be a list in decreasing order.
#
`SF/conjugate`:=proc(mu) local i,l;
  l:=nops(mu);
  if l=0 then [] else
    [l$mu[l],seq((-i)$(mu[-i]-mu[1-i]),i=1-l..-1)]
  fi
end:
#
# dominate(mu)    list all partitions <= mu in dominance order
# dominate(mu,n)  does the same but only for partitions with <= n rows
# dominate(mu,nu) returns true iff mu >= nu in dominance order
#
`SF/dominate`:=proc(mu) local n,nu,res,sat,i,j,m,nu0,lam;
  n:=convert(mu,`+`);
  if nargs>1 then
    if type(args[2],'list') then
      nu:=zip((x,y)->x-y,mu,args[2],0); m:=0;
      for i in nu do m:=m+i; if m<0 then RETURN(false) fi od;
      RETURN(true)
    else n:=min(n,args[2]) fi
  fi;
  if n<nops(mu) then res:=[] else res:=[mu] fi;
  for sat while sat<=nops(res) do
    nu:=res[sat]; m:=nops(nu); nu0:=[op(nu),0];
    for i to m do 
      if nu0[i]=nu0[i+1] then next
        elif nu0[i]>nu0[i+1]+1 then j:=i+1
        else for j from i+2 to m while nu0[j-1]=nu0[j] do od;
      fi;
      if j>min(n,m+1) then next
        elif j<=m then lam:=subsop(i=nu[i]-1,j=nu[j]+1,nu)
        else lam:=[op(subsop(i=nu[i]-1,nu)),1]
      fi;
      if not member(lam,res) then res:=[op(res),lam] fi;
    od
  od; res
end:
#
# Define a new basis b1 that is dual to some existing basis b2.
# Supply an optional scalar product as args[3] (default = zee).
#
# `SF/Bases`, `SF/dual`, and `SF/iprod` are subversively global
#
`SF/dual_basis`:=proc(b1,b2) local a,ip;
  if type(b1,'indexed') then a:=op(0,b1) else a:=b1 fi;
  if member(a[],`SF/Bases`) or member(a,`SF/Bases`) then
    ERROR(cat(`base `,a,` is already in use`)) fi;
  if nargs>2 then ip:=args[3] else ip:=SF['zee'] fi;
  assign([`SF/dual`[a]=`SF/verify`(b2), `SF/iprod`[a]=ip,
    evaln(`SF/Bases`)={op(`SF/Bases`),a[]}]);
  assign(cat(`to`,a)=proc() local f,B,b;
    B:=substring(procname,3..length(procname)); 
    b:=`SF/dual`[B]; f:=SF['top'](args);
    if type(args[nargs],'list') then
      `SF/dual_basis/to`([f,0],B,b,args[nargs])
    else
      f:=`SF/homog_cmps`(f,'p');
      convert(map(`SF/dual_basis/to`,f,B,b),`+`)
    fi
  end);
  'Okay';
end: 
#
# Convert a p-polynomial f0=[f,d] into a linear combination of B[mu]'s.
# If there is a fourth argument, it should be a list of partitions
# that contain the B-support of f. Otherwise, we make no assumptions
# about B-support except that f is homogeneous of degree d.
#
`SF/dual_basis/to`:=proc(f0,B,b) local f,sp,res,mu,v,j,d;
  if nargs>3 then sp:=args[4] else sp:=SF['Par'](f0[2]) fi;
  d:=SF['varset'](f0[1],'p');
  f:=[coeffs(f0[1],[seq(cat('p',j),j=1..d)],'v')];
  f:=[f,[v],d]; res:=0;
  for mu in sp do
    if type(b,'indexed') then v:=op(0,b)[op(mu)]
      else v:=convert(map((x,y)->cat(y,x),mu,b),`*`) fi;
    v:=SF['scalar'](f,v,0,b,`SF/iprod`[B]);
    res:=res+normal(v)*B[op(mu)];
  od; res
end:
#
# Assume 'a' is a basis defined via dual_basis.
# For each partition mu listed in sp, replace all instances of a[mu]
# in the symmetric function f with the corresponding p-polynomial.
#
# This requires an expensive computation for each a-variable, so we
# economize by doing this only once for each homogeneous degree that
# occurs among the a-variables indexed by sp.
#
`SF/dualize`:=proc(a,sp) local cmps,v,vars,t,f,g,b,pr;
  b:=`SF/dual`[a]; f:=args[3];
  if nops(sp)=0 then RETURN(f) fi;
  cmps:=convert([seq(a[op(v)]*t^convert(v,`+`),v=sp)],`+`);
  cmps:=map(indets,[coeffs(cmps,t)]);
  for vars in cmps do
    if degree(f,vars)=1 then
      f:=collect(f,vars);  # required only for Maple V.3
      pr:=(x,y)->normal(x/y);
      f:=`SF/dualize/hg`(f,vars,a,b,pr)
    else
      g:=convert([seq(v*t[op(v)],v=vars)],`+`);
      pr:=proc(x,y) if type(x,`+`) then
        map((X,Y)->normal(X/Y),x,y) else normal(x/y) fi end;
      g:=`SF/dualize/hg`(g,vars,a,b,pr);
      g:=collect(g,[seq(t[op(v)],v=vars)]);  # required for Maple V.3
      f:=subs({seq(v=coeff(g,t[op(v)]),v=vars)},f);
    fi
  od; f
end:
#
# Given a non-empty set of homogeneous a-variables, the dual basis b for a,
# and a symmetric function f that is *linear* in these a-variables,
# substitute p-polynomials for the occurrences of these a-variables in f.
# Future: might be worth it to remember-table the bv's.
#
`SF/dualize/hg`:=proc(f,vars,a,b,pr) local mu,g,v,sp,i,bv,scp,d;
  d:=convert(vars[1],`+`);
  if type(b,'indexed') then
    sp:=map((x,y)->[x,y[op(x)]],vars,op(0,b));
    scp:=`SF/dualize/ind`;
  else
    sp:=[seq([v,convert(map((x,y)->cat(y,x),[op(v)],b),`*`)],v=vars)];
    sp:=sp,[seq(cat(b,i),i=1..d)];
    scp:=`SF/dualize/mult`;
  fi;
  mu:=subs(0=NULL,[d]);
  g:=subs({seq(v=0,v=vars)},f);
  while mu<>NULL do
    v:=convert(map(x->cat('p',x),mu),`*`);
    bv:=`SF/apply`(b,v,'p');
    g:=g+pr(scp(f,bv,sp),`SF/iprod`[a](mu))*v;
    mu:=SF['nextPar'](mu)
  od; g
end:
#
# Compute the scalar product of f with the b-expression bv, given that
# the a-variables in f and b-variables in bv are listed by pairs in sp.
#
# Either b is a basis name of indexed type,...
#
`SF/dualize/ind`:=proc(f,bv,sp) local u;
  convert([seq(coeff(f,u[1])*coeff(bv,u[2]),u=sp)],`+`)
end:
#
# ...or b is a multiplicative basis name.
#
`SF/dualize/mult`:=proc(f,bv,sp,bvars) local u,g,tms,res,j;
  g:=[coeffs(bv,bvars,'tms')];
  tms:=[tms]; res:=0;
  for u in sp do
    if member(u[2],tms,'j') then res:=res+coeff(f,u[1])*g[j] fi;
  od; res
end:
#
# evalsf(f,a) generates a plethystic evaluation of f at a.
# Let a<j> denote the result of substituting x=x^j and p.i=p.(i*j) for
# each variable x and each power-sum p.i in the p-expansion of a.
# Then evalsf(f,a) is obtained by substituting p.j=a<j> in the power
# sum expansion of f for j=1,2,3,...
# The output is collected and normalized.
#
`SF/evalsf`:=proc() local f,df,a,da,expr,j,i;
  f:=SF['top'](args[1]); a:=SF['top'](args[2]);
  df:=SF['varset'](f,'p'); da:=SF['varset'](a,'p');
  expr:=indets(a) minus {seq(cat('p',i),i=1..da)};
  expr:=subs({seq(i=i^j,i=expr),seq(cat('p',i)='cat'('p',i*j),i=1..da)},a);
  expr:=subs({seq(cat('p',j)=eval(expr),j=1..df)},f);
  collect(expr,[seq(cat('p',i),i=1..da*df)],'distributed',normal)
end:
#
# homog_cmps(f) separates the symmetric function f into its homogeneous
#  components and returns a list [[f1,d1],[f2,d2],...], where f1,f2, are
#  the components and d1,d2,... are their degrees (unsorted).
# homog_cmps(f,b) does the same, assuming f is in base b.
#
`SF/homog_cmps`:=proc() local t,sp,i,f,x,b,mu,bases;
  if nargs>1 then bases:={args[2]} else bases:=`SF/Bases` fi;
  f:=args[1]; sp:=SF['varset'](f,bases);
  for x in bases do
    if type(x,'indexed') then
      b:=op(0,x);
      f:=subs({seq(b[op(mu)]=t^convert(mu,`+`)*b[op(mu)],mu=sp[b])},f)
    else
      f:=subs({seq(cat(x,i)=t^i*cat(x,i),i=1..sp[x])},f)
    fi
  od;
  f:=[coeffs(collect(f,t),t,'sp')];
  zip((x,y)->[x,degree(y)],f,[sp])
end:
#
# hooks(lambda)      is the list of hook-lengths in lambda.
# hooks(lambda,a)    is the Jack hook-product.
# hooks(lambda,q,t)  is the two-variable hook polynomial.
#
`SF/hooks`:=proc(mu,z,w) local nu,i,j,q,t;
  nu:=SF['conjugate'](mu);
  if nargs=1 then
    nu:=[seq(seq(mu[i]-i+nu[j]-j+1,j=1..mu[i]),i=1..nops(mu))];
    sort(nu,(x,y)->evalb(x>y))
  elif nargs=2 then
    nu:=[seq(seq(z*(mu[i]-j)+nu[j]-i+1,j=1..mu[i]),i=1..nops(mu))];
    convert(nu,`*`)
  else
    nu:=[seq(seq((1-q^(mu[i]-j)*t^(nu[j]-i+1),j=1..mu[i])),i=1..nops(mu))];
    subs(q=z,t=w,convert(nu,`*`))   # Maple is immature about 0^0
  fi
end:
#
# itensor(f,g) computes the inner tensor product (a.k.a. the internal
#  or Kronecker product)of symmetric functions f and g.
# itensor(f,g,b1,b2) does the same assuming f is in base b1 and g in b2.
# itensor(f,g,b) does the same, but converts the output to base b.
#  (The default base is 'p').
# itensor(f,g,b1,b2,b) combines both options.
#
`SF/itensor`:=proc() local b,d,j,f,g,tf,u,c,res,vars;
  if nargs>3 then b:=[args[3..4]] else b[1]:=NULL; b[2]:=NULL fi;
  f:=[seq(SF['top'](args[j],b[j]),j=1..2)];
  d:=SF['varset'](f,'p'); res:=0;
  vars:=[seq(cat('p',j),j=1..d)];
  f:=[[coeffs(f[1],vars,'tf')],[tf],[coeffs(f[2],vars,'tf')],[tf]];
  if nops(f[2])<=nops(f[4]) then
    g:=zip((x,y)->[x,y],f[1],f[2]); tf:=f[4]; f:=f[3]
  else
    g:=zip((x,y)->[x,y],f[3],f[4]); tf:=f[2]; f:=f[1]
  fi;
  for u in g do
    if member(u[2],tf,'j') then
      c:=SF['zee'](`SF/shape`(u[2],'p',d));
      res:=res+normal(c*u[1]*f[j])*u[2]
    fi;
  od;
  if modp(nargs,2)=1 then `SF/apply`(args[nargs],res,'p') else res fi;
end:
#
# When nrows=0, evaluate s[mu] as the determinant of a Jacobi-Trudi
#  matrix of either e's or h's, whichever version has fewer rows.
# Supply an optional third argument b='e' or b='h' to override the choice.
# When nrows>0, then all of the matrix variables b.n for n>nrows are
#  killed before the determinant is computed.
#
`SF/jt_det`:=proc(mu,nrows) local b,nu;
  if mu=[] then RETURN(1) fi;
  if nargs>2 then b:=args[3]
    elif nops(mu)>mu[1] then b:='e'
    else b:='h'
  fi;
  if b='e' then nu:=SF['conjugate'](mu) else nu:=mu fi;
  if nrows>0 then b:=b,nrows fi;
  expand(linalg['det'](SF['jt_matrix'](nu,[],b)));
end:  
#
# jt_matrix(lambda) produces the Jacobi-Trudi matrix corresponding to
#   the partition lambda.
# jt_matrix(lambda,mu) does the same, but for the skew shape defined by
#   the pair of partitions lambda,mu (mu = inner shape).
# jt_matrix(lambda,mu,b) does the same, but uses b as the base name for
#   the entries of the matrix.
# jt_matrix(lambda,mu,b,m) does the same, but kills all of the matrix
#   entries of the form b.k for k>m.
# Note that if lambda=mu=[], then the result should be a 0x0 matrix.
# But Maple (wrongly) disallows this, so we produce a 1x1 matrix.
#
`SF/jt_matrix`:=proc() local b,n,i,j,mu,nu;
  mu:=args[1];
  if nargs>1 then nu:=args[2] else nu:=[] fi;
  if nargs>2 then b:=args[3..nargs] else b:='h' fi;
  n:=max(1,nops(mu),nops(nu));
  mu:=[op(mu),0$n]; nu:=[op(nu),0$n];
  array([seq([seq(`SF/jt_matrix/ent`(mu[i]-i+j-nu[j],b),j=1..n)],i=1..n)])
end:
#
`SF/jt_matrix/ent`:=proc(k,b)
  if k<0 then 0
    elif k=0 then 1
    elif nargs>2 and k>args[3] then 0
    else cat(b,k)
  fi
end:
#
# Generate the next partition after mu in lexicographic order.
# If mu is the last partition (i.e., mu=[1,1,...,1]), return NULL.
# Note that this generates partitions in the same order as Par(n),
#  but is noticeably faster and more space-efficient for larger n.
#
# Example: mu:=[9]; while mu<>NULL do <something>; mu:=nextPar(mu) od;
#
`SF/nextPar`:=proc(mu) local i,k,m,r;
  if member(1,mu,'i') then i:=i-1 else i:=nops(mu) fi;
  if i=0 then NULL else
    k:=mu[i]-1; m:=iquo(nops(mu)-i+mu[i],k,'r');
    if r=0 then r:=NULL fi;
    [op(1..i-1,mu),k$m,r]
  fi
end:
#
# omega(f)       Apply the omega-automorphism to symmetric function f.
# omega(f,b1)    The same, but assume f is expressed in terms of basis b1.
# omega(f,b1,b2) The same, but convert the output to basis b2.
#
# With no output basis specified, the result may involve any convenient
# mix of functions from various bases.
#
`SF/omega`:=proc() local f,sp,j,b,a,b0,mu,b1,bases;
  f:=args[1]; b0:={'p','h','e','s[]'}; bases:=`SF/Bases`;
  if nargs>1 then b1:=`SF/verify`(args[2]); bases:={b1} fi;
  bases:=bases minus b0;
  sp:=SF['varset'](f,bases);
  for b in bases do
    a:=op(0,b); if sp[a]=[] then next fi;
    if assigned(`SF/dual`[a]) then
      f:=`SF/dualize`(a,sp[a],f); b1:='p'
    else
      f:=`SF/added`(a,sp[a],f); b1:='e'
    fi;
    if nargs>1 then f:=`SF/apply`(b1,f,b1) fi;
  od;
  b1:=subs({'e'='h','h'='e'},b1);
  sp:=SF['varset'](f,b0);
  f:=subs({seq(cat('e',j)=cat('h',j),j=1..sp['e']),
    seq(cat('h',j)=cat('e',j),j=1..sp['h']),
    seq(cat('p',j)=(-1)^(j-1)*cat('p',j),j=1..sp['p']),
    seq(s[op(mu)]=s[op(SF['conjugate'](mu))],mu=sp['s'])},f);
  if nargs>2 then `SF/apply`(args[3],f,b1) else f fi;
end:
#
# Par(n) returns a list of all partitions of n.
# Par(n,l) returns the partitions of n with length <=l.
# Par(n,k,l) returns the partitions of n with parts <=k, length <=l.
#
# Note: partition lists are lex-sorted first by largest part,
# then second-largest part, etc. This consistency is important for
# basis-conversion routines that involve Gram-Schmidt.
#
`SF/Par`:=proc(n)
  if nargs=1 then `SF/Par/sub`(n,n,n)
    elif nargs=2 then `SF/Par/sub`(n,n,args[2])
    else `SF/Par/sub`(args)
  fi
end:
#
`SF/Par/sub`:=proc(n,row,col) local i;
  if n=0 then [[]]
    elif col=0 then []
    else 
      [seq(op(map((x,y)->[y,op(x)],`SF/Par/sub`(n+i,-i,col-1),-i)),
        i=-min(row,n)..-iquo(n+col-1,col))]
  fi
end:
#
# plethysm(f,g) computes the plethysm f[g] of symmetric functions f,g.
# plethysm(f,g,b1,b2)   does the same assuming f is in base b1 and g in b2.
# plethysm(f,g,b)       does the same but converts the output to base b.
# plethysm(f,g,b1,b2,b) does both. The default output is base 'p'.
#
`SF/plethysm`:=proc() local d,i,f,b;
  if nargs>3 then b:=[args[3..4]] else b[1]:=NULL; b[2]:=NULL fi;
  f:=[seq(SF['top'](args[i],b[i]),i=1..2)];
  d:=map(SF['varset'],f,'p');
  f:=subs(map(`SF/plethysm/one`,[$1..d[1]],d[2],f[2]),f[1]);
  f:=collect(f,[seq(cat('p',i),i=1..d[1]*d[2])],'distributed',normal);
  if modp(nargs,2)=1 then `SF/apply`(args[nargs],f,'p') else f fi;
end:
#
# Compute the special plethysm f[g], where f = p.r and g has degree dg,
# and return the result as a substitution for f.
#
`SF/plethysm/one`:=proc(r,dg,g) local i;
  cat('p',r)=subs({seq(cat('p',i)=cat('p',r*i),i=1..dg)},g)
end:
#
# scalar(f,g) computes the scalar product of the symmetric functions f
#  and g with respect to the form for which the power sums are orthogonal
#  and <p[mu],p[mu]>=zee(mu).
# scalar(f,g,b1,b2) does the same, assuming f is in base b1 and g in b2.
# scalar(f,g,Z) or scalar(f,g,b1,b2,Z) does the same, relative to the
#  form for which the power sums are orthogonal and <p[mu],p[mu]>=Z(mu).
#
`SF/scalar`:=proc() local f,g,tf,ip,b,d,u,mu,j,res;
  if nargs=3 or nargs=5 then ip:=args[nargs] else ip:=SF['zee'] fi;
  if nargs>3 then b:=[args[3..4]] else b[1]:=NULL; b[2]:=NULL fi;
  f:=[seq(`SF/scalar/parse`(args[j],b[j]),j=1..2)];
  res:=0; d:=min(f[1][3],f[2][3]);
  if nops(f[1][2])<=nops(f[2][2]) then
    g:=zip((x,y)->[x,y],f[1][1],f[1][2]); tf:=f[2][2]; f:=f[2][1]
  else
    g:=zip((x,y)->[x,y],f[2][1],f[2][2]); tf:=f[1][2]; f:=f[1][1]
  fi;
  for u in g do
    if member(u[2],tf,'j') then
      mu:=`SF/shape`(u[2],'p',d);
      res:=res+ip(mu)*u[1]*f[j];
    fi
  od; res
end:
#
# Parse the data specified by f and a base/flag b:
# If b=0 assume that f=[coeffs,terms,d], represent the coefficients,
# terms, and max occurring p.d in some p-polynomial. Return f.
# Otherwise, assume f is a symmetric function in base b (or in unknown
# form if b=NULL). Convert f to base p and return [coeffs,terms,d].
#
`SF/scalar/parse`:=proc(f,b) local g,i,d,tms;
  if nargs>1 and b=0 then RETURN(f) fi;
  g:=SF['top'](args); d:=SF['varset'](g,'p');
  [[coeffs(g,[seq(cat('p',i),i=1..d)],'tms')],[tms],d]
end:
#
# sf2char(f) applies the inverse characteristic map to sym fun f.
# The result is expressed as a linear combination of characteristic
# functions cl[mu] for partitions mu.
# sf2char(f,b) does the same, assuming that f is in base b.
#
`SF/sf2char`:=proc() local f,i,d,res,cfs,term,mu;
  f:=SF['top'](args); d:=SF['varset'](f,'p');
  cfs:=[coeffs(f,[seq(cat('p',i),i=1..d)],'term')];
  term:=[term]; res:=0;
  for i to nops(cfs) do
    mu:=`SF/shape`(term[i],'p',d);
    res:=res+SF['zee'](mu)*cfs[i]*cl[op(mu)];
  od; res
end:
#
# skew(f,g)       applies the linear transformation f^* to g,
#  where f^* denotes the adjoint to multiplication by f.
# skew(f,g,ip)    does the same, relative to the inner product ip.
# skew(f,g,b1,b2) does the same, assuming f is in base b1 and g is in b2.
# skew(f,g,b1,b2,ip) combines both options.
# 
`SF/skew`:=proc() local f,ip,b,vars,i,g;
  if modp(nargs,2)=0 then ip:=SF['zee'] else ip:=args[nargs] fi;
  if nargs>3 then b:=[args[3..4]] else b[1]:=NULL; b[2]:=NULL fi;
  f:=[seq(SF['top'](args[i],b[i]),i=1..2)];
  f:=map((x,y)->map(`SF/skew/parse`,`SF/homog_cmps`(x,y)),f,'p');
  f:=convert([seq(op(map(`SF/skew/hg`,f[1],g,ip)),g=f[2])],`+`);
  vars:=[seq(cat('p',i),i=1..SF['varset'](f,'p'))];
  collect(f,vars,'distributed',normal);
end:
#
# Given f=[poly,d], where poly is a homogeneous p-polynomial
# of degree d, return the triple [coeffs,terms,d].
#
`SF/skew/parse`:=proc(f) local i,tf;
  [[coeffs(f[1],[seq(cat('p',i),i=1..f[2])],'tf')],[tf],f[2]]
end:
#
# Apply f^* to g, assuming f and g are homogeneous p-polynomials
# structured by the output of skew/parse.
#
`SF/skew/hg`:=proc(f,g,ip) local d,c,mu,v,res;
  d:=g[3]-f[3]; if d<0 then RETURN(0) fi;
  mu:=subs(0=NULL,[d]); res:=0;
  while mu<>NULL do
    v:=convert(map(x->cat('p',x),mu),`*`);
    c:=[f[1],map((x,y)->x*y,f[2],v),g[3]];
    c:=SF['scalar'](g,c,0,0,ip);
    res:=res+normal(c/ip(mu))*v;
    mu:=SF['nextPar'](mu);
  od; res
end:
#
# stdeg(f)   determine the degree of f w.r.t. the standard grading.
# stdeg(f,b) do the same, but assume f is in base b.
#
`SF/stdeg`:=proc() local f,bases,b,B,i,sp,t;
  if nargs>1 then bases:={args[2]} else bases:=`SF/Bases` fi;
  f:=args[1]; sp:=SF['varset'](f,bases);
  for b in bases do
    if type(b,'indexed') then B:=op(0,b);
      f:=subs({seq(B[op(i)]=t^convert(i,`+`)*B[op(i)],i=sp[B])},f)
    else
      f:=subs({seq(cat(b,i)=t^i*cat(b,i),i=1..sp[b])},f)
    fi
  od;
  degree(f,t)
end:
#
# subPar(mu)   returns all partitions that fit inside the diagram of mu.
# subPar(mu,n) does the same, restricted to the set of partitions of n.
#
`SF/subPar`:=proc(mu,n) local m,i,l,nu,j;
  l:=nops(mu);
  if nargs=1 then
    if mu=[] then RETURN([[]]) fi;
    for i to l-1 while mu[i]=mu[i+1] do od;
    if mu[i]>1 then j:=mu[i]-1 else j:=NULL fi;
    [seq([mu[1]$i,op(nu)],nu=`SF/subPar`([op(i+1..l,mu)])),
      op(`SF/subPar`(subsop(i=j,mu)))]
  else
    m:=convert(mu,`+`);
    if n>m or n<0 then RETURN([]) elif n=m then RETURN([mu]) fi;
    for i to l-1 while mu[i]=mu[i+1] do od; 
    if mu[i]>1 then j:=mu[i]-1 else j:=NULL fi; 
    [seq([mu[1]$i,op(nu)],nu=`SF/subPar`([op(i+1..l,mu)],n-i*mu[1])),
      op(`SF/subPar`(subsop(i=j,mu),n))]
  fi
end:
#
# theta(f,a)   Apply the automorphism p.j -> a*p.j to f.
# theta(f,q,t) Apply the automorphism p.j -> (1-q^j)/(1-t^j)*p.j to f.
#
`SF/theta`:=proc(g,q,t) local f,d,j,vars;
  f:=SF['top'](g); d:=SF['varset'](f,'p');
  vars:=[seq(cat('p',j),j=1..d)];
  if nargs=2 then
    f:=subs({seq(vars[j]=q*vars[j],j=1..d)},f)
  elif nargs=3 then
    f:=subs({seq(vars[j]=(1-q^j)/(1-t^j)*vars[j],j=1..d)},f)
  else
    ERROR(`wrong number of parameters`)
  fi;
  collect(f,vars,'distributed',normal);
end:
#
# toe(f)    convert the symmetric function f into an e-polynomial.
# toe(f,b)  does the same, assuming f is expressed solely in basis b.
# The final result is collected with respect to e1,e2,e3,...
# If b is not a predefined basis, f must be *linear* in the b[..]'s.
#
`SF/toe`:=proc() local f,bases,sp,i,mu,d,b,c,nrows;
  nrows:=`SF/getrows`(args); f:=args[1];
  bases:=`SF/getbase`(args) minus {'e'};
  sp:=SF['varset'](f,bases minus {'p','h'});
  for b in bases minus {'h','p','s[]'} do
    c:=op(0,b); if sp[c]=[] then next fi;
    if assigned(`SF/dual`[c]) then
      f:=`SF/dualize`(c,sp[c],f); bases:={op(bases),'p'}
    else
      f:=`SF/added`(c,sp[c],f)
    fi
  od;
  if member('s[]',bases) then b:='e';
    if nops(sp['s'])>1 and nrows=0 then b:=NULL; bases:={op(bases),'h'} fi;
    f:=subs({seq(s[op(mu)]=`SF/jt_det`(mu,nrows,b),mu=sp['s'])},f);
  fi;
  if nrows>0 then
    d:=SF['varset'](f,'e');
    f:=subs({seq(cat('e',i)=0,i=nrows+1..d)},f)
  fi;
  if member('p',bases) then f:=`SF/to_ehp`('p','e',f,nrows) fi;
  if member('h',bases) then f:=`SF/to_ehp`('h','e',f,nrows) fi;
  d:=SF['varset'](f,'e');
  collect(f,[seq(cat('e',i),i=1..d)],'distributed',normal);
end: 
#
# Assume that b1 and b2 are (distinct) members of {e,h,p}.
# to_ehp(b1,b2,f) takes all occurrences of the variables b1.i,i=1,2,... in
#  f and substitutes equivalent expressions in the variables b2.i,i=1,2,...
# Redesigned for v2.4: we now use convolution identities. Seems to be
#  faster and more space-efficient for high degrees.
# Also, it is *way* faster and slightly more space efficient to expand
#  one variable at a time when the ground field is Q, but better to
#  collect the whole at once over a rational function field.
#
`SF/to_ehp`:=proc(b1,b2) local i,f,d,g,sp,vars;
  f:=args[3]; d:=SF['varset'](f,{b1,b2});
  if d[b1]=0 then RETURN(f) fi;
  g:=`SF/to_ehp/subs`(b1,b2,d[b1],args[4..nargs]);
  sp:=[seq(cat(b1,i),i=1..d[b1])];
  vars:=[op(sp),seq(cat(b2,i),i=1..d[b2])];
  if type(f,'polynom'('rational',vars)) then
    vars:=indets(f) intersect {op(sp)};
    for i to d[b1] do
      if member(sp[i],vars) then f:=expand(subs(sp[i]=g[i],f)) fi
    od; f
  else
    subs({seq(sp[i]=g[i],i=1..d[b1])},f)
  fi
end:
#
# Create an array whose i-th term is the b2-expansion of b1.i, i=1..d.
#
`SF/to_ehp/subs`:=proc(b1,b2,d) local dr,i,k,vars,g,sig;
  if nargs>3 and args[4]>0 then dr:=min(d,args[4]) else dr:=d fi;
  if b1='e' or b2='e' then
    vars:=[seq((-1)^(i-1)*cat(b2,i),i=1..dr),0$(d-dr)]
  else
    vars:=[seq(cat(b2,i),i=1..dr),0$(d-dr)]
  fi;
  g:=array(0..d,[0=1]);
  if b2='p' then
    for k to d do
      g[k]:=expand(convert([seq(g[k-i]*vars[i],i=1..k)],`+`))/k
    od;
  elif b1='p' then
    if b2='h' then vars:=map(x->-x,vars); sig:=-1 else sig:=1 fi;
    for k to d do g[0]:=sig*k;
      g[k]:=expand(convert([seq(g[k-i]*vars[i],i=1..k)],`+`))
    od;
  else
    for k to d do
      g[k]:=expand(convert([seq(g[k-i]*vars[i],i=1..k)],`+`))
    od;
  fi; op(g) 
end:
#
# toh(f)    convert the symmetric function f into an h-polynomial.
# toh(f,b)  does the same, assuming f is expressed solely in basis b.
# The final result is collected with respect to h1,h2,h3,...
# If b is not a predefined basis, f must be *linear* in the b[..]'s.
#
`SF/toh`:=proc() local f,bases,sp,i,mu,d,b,c;
  f:=args[1]; bases:=`SF/getbase`(args) minus {'h'};
  sp:=SF['varset'](f,bases minus {'p','e'});
  for b in bases minus {'e','p','s[]'} do
    c:=op(0,b); if sp[c]=[] then next fi;
    if assigned(`SF/dual`[c]) then
      f:=`SF/dualize`(c,sp[c],f); bases:={op(bases),'p'}
    else
      f:=`SF/added`(c,sp[c],f); bases:={op(bases),'e'}
    fi
  od;
  if member('s[]',bases) then b:='h';
    if nops(sp['s'])>1 then b:=NULL; bases:={op(bases),'e'} fi;
    f:=subs({seq(s[op(mu)]=`SF/jt_det`(mu,0,b),mu=sp['s'])},f);
  fi;
  if member('p',bases) then f:=`SF/to_ehp`('p','h',f) fi;
  if member('e',bases) then f:=`SF/to_ehp`('e','h',f) fi;
  d:=SF['varset'](f,'h');
  collect(f,[seq(cat('h',i),i=1..d)],'distributed',normal);
end: 
#
## Various internal tools that are needed in multiple places
#
# Parse the argument list for an equation of the form 'nrows=<something>'.
# If found, return the <something> of the first such equation.
# Otherwise, return 0 as the default.
#
`SF/getrows`:=proc() local u;
  for u in [args] do
    if type(u,`=`) and op(1,u)='nrows' then RETURN(op(2,u)) fi
  od; 0
end:
#
# Parse the argument list (except args[1]) for a name.
# Verify the first one found and return it as a singleton set.
# If no names found, return `SF/bases` as the default.
#
`SF/getbase`:=proc() local i;
  for i from 2 to nargs do
    if type(args[i],'name') then RETURN({`SF/verify`(args[i])}) fi
  od; `SF/Bases`
end:
#
# Get the shape of a b-monomial f, given that b.d is the highest
# degree variable that might occur in f.
#
`SF/shape`:=proc(f,b,d) local k;
  [seq((-k)$degree(f,cat(b,-k)),k=-d..-1)]
end:
#
# Report the number of bytes allocated.
# We need this to work around feature-breakage that Maple imposes.
#
if `+`(0)=0 then # we are using Maple V.4 or later
  `SF/bytes`:=proc() kernelopts(bytesalloc) end
else
  `SF/bytes`:=proc() 4*status[2] end
fi:
#
# Put a sym function f into a normal form optimized for the ground field:
#  *If the coeffs are rational, do nothing.
#  *If the coeffs are polynomials over Q, expand f. This could waste a
# lot of memory, but the improvement in speed is spectacular.
#  *Otherwise, collect terms and apply normal to the coefficients.
#  *If the coefficients are not ratpolys, then floating-point arithmetic,
# radicals, or other obscenities are involved. In that case, make sure
# that the leading term v has been zapped (prevent infinite loops).
#
`SF/normal`:=proc(f,vars,v) local g;
  if type(f,'polynom'('rational',vars)) then RETURN(f)
    elif type(f,'polynom(rational)') then RETURN(expand(f)) fi;
  g:=collect(f,vars,'distributed',normal);
  if type(g,'ratpoly(rational)') then g
    elif type(g,`+`) then map(`SF/normal/zap`,g,vars,v)
    else `SF/normal/zap`(g,vars,v)
  fi
end:
#
# If the monomial part of f matches m, zap it.
# Otherwise pass it through unharmed.
#
`SF/normal/zap`:=proc(f,vars,m) local u;
  for u in vars do
    if degree(f,u)<>degree(m,u) then RETURN(f) fi
  od; 0
end:
#
# top(f)    convert the symmetric function f into a p-polynomial.
# top(f,b)  does the same, assuming f is expressed solely in basis b.
# The final result is collected with respect to p1,p2,p3,...
# If b is not a predefined basis, f must be *linear* in the b[..]'s.
#
`SF/top`:=proc() local f,bases,sp,i,mu,b,d,a;
  f:=args[1]; bases:=`SF/getbase`(args) minus {'p'};
  sp:=SF['varset'](f,bases minus {'h','e'});
  for b in bases minus {'h','e','s[]'} do
    a:=op(0,b); if sp[a]=[] then next fi;
    if assigned(`SF/dual`[a]) then
      f:=`SF/dualize`(a,sp[a],f)
    else
      f:=`SF/added`(a,sp[a],f); bases:=bases union {'e'}
    fi
  od;
  if member('s[]',bases) then
    f:=subs({seq(s[op(mu)]=`SF/jt_det`(mu,0),mu=sp['s'])},f);
    bases:={op(bases),'h','e'}
  fi;
  if member('h',bases) then f:=`SF/to_ehp`('h','p',f) fi;
  if member('e',bases) then f:=`SF/to_ehp`('e','p',f) fi;
  d:=SF['varset'](f,'p');
  collect(f,[seq(cat('p',i),i=1..d)],'distributed',normal);
end: 
#
# tos(f,<options>) converts f into a sum of Schur functions.
# The optional arguments can be given in any order:
# (1) an equation 'nrows=<integer>', where <integer> is a positive integer
#     that specifies that all calculations should take place in the ring
#     spanned by Schur functions with at most <integer> rows.
# (2) a basis name, indicating that f is expressed in terms of that basis.
# (3) list of partitions that support the Schur expansion of f.
# New in 2.4: we ignore (3). The new code determines the Schur-support of
# f dynamically so this option is pointless.
#
# Set infolevel[tos]:=2 in order to see information about the
# Schur function coefficients as they are developed.
#
`SF/tos`:=proc() local b,c,d,j,f,den,nrows,vars,mu,res,v;
  nrows:=`SF/getrows`(args);
  if nrows>0 then
    f:=SF['toe'](args); b:='e';
    d:=min(nrows,SF['stdeg'](f,b));
  else
    f:=SF['toh'](args); b:='h';
    d:=SF['stdeg'](f,b); nrows:=NULL;
  fi;
  vars:=[seq(cat(b,-j),j=-d..-1)];
  res:=0; f:=`SF/tos/numer`(f,vars,'den');
  while f<>0 do
    c:=tcoeff(f,vars,'v');
    mu:=`SF/shape`(v,b,d);
    f:=f-c*expand(linalg['det'](SF['jt_matrix'](mu,[],b,nrows)));
    f:=`SF/normal`(f,vars,v);
    if b='e' then mu:=SF['conjugate'](mu) fi;
    res:=normal(c/den)*s[op(mu)]+res;
    userinfo(2,tos,c,mu,nops(f));
  od; res
end:
#
# If f is a sym poly with ratpoly coefficients, extract a least common
# denominator for the coefficients, assign it to args[3], and
# return the result of rescaling (and normalizing) f by this factor.
# Similar to 'primpart', but *way* more space efficient.
# If the coeffs are not ratpolys, do nothing (use 1 as the denominator).
#
`SF/tos/numer`:=proc(f,vars) local cfs,tms,den;
  if type(f,'polynom'('rational',vars)) then
    den:=denom(f); assign(args[3],den); den*f
  elif type(f,'ratpoly(rational)') then
    cfs:=[coeffs(f,vars,'tms')];
    den:=lcm(op(map(denom,cfs))); assign(args[3],den);
    cfs:=map((x,y)->normal(y*x),cfs,den);
    convert(zip((x,y)->x*y,cfs,[tms]),`+`)
  else
    assign(args[3],1); f
  fi
end:
#
# varset(f,<bases>) will return a table whose entries describe the
# sets of variables from <bases> that occur in symmetric function f.
#
# <bases> may be a list or set of string names and indexed names. For
# an indexed name such as 's[]', there will be a table entry indexed
# by 's' consisting of a list of partitions that indicate the support
# of this basis in f. For a string name such as 'p', there will be a
# table entry indexed by 'p' equal to the largest n s.t. p.n occurs in f.
#
# We assume that string names of bases are single characters.
#
# If <bases> is a name (not a list or set), then the entry of the above
# table corresponding to this name, not the table itself, is returned.
# If the second argument is omitted, the default is <bases>=`SF/Bases`.
#
`SF/varset`:=proc(f) local inds,strs,one_base,b,v,res,i,digits;
  inds:={}; strs:={}; res:=table(); one_base:=false;
  if nargs=1 then b:=`SF/Bases` else b:=args[2] fi;
  if type(b,'name') then one_base:=true; b:={b} fi;
  for v in b do
    if type(v,'indexed') then
      inds:={op(inds),op(0,v)}; res[op(0,v)]:=NULL
      else strs:={op(strs),v}; res[v]:=0
    fi
  od;
  if nops(strs)>0 then digits:=table([seq(cat(i)=i,i=0..9)]) fi;
  for v in indets(f,'name') do
    if type(v,'indexed') then b:=op(0,v);
      if member(b,inds) then res[b]:=res[b],[op(v)] fi
    else b:=substring(v,1..1);
      if member(b,strs) then
        res[b]:=max(res[b],`SF/varset/deg`(v,digits)) fi
    fi
  od;
  for b in inds do res[b]:=[res[b]] od;
  if one_base then res[op(strs),op(inds)] else op(res) fi;
end:
#
# If deleting the first character of string v yields a positive
# integer n then return n; otherwise, return 0.
#
`SF/varset/deg`:=proc(v,digits) local i,j,n;
  n:=0;
  for i from 2 to length(v) do
    j:=substring(v,i..i);
    if not assigned(digits[j]) then RETURN(0) fi;
    n:=10*n+digits[j]
  od; n
end:  
#
# Check whether b or b[] is a known basis.
#
`SF/verify`:=proc(b)
  if member(b,`SF/Bases`) then b
    elif member(b[],`SF/Bases`) then b[]
    else ERROR(cat(b,` is not a known basis`))
  fi
end:
#
# `SF/apply`(b,f,...) verifies the existence of `to`.b, then 
#    applies it to f,...
#
`SF/apply`:=proc() local b,pr;
  b:=`SF/verify`(args[1]);
  if type(b,'indexed') then b:=op(0,b) fi;
  pr:=cat(`to`,b);
  if type(pr,'procedure') then pr(args[2..nargs])
    else SF[pr](args[2..nargs]) fi;
end:
# 
# zee(lambda) = the order of the S_n centralizer of a permutation of
#  cycle type lambda; i.e., 1^(m1)*m1!*2^(m2)*m2!*...
#
# zee(lambda,a)   = zee(lambda)*a^nops(lambda) 
# zee(lambda,q,t) = zee(lambda)*prod((1-q^(lambda_i))/(1-t^(lambda_i)))
#
`SF/zee`:=proc(mu) local res,m,i;
  m:=1; res:=convert(mu,`*`);
  if nargs=2 then
    res:=res*args[2]^nops(mu)
  elif nargs=3 then 
    res:=res*convert([seq((1-args[2]^i)/(1-args[3]^i),i=mu)],`*`)
  fi;
  for i from 2 to nops(mu) do
    if mu[i]<mu[i-1] then m:=1 else m:=m+1 fi;
    res:=res*m;
  od; res
end:
#
map(x->assign(SF[x],cat(`SF/`,x)),
  ['('Par')', '('add_basis')', '('char2sf')', '('conjugate')',
   '('dominate')', '('dual_basis')', '('evalsf')', '('hooks')',
   '('itensor')', '('jt_matrix')', '('nextPar')', '('omega')',
   '('plethysm')', '('scalar')', '('sf2char')', '('skew')', '('stdeg')',
   '('subPar')', '('theta')', '('toe')', '('toh')', '('top')',
   '('tos')', '('varset')', '('zee')']):
#
printf(`SF 2.4v loaded. Run 'withSF()' to use abbreviated names.\n`);