(*
 * SymmetricFunctions.m
 * 
 * This Mathematica package implements the five usual bases of the ring
 * of symmetric functions (polynomials): elementary, homogeneous, power
 * sum, monomial, and Schur functions.
 * 
 * This version: Sun Dec  5 12:38:26 CST 2004.               %DATE_TAG%
 * 
 * The current version of this package and some possibly useful
 * documentation is available from http://www.math.umn.edu/~drake/.
 * 
 * 
 *   (c) 2004 Dan Drake <drake@math.umn.edu>
 * 
 *   This program is free software; you can redistribute it and/or
 *   modify it under the terms of the GNU General Public License as
 *   published by the Free Software Foundation; either version 2 of the
 *   License, or (at your option) any later version.
 * 
 *   This program is distributed in the hope that it will be useful, but
 *   WITHOUT ANY WARRANTY; without even the implied warranty of
 *   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
 *   General Public License for more details.
 * 
 *   The GNU GPL is available at http://www.gnu.org/licenses/gpl.html.
 * 
 *)

(* 
 * TODO:
 *   skew functions - DONE for Schur! (are there others?)
 *   Littlewood-Richardson coefficients
 *   omega involution
 *   "variable-less" versions: just work with e[3], p[{2,1,1}], etc
 *   look at IntegerPartitions.m and see if that would make the h's
 *     faster
 *)

BeginPackage["SymmetricFunctions`", {"DiscreteMath`Combinatorica`",
  "Algebra`SymmetricPolynomials`"}]

Unprotect[p, e, h, m, s]

p::usage = "p[varlist, x] returns the power-sum symmetric function for
the variables in varlist. x may be an integer or a partition."

e::usage = "e[varlist, x] returns the elementary symmetric function for
the variables in varlist. x may be an integer or a partition; for
integers, e[varlist,x] returns SymmetricPolynomial[varlist,x]. See
Algebra`SymmetricPolynomials`."

h::usage = "h[varlist, x] returns the homogeneous symmetric function for
the variables in varlist. x may be an integer or a partition."

m::usage = "m[varlist, x] returns the monomial symmetric function for
the variables in varlist. x may be an integer or a partition; if x is an
integer, the function returned is for the partition with a single part
'x'."

s::usage = "s[varlist, x] returns the Schur function for
the variables in varlist. x may be an integer or a partition; if x is an
integer, the function returned is for the partition with a single part
'x'.

s[varlist, x, y] returns the skew Schur function for x/y. If x or y is
an integer, it is treated as above."

schurPS::usage = "schurPS[n,p] returns the principal specialization of
the Schur function for the partition p, using q, q^2, ..., q^n for the
variables."

SymmetricFunctionQ::usage = "symmetricQ[f, vars] returns True if f is a
symmetric polynomial in the variables in the list vars, and False
otherwise."

zlambda::usage = "zlambda[lambda] returns z_lambda for a partition
lambda. n!/z_lambda equals the number of permutations of type lambda in
the symmetric group on n letters."

sfip::usage = "sfip[a, b, vars] returns the standard symmetric function
inner product of a and b."

Begin["`Private`"]

(*
 * Power sum symmetric functions
 *)
