This interesting problem may be solved by Power Group Enumeration.
The challenging part is that the sequence is almost impossible to
compute beyond the initial five terms by total enumeration. The
computation is done by a particularly simple form of PGE as presented
by e.g. Harary and Palmer and (in a different publication)
Fripertinger. The present case does not involve generating functions.
and can in fact be done with Burnside's lemma. The scenario is that we
have two groups, one permuting the slots and another the elements
being placed into these slots. The easy part is that $n$ letters go
into the slots and these letters have the symmetric group $S_n$ acting
on them. What is somewhat more challenging is the cycle index of the
group permuting the slots but only moderately so. We implement the
concept of a multiset by partitioning the $n$ slots according to one
of the $p(n)$ partitions of $n.$ (These are not combinations. We use
combinations in the total enumeration routine but not here.) The
partition corresponds to a multiset of sequences where the lengths of
the sequences are given by the constitutents of the partition.
E.g. the partition $2+2+3$ of $n=7$ corresponds to two sequences of
length two and one of length three. It is now easy to describe the
action of the group acting on the slots (there is one such group for
each partition): it permutes the $q$ constituents of the partition of
the same length $m$ aocording to the symmetric group $S_q$. This
means it must move these slots in parallel since the sequences form
inseparable blocks. Therefore a cycle from $S_q$ must have its
variables $a_p$ replaced by $a_p^m$ as there is one copy of $a_p$
moving the first elements of the blocks of equal length, another for
the second elements and so on until the copy of $a_p$ that moves
elements at position $m.$ These movements must happen
simultaneously. We then obtain the desired cycle index corresponding
to a given partition by combining the action on sequences of the same
length for all lengths that appear (multiplication of cycle
indices). This effectively implements the concept of a multiset of
sequences. We can then compute the number of configurations by
Burnside's lemma which says to average the number of assignments fixed
by the elements of the power group. But this number is easy to
compute. Suppose we have a permutation $\alpha$ from the slot
permutation group $Q$ and a permutation $\beta$ from $S_n.$ If we
place the appropriate number of complete, directed and consecutive
copies of a cycle from $\beta$ on a cycle from $\alpha$ then this
assignment is fixed under the power group action for $(\alpha,\beta)$,
and this is possible iff the length of the cycle from $\beta$ divides
the length of the cycle from $\alpha$. The process yields as many
assignments as the length of the cycle from $\beta.$ In what follows
we have implemented one extra optimization which was absent from
earlier contributions like those in the links shown at the end. This
is that instead of flattening $\alpha$ and $\beta$ into lists of
single cycles and iterating over these singletons to discover the
number of coverings of a cycle from $\alpha$ by cycles from $\beta$ we
have made use of the fact that the existence of a covering only
depends on the lengths of the two cycles and hence we may process sets
of cycles having the same length all at once, a significant
savings. This is achieved by multiplying the number of coverings by
the number of cycles of the same length in $\beta$ as the one used in
the covering and raising the number of coverings to the appropriate
power that represents the number of cycles of the current length in
$\alpha$, which reflects the fact that we may freely combine any
admissible covering of a cycle from $\alpha$ with any other.
We now share the results and the program that was used to compute
them. The sequence goes
$$1, 4, 13, 57, 252, 1322, 7361, 45057, 294392, 2054394, 15172872,
\\ 118175823, 966300054, 8268640847, 73825951226,
\\ 686049132714, 6620780612228, \ldots$$
This is the PGE code that was used.
with(combinat);
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_flatten_termA :=
proc(varp)
local terml, d, cf, v;
terml := [];
cf := varp;
for v in indets(varp) do
d := degree(varp, v);
terml := [op(terml), [op(1,v), d]];
cf := cf/v^d;
od;
[cf, terml];
end;
pet_cycleind_mset :=
proc(msdata)
local msrep, res, cycs, num, slist;
msrep := convert(msdata, `multiset`);
res := 1;
for cycs in msrep do
num := op(2, cycs);
slist := [seq(a[q]=a[q]^op(1, cycs), q=1..num)];
res := res
*subs(slist, pet_cycleind_symm(num));
od;
expand(res);
end;
set_seq :=
proc(n)
option remember;
local part, idx_slots, idx_syms, res, a, b,
flat_a, flat_b, cyc_a, cyc_b, len_a, len_b, p, q;
if n > 1 then
idx_syms := pet_cycleind_symm(n);
else
idx_syms := [a[1]];
fi;
res := 0;
for part in partition(n) do
idx_slots := pet_cycleind_mset(part);
if not type(idx_slots, `+`) then
idx_slots := [idx_slots];
fi;
for a in idx_slots do
flat_a := pet_flatten_termA(a);
for b in idx_syms do
flat_b := pet_flatten_termA(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*op(2, cyc_b);
fi;
od;
p := p*q^op(2, cyc_a);
od;
res := res + p*flat_a[1]*flat_b[1];
od;
od;
od;
res;
end;
set_seq_enum :=
proc(n)
option remember;
local orbits, orbit, ind, slist, perm, part, comb, item,
entlst, sofar, sb, len;
if n=1 then return 1 fi;
slist := [];
for perm in permute(n) do
slist :=
[op(slist), [seq(q-1=perm[q]-1, q=1..n)]];
od;
orbits := table();
for ind from n^n to 2*n^n-1 do
item := convert(ind, base, n)[1..n];
for part in partition(n) do
for comb in permute(part) do
entlst := []; sofar := 0;
for len in comb do
entlst :=
[op(entlst), item[1+sofar..len+sofar]];
sofar := sofar + len;
od;
orbit := [];
for sb in slist do
orbit :=
[op(orbit), sort(subs(sb, entlst))];
od;
orbits[convert(orbit, `set`)] := 1;
od;
od;
od;
# nops([indices(orbits, `nolist`)]);
numelems(orbits);
end;
PGE recently appeared at this MSE link
I and this MSE
link II.
