#
# Symmetric Functions package: SF version 2.3, vanilla edition.
# This version/edition requires Maple V R1, R2, R3, R4, R5, or Maple 6.
#
# 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) 2001 by John R. Stembridge
#
#########################################################################
#
SF:=table():
e:='e': p:='p': h:='h': s:='s':
`SF/Bases`:={e,h,p,s[]}:
#
# assign short names, printing warnings if conflicts occur.
#
withSF:=proc() local install,f,lpr;
  if [op(I)]<>[1] then lpr:=lprint else # maple6 has a broken lprint
    lpr:=proc(x,y) printf(`%0.70s\n`,cat(x,`  `,y)) end fi;
  install:=proc(x,wrn)
    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 wrn(`Warning: new definition for`,x) fi;
      assign(x,SF[x])
    fi; x
  end;
  if nargs>0 then map(install,[args],lpr) else
    f:=proc() map(op,[args]) end; # hack the names w/o full evaluation!
    map(install,f(indices(SF)),lpr)
  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.
#
`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 RETURN([[]]) fi;
  if col=0 then RETURN([]) fi;
  [seq(op(map(proc(x,y) [y,op(x)] end,`SF/Par/sub`(n+i,-i,col-1),-i)),
    i=-min(row,n)..-iquo(n+col-1,col))];
end:
#
# add_basis(<name>,<inner_product>,<normalization>)
#  <name> = variable name for indexing terms of the new basis
#  <inner_product> = a procedure that accepts partitions as input, and
#    returns the scalar product of the corresponding power-sum with itself 
#  (optional) <normalization> = a procedure that accepts partitions as input,
#    returns the desired %leading% term of the basis element indexed by the
#    given partition.
#  The default normalization is to take the leading coefficients to be 1.
#
#  `SF/Bases`, `SF/iprod`, and `SF/lcoeff` are subversively global
#
`SF/add_basis`:=proc(b,ip) local a;
  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;
  assign(evaln(`SF/Bases`)={op(`SF/Bases`), a[]});
  assign(`SF/iprod`[a]=ip);
  if nargs=2 then assign(`SF/lcoeff`[a]=1) else
    assign(`SF/lcoeff`[a]=args[3]) fi;
  assign(cat(`to`,b)=proc() local c,d,f,B,i,sp,vars;
    if not type(args[nargs],'list') then
      f:=map(procname,[`SF/homog_cmps`(args[1])],args[2..nargs],[]);
      RETURN(convert(f,`+`))
    fi;
    B:=substring(procname,3..length(procname)); 
    f:=SF['toe'](args[1..nargs-1]); d:=SF['stdeg'](f,'e');
    vars:={seq(cat('e',i),i=1..d)}; sp:=args[nargs];
    if sp=[] then sp:=SF['Par'](d,degree(f,vars),d) fi;
    f:=convert([seq(c[i]*`SF/toe/extra`(B,sp[i]),i=1..nops(sp))],`+`)-f; 
    f:=collect(f,vars,'distributed');
    f:=readlib(`solve/linear`)({coeffs(f,vars)},{seq(c[i],i=1..nops(sp))});
    subs(f,convert([seq(c[i]*B[op(sp[i])],i=1..nops(sp))],`+`));
  end);
  'Okay';
end: 
#
# char2sf(<classfn>) will apply the characteristic map to <classfn>, producing
# a symmetric function as the output. <classfn> must be expressed as a linear
# combination of characteristic functions cl[<mu>] for various partitions <mu>.
# The output, by default, will be a p-polynomial.
# Use char2sf(<classfn>,b) to specify another base.
#
`SF/char2sf`:=proc(char) local res,mu;
  res:=SF['varset'](char,'cl[]');
  res:=convert([seq(coeff(char,cl[op(mu)])*
    convert(map(x->cat('p',x),mu),`*`)/SF['zee'](mu),mu=res)],`+`);
  if nargs>1 then `SF/verify/apply`(args[2],res,'p') else res fi;
end: 
#
# conjugate(lambda) returns the conjugate of partition lambda.
# lambda must be a list in decreasing order.
#
`SF/conjugate`:=proc(mu) local i,l;
  l:=nops(mu); if l=0 then RETURN([]) fi;                                         [l$mu[l],seq((-i)$(mu[-i]-mu[1-i]),i=1-l..-1)]
end:
#
# dominate(lambda) return a list of partitions <= lambda in dominance order
# dominate(lambda,n) does the same but only for partitions with <= n rows
#
`SF/dominate`:=proc(mu) local n,nu,res,sat,i,j,m,nu0,lam;
  n:=convert(mu,`+`);
  if nargs>1 then n:=min(n,args[2]) fi;
  if n<nops(mu) then RETURN([]) fi;
  res:=[mu]; sat:=0;
  while sat<nops(res) do
    nu:=res[sat+1]; 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;
    sat:=sat+1;
  od;
  res
end:
#
# dual_basis(<name1>,<name2>,<inner_prod>)
#  <name1> = variable name for indexing terms of the new basis
#  <name2> = name of the basis whose dual is desired 
#  <inner_prod> = a procedure that accepts partitions as input, and
#    returns the scalar product of the corresponding power-sum with itself 
#
# `SF/Bases`, `SF/dual`, and `SF/iprod` are subversively global
#
`SF/dual_basis`:=proc(b1,b2) local a;
  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;
  assign(`SF/dual`[a]=`SF/verify`(b2));
  assign(evaln(`SF/Bases`)={op(`SF/Bases`), a[]});
  if nargs>2 then assign(`SF/iprod`[a]=args[3])
    else assign(`SF/iprod`[a]=SF['zee']) fi;
  assign(cat(`to`,a)=proc() local c,d,f,B,i,j,sp,b,bb,res,n;
    if not type(args[nargs],'list') then
      f:=map(procname,[`SF/homog_cmps`(args[1])],args[2..nargs],[]);
      RETURN(convert(f,`+`))
    fi;
    B:=substring(procname,3..length(procname)); 
    f:=SF['top'](args[1..nargs-1]); d:=SF['stdeg'](f,'p');
    sp:=args[nargs]; if sp=[] then sp:=SF['Par'](d) fi;
    n:=nops(sp); b:=`SF/dual`[B];
    if type(b,'indexed') then bb:=op(0,b);
      res:=convert([seq(c[i]*bb[op(sp[i])],i=1..n)],`+`)
    else
      res:=convert([seq(c[i]*convert([seq(cat(b,j),j=sp[i])],`*`),
        i=1..n)],`+`)
    fi;
    f:=SF['scalar'](f,res,'p',b,`SF/iprod`[B]);
    f:=collect(f,[seq(c[i],i=1..n)],'distributed',normal);
    subs({seq(c[i]=B[op(sp[i])],i=1..n)},f);
  end);
  'Okay';
end: 
#
# evalsf(f,a) generates a plethystic evaluation of f at a as follows:
# 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 appearing in the 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,...
#
`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);
  if da=0 then expr else
    collect(expr,[seq(cat('p',i),i=1..da*df)],'distributed',normal)
  fi;