p[l_List, k_Integer?Negative] := 0
p[l_List, k_Integer?NonNegative] := Map[(#1)^k &, Apply[Plus, l]]
p[l_List, q_?PartitionQ] := Apply[Times, Map[p[l, #1] &, q]]

(*
 * Elementary symmetric functions, same as SymmetricPolynomial
 *)
e[l_List, k_Integer?Negative] := 0
e[l_List, 0] := 1
e[l_List, k_Integer?NonNegative] := SymmetricPolynomial[l, k]
e[l_List, p_?PartitionQ] := Apply[Times, Map[e[l, #1] &, p]]

(*
 * Homogeneous symmetric functions
 *)
h[l_List, k_Integer?Negative] := 0
h[l_List, 0] := 1
h[l_List, k_Integer?NonNegative] := 
  Apply[Plus, Map[m[l, #1] &, Partitions[k]]]
h[l_List, p_?PartitionQ] := Apply[Times, Map[h[l, #1] &, p]]

(*
 * Monomial symmetric functions
 *)
m[l_List, k_Integer?Negative] := 0
m[l_List, 0] := 1
m[l_List, p_?PartitionQ] := 
  Module[{i, j, k, subsets, exponents, biglist},
    subsets = KSubsets[l, Length[p]];
    exponents = Permutations[p];
    biglist = Table[
        Product[(subsets[[i]][[k]])^(exponents[[j]][[k]]), {k, 1, Length[p]}],
        {j, 1, Length[exponents]}, {i, 1, Length[subsets]}];
    Nest[Apply[Plus, #1] &, biglist, 2]
    ]
m[l_List, k_Integer?Positive] := p[l, k]

(*
 * Schur functions via Jacobi - Trudi determinants
 *)
s[l_List, k_Integer?Negative] := 0
s[l_List, 0] := 1
s[l_List, p_?PartitionQ] := Module[{q = TransposePartition[p]},
  Det[Table[e[l, q[[i]] - i + j], {i, 1, Length[q]}, {j, 1, Length[q]}]]
  ]
s[l_List, k_Integer?Positive] := s[l, {k}]

(*
 * Skew Schur functions; just modify the determinant (from Macdonald)
 *)
s[l_List, p_, q_] := Which[
  IntegerQ[p] && IntegerQ[q], s[l,{p}, {q}],
  IntegerQ[p] && PartitionQ[q], s[l,{p}, q],
  PartitionQ[p] && IntegerQ[q], s[l, p, {q}],
  PartitionQ[p] && PartitionQ[q], Module[
    {pp = TransposePartition[p], qq = TransposePartition[q]},
    qq = PadRight[qq, Max[Length[pp], Length[qq]]];
    Det[Table[e[l, pp[[i]] - qq[[j]] - i + j],
      {i, 1, Length[pp]}, {j, 1, Length[pp]}]]
    ]
  ]

(*
 * So nasty hobbitses can't redefine our symmetric functions.
 *)
Protect[e, h, p, m, s]

(*
 * Principal specialization of Schur or skew Schur functions
 *)
SchurPS[n_Integer?NonNegative, p_] :=
  s[Table[q^k, {k,1,n}], p]
SchurPS[n_Integer?NonNegative, p_, q_] :=
  s[Table[q^k, {k,1,n}], p, q]

SymmetricFunctionQ[f_, vars_List] :=
  If[SymmetricReduction[f,vars][[2]] == 0, True, False]

zlambda[p_?PartitionQ] := Apply[Times, Table[If[Count[p, k] > 0,
  Count[p, k]! k^Count[p, k], 1], {k, 1, Max[p]}]]

(*
 * This computes the standard inner product on symmetric functions;
 * defined by 
 * 
 *   <p_a, p_b> = [a = b] z_a
 * 
 * where a and b are partitions. I should try and see if I can get a
 * version working where you don't have to pass it the variables, since
 * the answer doesn't depend on them.
 *)
sfip[a_, b_, vars_] := 
  Module[{u = convert2p[a, vars], v = convert2p[b, vars]}, 
    sfipexpand[Map[powersum2partition, 
        Switch[Head[u], p | Times | Power, {u}, Plus, Apply[List, u]]], 
      Map[powersum2partition, 
        Switch[Head[v], p | Times | Power, {v}, Plus, Apply[List, v]]]]
        ]

(*
 * Returns the symmetric function f as a polynomial in power sums p[k].
 * Since p is protected and p[k] undefined, the p[k]'s won't evaluate
 * and we can work with them that way.
 *)
convert2p[f_, vars_] := 
  SymmetricReduction[f, vars, Table[ep[i], {i, 1, Length[vars]}]][[1]] // 
    Expand

(* 
 * Returns the nth elementary symmetric function as a polynomial in p's.
 *)
ep[n_] := 
  Det[Table[Which[
    i-j+1 > 0, p[i-j+1],
    i-j+1 == 0, i,
    True, 0], {i, 1, n}, {j, 1, n}]]/n!

(*
 * This does the expand-by-linearity work of the inner product. You feed
 * it a and b, lists of two-element lists as described below, and it
 * expands by linearity exactly how you think it should.
 *)
sfipexpand[a_List, b_List] := Apply[Plus, Flatten[Outer[ip, a, b, 1]]]

(*
 * Given a monomial of p's, this returns a two-element list {k,p} where
 * k is the coefficient and p is the partition corresponding to the
 * monomial. For example, 3 p[3]^2 p[1] becomes {3, {3,3,1}}. 
 * 
 * We use the Switch[Head[x] stuff because, if given just p[k] or p[k]^n
 * and we 'Apply[List' it, we lose information (Apply[List,p[3]] is 3,
 * for instance). Everything works fine when given something with a Head
 * of Times, like p[1]p[2], since then Apply[List splits it up properly.
 * 
 * I've put the line for matching real coefficients last because the
 * Switch terminates as soon as it matches something, and it's not so
 * common to be getting real coefficients. This is very slightly more
 * efficient.
 *)
powersum2partition[x_] := Module[{a = {1, {}}, y, z},
  y = Switch[Head[x], p | Power, {x}, Times, Apply[List, x]];
  Map[Switch[Head[#1],
    Integer, a[[1]] *= #1,
    Rational, a[[1]] *= #1,
    p, AppendTo[a[[2]], Apply[Identity, #1]],
    Power, z = Apply[List, #1]; 
      a[[2]] = multiappend[Apply[Identity, z[[1]]], a[[2]], z[[2]]],
    Real, a[[1]] *= #1] &, y];
  Return[a]
  ]

(* 
 * The core of the inner product: we've represent p's by two-element
 * lists as mentioned above, and if they're equal, ip returns the
 * product of the coefficient times zlambda. Otherwise it returns zero.
 *)
ip[a_, b_] /; PartitionQ[a[[2]]] && PartitionQ[b[[2]]] := 
  If[Sort[a[[2]]] == Sort[b[[2]]], zlambda[a[[2]]] a[[1]]b[[1]], 0]

(*
 * A little helper function; returns the list l with n copies of x
 * appended to it.
 *)
multiappend[x_, l_List, n_Integer?NonNegative] :=
  Join[l, Table[x, {i,1,n}]]

End[]

EndPackage[]