with(combinat);

pet_cycleind_symm :=
proc(n)
local p, s;
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;


pet_cycleind_coll :=
proc(n)
option remember;
local idx, symm, term, flat, len;

    idx := 0;

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

    for term in symm do
        flat := pet_flatten_term(term);

        len :=
        lcm(seq(q, q in map(cyc->op(1, cyc), flat[2])));

        idx := idx + flat[1]*a[len]^(n!/len);
    od;

    idx;
end;

coll :=
proc(n)
option remember;
local idx_slots, idx_sets, res, a, b,
    flat_a, flat_b, cyc_a, cyc_b, len_a, len_b, p, q;

    if n > 1 then
        idx_slots := pet_cycleind_symm(n);
        idx_sets := pet_cycleind_coll(n);
    else
        idx_slots := [a[1]];
        idx_sets := [a[1]];
    fi;

    res := 0;

    for a in idx_slots do
        flat_a := pet_flatten_term(a);

        for b in idx_sets do
            flat_b := pet_flatten_term(b);

            p := 1;
            for cyc_a in flat_a[2] do
                len_a := op(1, cyc_a);
                q := 0;

                for cyc_b in flat_b[2] do
                    len_b := op(1, cyc_b);

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

                p := p*q;
            od;

            res := res + p*flat_a[1]*flat_b[1];
        od;
    od;

    res;
end;

coll_enum :=
proc(n)
    option remember;
    local iter, orbits, orbit, perms, perm, slist;

    orbits := {};

    iter :=
    proc(all, idx, sofar)
        if nops(sofar) = n then
            orbit := {};

            for perm in all do
                slist :=
                [seq(q = perm[q], q=1..n)];
                orbit :=
                {op(orbit),
                 convert(subs(slist, sofar),
                         `multiset`)};
            od;

            orbits := {op(orbits), orbit};
            return;
        fi;

        if idx > n! then return fi;

        iter(all, idx, [op(sofar), all[idx]]);
        iter(all, idx+1, sofar);
    end;

    iter(permute(n), 1, []);

    nops(orbits);
end;