end:
#
# 
# homog_cmps(f)   separates the symmetric function f into its homogeneous
#  components. The components are not sorted by degree.
# homog_cmps(f,b) same, but assume 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;
  coeffs(collect(f,t),t);
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() local mu,nu,i,j,z,w;
  mu:=args[1]; nu:=SF['conjugate'](mu);
  if nargs=1 then
    mu:=[seq(seq(mu[i]-i+nu[j]-j+1,j=1..mu[i]),i=1..nops(mu))];
    sort(mu,proc(x,y) evalb(x>y) end)
  elif nargs=2 then
    z:=args[2];
    mu:=[seq(seq(z*(mu[i]-j)+nu[j]-i+1,j=1..mu[i]),i=1..nops(mu))];
    convert(mu,`*`)
  else
    mu:=[seq(seq((1-z^(mu[i]-j)*w^(nu[j]-i+1),j=1..mu[i])),i=1..nops(mu))];
    subs(z=args[2],w=args[3],convert(mu,`*`));
  fi;
end:
#
# itensor(f,g)         compute the inner tensor product (a.k.a. the internal
# itensor(f,g,b1,b2)   product) of symmetric functions f and g. Use b1 and b2
# itensor(f,g,b)       to specify bases used for f and g. Use b to specify
# itensor(f,g,b1,b2,b) the desired output basis (Default = 'p'.)
#
`SF/itensor`:=proc() local b,d,i,j,k,f,cfs,tms,res,vars,mu;
  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)];
  if nops(f[1])>nops(f[2]) then f:=[f[2],f[1]] fi; # a dirty hack
  d:=SF['varset'](f[1]*f[2],'p');
  vars:=[seq(cat('p',i),i=1..d)];
  cfs:=[seq([coeffs(f[i],vars,tms[i])],i=1..2)];
  tms:=[[tms[1]],[tms[2]]]; res:=0;
  for i to nops(tms[1]) do
    if member(tms[1][i],tms[2],'j') then
      mu:=[seq((-k)$(degree(tms[1][i],vars[-k])),k=-d..-1)];
      res:=res+SF['zee'](mu)*normal(cfs[1][i]*cfs[2][j])*tms[1][i];
    fi;
  od;
  if modp(nargs,2)=1 then `SF/verify/apply`(args[nargs],res,'p') else res fi;
