Computing certain polynomials related to symmetric functions and $\lambda$-rings in Sage The definition of a $\lambda$-ring (https://en.wikipedia.org/wiki/%CE%9B-ring) makes use of certain "universal" polynomials $P_n$ and $P_{n,m}$, which basically give you the formulas for computing expressions such as (respectively) $\lambda^n(xy)$ or $\lambda^n(\lambda^m(x))$.
There is no need to know anything about $\lambda$-rings to compute those polynomials, they have elementary descriptions in terms of symmetric functions. Namely, let $(x_i)_{i\in \mathbb{N}}$ and $(y_i)_{i\in \mathbb{N}}$ be two families of indeterminates, and let $e_k$ (resp. $f_k$) be the $k$th elementary symmetric function in the $x_i$ (resp. the $y_i$). Then $P_{n,m}(T_1,\dots,T_{nm})\in \mathbb{Z}[T_1,\dots,T_{nm}]$ is the unique polynomial such that $P_{n,m}(e_1,\dots,e_{nm})$ is the coefficient of $t^n$ in
$$\prod_{1\leqslant i_1 < \dots < i_m\leqslant nm}(1+t\cdot x_{i_1}\cdots x_{i_m}),$$
and $P_n(S_1,\dots,S_n,T_1,\dots,T_n)\in \mathbb{Z}[S_1,\dots,S_n,T_1,\dots,T_n]$ is the unique polynomial such that $P_n(e_1,\dots,e_n,f_1,\dots,f_n)$ is the coefficient of $t^n$ in
$$\prod_{1\leqslant i,j\leqslant n}(1+t\cdot x_iy_j).$$
It turns out these polynomials are very annoying to actually compute by hand except for very small cases, and I haven't been able to find a reference with enough examples (usually only 3 or 4 examples are given in a book).
Is there a function in some computer algebra system (preferably Sage as I have access to it) to compute these polynomials? They do not see a huge amount of usage, but they are not completely unknown either, so I thought someone might have already implemented that.
If not, I would be very grateful if someone could help me write code to compute them (preferably in Sage, again). I'm not very proficient with those algebra systems, but I figure that it should be relatively easy for anyone who is used to them (especially with functions regarding symmetric polynomials).
 A: Here is a direct computation of $P_{n,m}$ from the definition, in SageMath:
def P(n, m):
    Rx = PolynomialRing(QQ, n*m, names='x')
    x = Rx.gens()
    Rt.<t> = PolynomialRing(Rx)
    from itertools import combinations
    f = prod(1 + t*prod(x[i[k]] for k in range(m)) for i in combinations(range(n*m), m))
    e = SymmetricFunctions(QQ).elementary()
    f_e = e.from_polynomial(f[n])
    RT = PolynomialRing(QQ, n*m, names=['T{}'.format(k+1) for k in range(n*m)])
    f_T = sum(c*prod(RT.gen(k-1) for k in m) for (m,c) in f_e)
    return f_T

For example:
sage: P(1,1)
T1
sage: P(1,2)
T2
sage: P(2,1)
T2
sage: P(2,2)
T1*T3 - T4
sage: P(2,3)
T2*T4 - T1*T5 + T6
sage: P(3,2)
T1^2*T4 + T3^2 - 2*T2*T4 - T1*T5 + T6


And here is $P_n$ directly from the definition, in SageMath:
def PP(n):
    Rx = PolynomialRing(QQ, n, names='x')
    Ry = PolynomialRing(Rx, n, names='y')
    from itertools import product
    Rt.<t> = PolynomialRing(Ry)
    f = prod(1 + t*Rx.gen(i)*Ry.gen(j) for i,j in product(range(n),repeat=2))
    SymRx = SymmetricFunctions(Rx).elementary()
    SymRx._prefix = 'f'
    f_f = SymRx.from_polynomial(f[n])
    Sym = SymmetricFunctions(QQ).elementary()
    RST = PolynomialRing(QQ, 2*n, names=['S{}'.format(k+1) for k in range(n)] + ['T{}'.format(k+1) for k in range(n)])
    S = RST.gens()[:n]
    T = RST.gens()[n:]
    result = 0
    for m, c in f_f:
        prefactor = prod(T[k-1] for k in m)
        for m2, c2 in Sym.from_polynomial(c):
            result += prefactor*c2*prod(S[k-1] for k in m2)
    return result

For example:
sage: PP(1)
S1*T1
sage: PP(2)
S2*T1^2 + S1^2*T2 - 2*S2*T2
sage: PP(3)
S3*T1^3 + S1*S2*T1*T2 + S1^3*T3 - 3*S3*T1*T2 - 3*S1*S2*T3 + 3*S3*T3

