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;

q :=
proc(n)
option remember;
local idx_colors, res, term_a,
    v_a, v_b, inst_a, inst_b, len_a, len_b, p, q;

    if n = 1 then
        idx_colors := [a[1]];
    else
        idx_colors := pet_cycleind_symm(n);
    fi;

    res := 0;

    for term_a in idx_colors do
        p := 1;

        for v_a in indets(term_a) do
            len_a := op(1, v_a);
            inst_a := degree(term_a, v_a);

            q := 0;

            for v_b in indets(term_a) do
                len_b := op(1, v_b);
                inst_b := degree(term_a, v_b);

                if len_a mod len_b = 0 then
                    q := q + len_b*inst_b;
                fi;
            od;

            p := p*q^inst_a;
        od;

        res := res + lcoeff(term_a)*p;
    od;

    res;
end;