end:
#
# jt_matrix(lambda)     will produce the Jacobi-Trudi matrix for the
#   partition lambda.
# jt_matrix(lambda,mu)  will do the same for the skew shape defined by the
#   pair of partitions lambda,mu (mu = inner shape).
#
`SF/jt_matrix`:=proc() local b,n,i,j,mu,nu,entry;
  mu:=args[1];
  if nargs>1 then nu:=args[2] else nu:=[] fi;
  if nargs>2 then b:=args[3] else b:='h' fi;
  n:=max(1,nops(mu),nops(nu));
  mu:=[op(mu),0$n]; nu:=[op(nu),0$n];
  entry:=proc(x,y) if x>0 then cat(y,x) elif x=0 then 1 else 0 fi end;
  array([seq([seq(entry(mu[i]-i+j-nu[j],b),j=1..n)],i=1..n)]);
end:
#
# omega(f)       Apply the omega-automorphism to symmetric function f. If a
# omega(f,b1)    second argument is specified, omega assumes that f is
# omega(f,b1,b2) expressed in terms of the basis b1. If a third argument is
#   present, the output is converted to base b2; otherwise, omega() will use
#   any basis or mix of bases that it finds convenient.
#
`SF/omega`:=proc() local f,sp,j,b,a,b0,mu,b1,bases;
  f:=args[1]; b0:={'p','h','e','s[]'};
  if nargs=1 then bases:=`SF/Bases`
    else b1:=`SF/verify`(args[2]); bases:={b1}
  fi;
  sp:=SF['varset'](f,bases minus b0);
  for b in bases minus b0 do a:=op(0,b);
    if assigned(`SF/dual`[a]) then b1:='p' else b1:='e' fi;
    if nargs=1 then
      f:=subs({seq(a[op(mu)]=SF[cat(`to`,b1)](a[op(mu)],b),mu=sp[a])},f)
      else f:=SF[cat(`to`,b1)](f,b)
    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'])},f);
  f:=subs({seq(cat('p',j)=(-1)^(j-1)*cat('p',j),j=1..sp['p'])},f);
  f:=subs({seq(s[op(mu)]=s[op(SF['conjugate'](mu))],mu=sp['s'])},f);
  if nargs>2 then `SF/verify/apply`(args[3],f,b1) else f fi;
