pet_cycleind_symm :=
proc(n)
option remember;

    if n=0 then return 1; fi;

    expand(1/n*add(a[l]*pet_cycleind_symm(n-l), l=1..n));
end;


pet_flatten_term :=
proc(varp)
local terml, d, cf, v;

    terml := [];

    cf := varp;
    for v in indets(varp) do
        d := degree(varp, v);
        terml := [op(terml), seq(v, k=1..d)];
        cf := cf/v^d;
    od;

    [cf, terml];
end;

T :=
proc(n, k)
option remember;
local q, idx, recurse, contrib, term, flat,
    res, sols, sol;

    if n = 1 then return 1 fi;

    recurse :=
    proc(l, m, sofar, pos)
    local fact, mult;
        fact := op(1, l[pos]);

        if pos = nops(l)  then
            if m mod fact = 0 then
                contrib := contrib +
                sofar * T(m/fact, k);
            fi;

            return;
        fi;

        for mult to floor(m/fact) do
            if m-mult*fact > 0 then
                recurse(l, m-mult*fact,
                        sofar * T(mult, k),
                        pos + 1);
            fi;
        od;
    end;
    
    res := T(n-1, k);

    for q from 2 to k do
        idx := pet_cycleind_symm(q);

        for term in idx do
            flat := pet_flatten_term(term);

            contrib := 0;
            recurse(flat[2], n-1, 1, 1);
            res := res + flat[1]*contrib;
        od;
    od;

    res;
end;