Symmetric polynomial functions

sym.py

from sage.combinat.sf.sfa import SFAElementary
from sage.rings.all import PolynomialRing
from sage.misc.misc import prod

def is_symmetric(p):
    vars = p.variables()
    if len(vars) == 1: return True

    permutation = {}
    for i in range(len(vars)-1):
        permutation[str(vars[i])] = vars[i+1]
    permutation[str(vars[len(vars)-1])] = vars[0]

    if(p != p(**permutation)): return False

    permutation = {str(vars[0]): vars[1], str(vars[1]) : vars[0]}
    return p == p(**permutation)

def symmetrize(p):
    if not is_symmetric(p): raise Error(str(p) + " is not a symmetric polynomial.")
    vars = p.variables()
    nvars = len(vars)

    e = SFAElementary(p.base_ring())
    sigmas = [e([i]).expand(nvars, alphabet=vars) for i in range(1, nvars+1)]

    R = PolynomialRing(p.base_ring(), ['sigma_%s' % i for i in range(1, nvars+1)])
    sigma_vars = R.gens()

    def sym(f):
        if f == 0: return 0
        c = f.lc()
        degrees = f.lm().degrees()

        exps = [degrees[i]-degrees[i+1] for i in range(len(degrees)-1)]
        exps.extend(degrees[-1:])

        g = prod([sigma_vars[i]**exps[i] for i in range(len(exps))])
        gp = prod([sigmas[i]**exps[i] for i in range(len(exps))])

        return c*g + sym(f - c*gp)
    return sym(p)