end:
#
# plethysm(f,g)         compute the plethysm f[g] of symmetric functions f,g.
# plethysm(f,g,b1,b2)   Use b1 and b2 to specify bases used for f and g.
# plethysm(f,g,b)       Use b to specify the desired output basis.
# plethysm(f,g,b1,b2,b) The default for b is 'p'.
#
`SF/plethysm`:=proc() local d,i,f,g,b,res;
  if nargs>3 then b:=[args[3..4]] else b[1]:=NULL; b[2]:=NULL fi;
  f:=SF['top'](args[1],b[1]); 
  d:=SF['varset'](f,'p');
  if f=cat('p',d) then
    g:=SF['top'](args[2],b[2]);
    res:=subs({seq(cat('p',i)=cat('p',d*i),i=1..SF['varset'](g,'p'))},g);
  else
    if b[1]<>NULL then b:=subsop(1='p',b) fi;
    res:=subs({seq(cat('p',i)=
      `SF/plethysm`(cat('p',i),args[2],b[1],b[2]),i=1..d)},f);
  fi;
  d:=SF['varset'](res,'p');
  res:=collect(res,[seq(cat('p',i),i=1..d)],'distributed',normal);
  if modp(nargs,2)=1 then `SF/verify/apply`(args[nargs],res,'p') else res fi;
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).
#
# Options: scalar(f,g,'b1','b2') will compute scalar(f,g) under the
# assumption that f is in base b1 and g is in base b2. 
# Remember that if b1 is an indexed basis (other than s[]), then f must be
# a linear combination of the b1[mu]'s.)
#
# If the last argument (the third or fifth) is a procedure f() that accepts
# partitions as input, then the scalar product will be computed with respect
# to the form for which power sums are orthogonal and <p[mu],p[mu]>=f(mu).
#
`SF/scalar`:=proc() local f,ip,b,d,vars,tms,cfs,mu,i,j,k,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['top'](args[i],b[i]),i=1..2)];
  if nops(f[1])>nops(f[2]) then f:=[f[2],f[1]] fi; # a dirty hack
  d:=SF['varset'](f[1]*f[2],'p');
  vars:=[seq(cat('p',i),i=1..d)];
  cfs:=[seq([coeffs(f[i],vars,tms[i])],i=1..2)];
  tms:=[[tms[1]],[tms[2]]]; res:=0;
  for i to nops(tms[1]) do
    if member(tms[1][i],tms[2],'j') then
      mu:=[seq((-k)$(degree(tms[1][i],vars[-k])),k=-d..-1)];
      res:=res+ip(mu)*cfs[1][i]*cfs[2][j];
    fi;
  od;
  res
end:
#
# sf2char(f) will apply the (inverse) characteristic map to the symmetric
#   function f. The result is expressed as a linear combination of
#   characteristic functions cl[<mu>] for various partitions <mu>.
# sf2char(f,b) does the same, assuming that f is in base b.
#
`SF/sf2char`:=proc() local poly,i,j,d,res,cfs,term,mu;
  poly:=SF['top'](args);
  d:=SF['varset'](poly,'p');
  cfs:=[coeffs(poly,[seq(cat('p',i),i=1..d)],'term')];
  term:=[term]; res:=0;
  for i from 1 to nops(cfs) do;
    mu:=seq((d-j)$(degree(term[i],cat('p',d-j))),j=0..d-1);
    res:=res+SF['zee']([mu])*cfs[i]*cl[mu];
  od;
  res;
end: 
#
# skew(f,g)               Apply the linear transformation f^* to g, where f^*
# skew(f,g,<iprod>)       denotes the adjoint to multiplication by f. If b1
# skew(f,g,b1,b2)         and b2 are specified, assume that f and g are in
# skew(f,g,b1,b2,<iprod>) bases b1 and b2, resp. The output is expressed in
#  base p. A third or fifth argument can be a procedure specifying the scalar
#  product for the adjoint. The default scalar product is zee().
# 
`SF/skew`:=proc() local f,ip,b,d,i,res,g1,g2;
  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)];
  res:=[seq(seq(`SF/skew/hg`(g1,g2,ip),g1=[`SF/homog_cmps`(f[1],'p')]),
    g2=[`SF/homog_cmps`(f[2],'p')])];
  res:=convert(res,`+`); d:=SF['varset'](res,'p');
  collect(res,[seq(cat('p',i),i=1..d)],'distributed',normal);
