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;