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 = 0 then return 1 fi;
    if k = 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 from 0 to floor(m/fact) do
            recurse(l, m-mult*fact,
                    sofar * T(mult, k),
                    pos + 1);
        od;
    end;

    idx := pet_cycleind_symm(k);
    res := 0;
    
    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;

    res;
end;