end:
#
`SF/skew/hg`:=proc(f,g,ip) local d,c,i,sp,n,res,pm;
  d:=SF['stdeg'](g,'p')-SF['stdeg'](f,'p');
  if d<0 then RETURN(0) fi;
  sp:=SF['Par'](d); n:=nops(sp);
  pm:=[seq(convert(map(x->cat('p',x),i),`*`),i=sp)];
  res:=convert([seq(c[i]*pm[i]/ip(sp[i]),i=1..n)],`+`);
  res:=SF['scalar'](g,res*f,'p','p',ip);
  res:=collect(res,[seq(c[i],i=1..n)],'distributed',normal);
  subs({seq(c[i]=pm[i],i=1..n)},res);
end:
#
# stdeg(f)   determine the degree of f with respect to 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) local n,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,`+`); n:=args[2];
    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 symmetric function f.
# theta(f,q,t) Apply the automorphism p.j -> (1-q^j)/(1-t^j)*p.j to the
#   symmetric function f.
#
`SF/theta`:=proc(g,q,t) local f,d,j;
  f:=SF['top'](g); d:=SF['varset'](f,'p');
  if nargs=2 then
    f:=subs({seq(cat('p',j)=q*cat('p',j),j=1..d)},f)
  elif nargs>2 then
    f:=subs({seq(cat('p',j)=(1-q^j)/(1-t^j)*cat('p',j),j=1..d)},f)
  else
    ERROR(`wrong number of parameters`)
  fi;
  map(normal,f);
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,...
# For d<15, it is optimized for speed; for d>=15, it is slower but tries
#  to minimize space.
#
`SF/to_ehp`:=proc(b1,b2) local d,d2,z,i,f,hack,vars,inds,g;
  f:=args[3]; d:=SF['varset'](f,b1);
  if d=0 then RETURN(f) fi;
  if nargs>3 then d2:=min(d,args[4]) else d2:=d fi;
  if b1='e' then if b2='p' then
      g:=exp(convert([seq(-cat('p',i)*(-z)^i/i,i=1..d2)],`+`))
    else
      g:=1/(1+convert([seq(cat('h',i)*(-z)^i,i=1..d2)],`+`))
    fi
  elif b1='p' then if b2='e' then
      g:=-ln(1+convert([seq(cat('e',i)*(-z)^i,i=1..d2)],`+`))
    else
      g:=ln(1+convert([seq(cat('h',i)*z^i,i=1..d2)],`+`))
    fi;
    g:=z*diff(taylor(g,z,d+1),z)
  elif b2='e' then
      g:=1/(1+convert([seq(cat('e',i)*(-z)^i,i=1..d2)],`+`))
    else
      g:=exp(convert([seq(cat('p',i)*z^i/i,i=1..d2)],`+`))
  fi;
  hack:=taylor(g,z,d+1);
  if d<15 then
    f:=subs({seq(cat(b1,i)=coeff(hack,z,i),i=1..d)},f)
  else
    vars:=seq(cat(b1,i),i=1..d);
    inds:=indets(f) intersect {vars};
    vars:=[seq(cat(b2,i),i=1..d2),vars];
    for i to d do
      if not member(cat(b1,i),inds) then next fi;
      f:=subs(cat(b1,i)=expand(coeff(hack,z,i)),f);
      f:=collect(f,vars,'distributed');
    od
  fi; f
end:
#
`SF/to_ehp/det`:=proc(nu,n) local i,g,b;
  if nu=[] then RETURN(1) fi;
  if nargs>2 then b:=args[3] elif nops(nu)>nu[1] then b:='e' fi;
  if b='e' then g:=SF['conjugate'](nu),[],'e' else g:=nu fi;
  g:=SF['jt_matrix'](g);
  if n>0 then g:=subs({seq(cat('e',i)=0,i=n+1..nops(nu)+nu[1]-1)},op(g)) fi;
  linalg['det'](g);
