Combinatorial counting with symmetry Let $A$ be a set of objects where $|A|=n$. We want to count all the possible ways that we can arrange these objects into $n$ bags with exactly $n$ objects in each. We can reuse any object, however, no repetition is allowed inside the bags. 
With $A=\{a,b,c\}$, for example, $[(a,b,c), (a,b,c), (b,c,a)]$ is a valid outcome.
Obviously there are $(n!)^n$ ways to do this. 
Now we want to add two extra constraints:


*

*The order of bags is not important. 


For example, $[(a,b,c), (a,b,c), (b,c,a)]$ would be identical to $[(b,c,a), (a,b,c), (a,b,c)]$.


*

*The label of objects inside the bags do not matter. Only the relative positions are important.


For example, $[(a,b,c), (a,b,c), (b,c,a)]$ would be identical to $[(c,b,a), (c,b,a), (b,a,c)]$ and is identical to $[(a,c,b), (a,c,b), (c,b,a)]$ etc.
Questions are:


*

*How many ways can we set these bags given the above constrains ?

*Is there any algorithm to output all these possible combinations?

 A: This appears  to be  an interesting problem  which has  a surprisingly
simple answer and can be attacked using Power Group Enumeration as
described in quite some detail at the following
MSE link.

I will not  describe the algorithm here as all the  details are in the
cited  post.   In the  present  case we  have  for  the objects  being
distributed into the slots the permutations of $n$ elements. The group
acting on the slots is the  symmetric group on $n$ elements. The group
acting on the $n!$  permutations simultaneously is the symmetric group
whose  constitutent  permutations create  a  permutation  of the  $n!$
permutations by acting on the elements of each.

We may therefore apply the algorithm that I cited and it does not need
to be  modified.  The  only missing  piece is the  cycle index $Z(G)$ of the
action of  the symmetric  group on the  permutations. This is  easy to
compute however  as a permutation  $\gamma$ acting on the  elements of
the  $n!$ permutations that  go into  the slots  partitions everything
into  cycles whose  length is  the  LCM $q$  of the  cycle lengths  of
$\gamma.$ (All elements of the permutation move simultaneously when we
apply $\gamma.$) Think of placing  a marker indicating the position of
a  given element beside the elements of the cycles of $\gamma.$  These markers move  in parallel and  the first  time  all markers  return to  their original position is after $q$ steps.

Implementing this  in Maple yields  the following code. (We  have also
included a  total enumeration  routine to verify  the count  for small
$n.$)

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;

The Power Group Enumeration code gave  the following sequence (as opposed to
total enumeration which is practicable only up to $n=4$):
$$1, 2, 10, 762, 1876255, 274382326290, 3265588553925722827,
\\ 4299566944396584777543664576, 
\\ 828675148077536475804944305151462053905, 
\\ 30068353582978459601855528390398866877243129478172220,
\ldots$$
Remark. The code for the case $n=1$ given above can be simplified.
I leave this to the next version.
