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;