end:
#
# toe(f) converts the symmetric function f into an e-polynomial.
# toe(f,b), does the same, assuming that f is expressed solely in terms of
#   the base b. The final result is collected with respect to e1,e2,e3,...
# If b is not a predefined basis, f must be a *linear* function of the b[..]'s. 
#
`SF/toe`:=proc() local poly,bases,sp,i,mu,d,b,c,nrows,pref;
  poly:=args[1]; bases:=`SF/Bases`; nrows:=0;
  if type(args[nargs],`=`) then
    nrows:=op(2,args[nargs]);
    if nargs>2 then bases:={`SF/verify`(args[2])} fi;
  elif nargs>1 then bases:={`SF/verify`(args[2])}
  fi;
  bases:=bases minus {'e'};
  sp:=SF['varset'](poly,bases minus {'p','h'});
  for b in bases minus {'h','p','s[]'} do;
    c:=op(0,b); if sp[c]=[] then next fi;
    if not assigned(`SF/dual`[c]) then
      poly:=subs({seq(c[op(mu)]=`SF/toe/extra`(c,mu),mu=sp[c])},poly)
    elif nops(bases)=1 then 
      poly:=SF['top'](poly,b); bases:=bases union {'p'};
    else
      poly:=subs({seq(c[op(mu)]=SF['top'](c[op(mu)],b),mu=sp[c])},poly);
      bases:=bases union {'p'}
    fi;
  od;
  if member('s[]',bases) then pref:='e';
    if nops(sp['s'])>1 then pref:=NULL; bases:=bases union {'h'} fi;
    poly:=subs({seq(s[op(mu)]=
      `SF/to_ehp/det`(mu,nrows,pref),mu=sp['s'])},poly);
  fi;
  if nrows>0 then d:=SF['varset'](poly,'e');
    poly:=subs({seq(cat('e',i)=0,i=nrows+1..d)},poly)
    else nrows:=NULL
  fi;
  if member('p',bases) then poly:=`SF/to_ehp`('p','e',poly,nrows) fi;
  if member('h',bases) then poly:=`SF/to_ehp`('h','e',poly,nrows) fi;
  d:=SF['varset'](poly,'e');
  collect(poly,[seq(cat('e',i),i=1..d)],'distributed',normal);
end: 
#
`SF/toe/extra`:=proc(b,mu) local n,sp,res,pres,vars,i,c,j,eq;
  option remember;
  c[1]:=`SF/lcoeff`[b](mu); n:=convert(mu,`+`);
  if mu=[] then RETURN(c[1]) fi;
  sp:=SF['Par'](n,mu[1],n); member(mu,sp,'i');
  sp:=[seq(convert(map(x->cat('e',x),SF['conjugate'](sp[j])),`*`),
    j=i..nops(sp))];
  res:=convert([seq(c[i]*sp[i],i=1..nops(sp))],`+`);
  pres:=SF['top'](res,'e'); vars:={seq(c[i],i=2..nops(sp))};
  eq:={seq(collect(SF['scalar'](sp[i],pres,'e','p',`SF/iprod`[b]),
    vars,'distributed',normal),i=2..nops(sp))};
  subs(normal(readlib(`solve/linear`)(eq,vars)),res);
end:
#
# toh(f) converts the symmetric function f into an h-polynomial.
# toh(f,b), does the same, assuming that f is expressed solely in terms of
#   the base b. The final result is collected with respect to h1,h2,h3,...
# If b is not a predefined basis, f must be a *linear* function of the b[..]'s.
#
`SF/toh`:=proc() local poly,bases,sp,i,mu,d,b,c,pref;
  if nargs=1 then bases:=`SF/Bases` else bases:={`SF/verify`(args[2])} fi;
  poly:=args[1]; bases:=bases minus {'h'};
  sp:=SF['varset'](poly,bases minus {'p','e'});
  for b in bases minus {'e','p','s[]'} do;
    c:=op(0,b); if sp[c]=[] then next fi;
    if not assigned(`SF/dual`[c]) then
      poly:=subs({seq(c[op(mu)]=SF['toe'](c[op(mu)],b),mu=sp[c])},poly);
      bases:=bases union {'e'}
    elif nops(bases)=1 then 
      poly:=SF['top'](poly,b); bases:=bases union {'p'};
    else
      poly:=subs({seq(c[op(mu)]=SF['top'](c[op(mu)],b),mu=sp[c])},poly);
      bases:=bases union {'p'}
    fi;
  od;
  if member('s[]',bases) then pref:='h';
    if nops(sp['s'])>1 then pref:=NULL; bases:=bases union {'e'} fi;
    poly:=subs({seq(s[op(mu)]=`SF/to_ehp/det`(mu,0,pref),mu=sp['s'])},poly);
  fi;
  if member('p',bases) then poly:=`SF/to_ehp`('p','h',poly) fi;
  if member('e',bases) then poly:=`SF/to_ehp`('e','h',poly) fi;
  d:=SF['varset'](poly,'h');
  collect(poly,[seq(cat('h',i),i=1..d)],'distributed',normal);
