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;

cube_edge_cind :=
1/24*(a[1]^12 + 8*a[3]^4 + 6*a[4]^3 + 3*a[2]^6 + 6*a[1]^2*a[2]^5);

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

    if n = 1 then return 1 fi;
    idx_cols := pet_cycleind_symm(n);

    res := 0;

    for term_a in cube_edge_cind do
        for term_b in idx_cols 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_b) do
                    len_b := op(1, v_b);
                    inst_b := degree(term_b, 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)*lcoeff(term_b)*p;
        od;
    od;

    res;
end;