with(combinat); with(numtheory); 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_cycleind_cards := proc() option remember; subs([seq(a[q]=a[q]^13, q=1..4)], pet_cycleind_symm(4)); end; q := proc(N) option remember; local idx_slots, idx_cards, res, term_a, term_b, v_a, inst_a, inst_b, len_a, p; if N = 0 then return 1 fi; if N = 1 then idx_slots := [a[1]]; else idx_slots := pet_cycleind_symm(N); fi; idx_cards := pet_cycleind_cards(); res := 0; for term_a in idx_slots do for term_b in idx_cards do p := 1; for v_a in indets(term_a) do len_a := op(1, v_a); inst_a := degree(term_a, v_a); inst_b := degree(term_b, v_a); if inst_b >= inst_a then p := p*binomial(inst_b, inst_a) *inst_a!*len_a^inst_a; else p := 0; break; fi; od; res := res + lcoeff(term_a)*lcoeff(term_b)*p; od; od; res; end;