end: 
#
# top(f) converts the symmetric function f into an p-polynomial.
# top(f,b), does the same, assuming that f is expressed solely in terms of
#   the base b. The final result is collected with respect to p1,p2,p3,...
# If b is not a predefined basis, f must be a *linear* function of the b[..]'s.
#
`SF/top`:=proc() local poly,bases,sp,i,mu,b,d,a;
  if nargs=1 then bases:=`SF/Bases` else bases:={`SF/verify`(args[2])} fi;
  poly:=args[1]; bases:=bases minus {'p'};
  sp:=SF['varset'](poly,bases minus {'h','e'});
  for b in bases minus {'h','e','s[]'} do;
    a:=op(0,b); if sp[a]=[] then next fi;
    if not assigned(`SF/dual`[a]) then
      poly:=subs({seq(a[op(mu)]=SF['toe'](a[op(mu)],b),mu=sp[a])},poly);
      bases:=bases union {'e'}
    elif nops(bases)=1 then
      poly:=convert(map(`SF/top/dual`,[`SF/homog_cmps`(poly,b)],a),`+`)
    else
      poly:=subs({seq(a[op(mu)]=`SF/top/dual`(a[op(mu)],a),mu=sp[a])},poly)
    fi
  od;
  if member('s[]',bases) then
    poly:=subs({seq(s[op(mu)]=`SF/to_ehp/det`(mu,0),mu=sp['s'])},poly);
    bases:=bases union {'h','e'}
  fi;
  if member('h',bases) then poly:=`SF/to_ehp`('h','p',poly) fi;
  if member('e',bases) then poly:=`SF/to_ehp`('e','p',poly) fi;
  d:=SF['varset'](poly,'p');
  collect(poly,[seq(cat('p',i),i=1..d)],'distributed',normal);
end: 
#
`SF/top/dual`:=proc(f,a) local c,d,i,sp,pm,res,b,du,mu,cfs,tms;
  b:=`SF/dual`[a]; d:=SF['stdeg'](f,a[]); sp:=SF['Par'](d);
  pm:=[seq(convert(map(x->cat('p',x),mu),`*`),mu=sp)];
  du:=[seq(c[i]*pm[i]/`SF/iprod`[a](sp[i]),i=1..nops(sp))];
  du:=`SF/verify/apply`(b,convert(du,`+`),'p');
  if type(b,'indexed') then b:=op(0,b);
    res:=convert([seq(coeff(f,a[op(mu)])*coeff(du,b[op(mu)]),
      mu=SF['varset'](f,a[]))],`+`)
  else
    cfs:=[coeffs(du,[seq(cat(b,i),i=1..d)],tms)];
    tms:=[tms]; res:=0;
    for mu in SF['varset'](f,a[]) do
      member(convert(map((x,y)->cat(y,x),mu,b),`*`),tms,'i');
      res:=res+coeff(f,a[op(mu)])*cfs[i];
    od
  fi;
  res:=collect(res,[seq(c[i],i=1..nops(sp))],'distributed',normal);
  subs({seq(c[i]=pm[i],i=1..nops(sp))},res);
end:
#
# tos(f,<options>) will convert the symmetric function f into a sum of
# Schur functions. <options> is a sequence of zero or more of the
# following expressions (in any order):
# (1) a list of partitions that support the Schur expansion of f.
# (2) 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.
# (3) a name 'b' that indicates what basis f is expressed in. 
#
`SF/tos`:=proc() local c,j,d,f,sp,n,nrows,inpt,opt,inds;
  inpt:=args[1]; nrows:=0;
  for opt in [args[2..nargs]] do;
    if type(opt,'list') then sp:=opt
      elif type(opt,`=`) then nrows:=op(2,opt)
      elif type(opt,'name') then inpt:=inpt,opt
    fi;
  od;
  if not assigned(sp) then
    f:=map(`SF/tos`,[`SF/homog_cmps`(args[1])],args[2..nargs],[]);
    RETURN(convert(f,`+`))
  elif nrows>0 then
    f:=SF['toe'](inpt,'nrows'=nrows);
    d:=SF['stdeg'](f,'e'); nrows:=min(nrows,d);
    inds:={seq(cat('e',j),j=1..nrows)};
    if sp=[] then sp:=SF['Par'](d,degree(f,inds),nrows)
      else sp:=map(proc(x,y) if nops(x)<y then x fi end,sp,nrows+1) fi;
    n:=nops(sp);
    inpt:=[seq(c[j]*`SF/to_ehp/det`(sp[j],nrows,'e'),j=1..n)];
  else
    f:=SF['toh'](inpt); d:=SF['stdeg'](f,'h');
    inds:={seq(cat('h',j),j=1..d)};
    if sp=[] then sp:=SF['Par'](d,degree(f,inds)) fi;
    n:=nops(sp);
    inpt:=[seq(c[j]*linalg['det'](SF['jt_matrix'](sp[j])),j=1..n)];
  fi;
  inpt:=convert(inpt,`+`)-f;
  inpt:={coeffs(collect(inpt,inds,'distributed'),inds)};
#  inpt:=readlib(`solve/linear`)(inpt,{seq(c[j],j=1..n)});
  inpt:=SolveTools:-Linear(inpt,{seq(c[j],j=1..n)});
  f:=subs(inpt,[seq(c[j]*s[op(sp[j])],j=1..n)]);
  convert(map(normal,f),`+`);
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 each
# indexed name (e.g., '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 each string name (e.g.,'p'), there will be a table
# entry indexed by 'p' equal to the largest n s.t. p.n occurs in f.
#
# If <bases> is a single 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,x,res;
  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 x in b do
    if type(x,'indexed') then
      inds:={op(inds),op(0,x)}; res[op(0,x)]:=NULL
      else strs:={op(strs),x}; res[x]:=0
    fi
  od;
  for x in indets(f,'name') do
    if type(x,'indexed') then b:=op(0,x);
      if member(b,inds) then res[b]:=res[b],[op(x)] fi
    else b:=substring(x,1..1);
      if member(b,strs) then res[b]:=max(res[b],`SF/varset/deg`(x,b)) fi
    fi
  od;
  for x in inds do res[x]:=[res[x]] od;
  if one_base then res[op(strs),op(inds)] else op(res) fi;
end:
#
`SF/varset/deg`:=proc(v,b) local i,n,pos,digits;
  n:=0; digits:=[`0`,`1`,`2`,`3`,`4`,`5`,`6`,`7`,`8`,`9`];
  for i from length(b)+1 to length(v) while 
    member(substring(v,i..i),digits,'pos') do n:=10*n+pos-1 od;
  if i>length(v) then n else 0 fi;
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/verify/apply`(b,f,...) verifies the existence of `to`.b, then 
#    applies it to f,...
#
`SF/verify/apply`:=proc() local b;
  b:=`SF/verify`(args[1]);
  if type(b,'indexed') then b:=op(0,b) fi;
  if type(cat(`to`,b),'procedure') then cat(`to`,b)(args[2..nargs])
    else SF[cat(`to`,b)](args[2..nargs]) fi;
end:
#
# zee(lambda)=the order of the centralizer in S_n of a permutation of
# cycle type lambda  =  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')', '('omega')', '('plethysm')', '('scalar')', '('sf2char')',
  '('skew')', '('stdeg')', '('subPar')', '('theta')', '('toe')', '('toh')',
  '('top')', '('tos')', '('varset')', '('zee')']):
if [op(I)]=[1] then # we are in maple6 and lprint is broken
  printf(`%0.70s\n`,
         `SF 2.3v loaded. Run 'withSF()' to use abbreviated names`)
else
  lprint(`SF 2.3v loaded. Run 'withSF()' to use abbreviated names`)
fi;