Consult the following link on how to compute the cycle index $Z(G)$ of
the edge permutation group of $K_4$: MSE
link. (We would not
be adding anything here as we would essentially quote it verbatim.) It
was found that
$$Z(G) = \frac{1}{24}
\left(a_1^6 + 8 a_3^2 + 9 a_1^2 a_2^2 + 6 a_2 a_4\right).$$
As a sanity check let us compute non-isomorphic colorings of $K_4$
using at most $N$ colors. We get using Burnside the count
$$\frac{1}{24}
\left(N^6 + 8 N^2 + 9 N^4 + 6 N^2\right)
= \frac{1}{24}
\left(N^6 + 9 N^4 + 14 N^2\right)$$
which yields the sequence
$$1, 11, 66, 276, 900, 2451, 5831, 12496, 24651, 45475, \ldots $$
which points us to OEIS A063842 where this
sequence is described as the number of multigraphs rather than
colorings, a claim that has to be investigated. (Remark, somewhat
later. This issue has now been fixed.) With this in mind we now
compute (use PET rather than Burnside for generating functions)
$$[z^n] Z(G)\left(\frac{1}{1-z}\right)$$
which has OGF as requested in the problem statement
$$\frac{1}{24}\frac{1}{(1-z)^6} + \frac{1}{3}\frac{1}{(1-z^3)^2}
+ \frac{3}{8}\frac{1}{(1-z)^2}\frac{1}{(1-z^2)^2}
+ \frac{1}{4}\frac{1}{1-z^2}\frac{1}{1-z^4}.$$
The repertoire that goes into PET represents the fact that there is
exactly one multi-edge consisting of $q$ single edges, yielding an OGF
of $\sum_{q\ge 0} z^q = 1/(1-z),$ which includes the off edge with
weight zero (no edge present between pair of vertices).
This yields the sequence
$$1, 3, 6, 11, 18, 32, 48, 75, 111, 160, 224, 313,
420, 562, 738, \ldots$$
which points us to OEIS A003082. This has
a perfect match to the problem definition and it looks like the first
entry needs to be qualified or possibly even corrected.
The values look good, e.g. the six multigraphs with three edges
are: a path, a tree with three leaves, a triangle, a double edge and a
single edge, not connected, a double edge and a single edge attached
at one of the vertices of the double one and a triple edge between two
vertices.
Extracting coefficients is not terribly rewarding here but we may
use the Maple code from the following MSE
link to get the
set of polynomials (period due to roots of unity is
$\mathrm{lcm}(2,3,4) = 12):$
5 4 13 3 2 1309 53
n -> 1/2880 n + 1/192 n + --- n + 7/32 n + ---- n + ---
288 2880 192
5 4 13 3 2 203 13
n -> 1/2880 n + 1/192 n + --- n + 7/32 n + --- n + --
288 360 24
5 4 13 3 2 181 39
n -> 1/2880 n + 1/192 n + --- n + 7/32 n + --- n + --
288 320 64
5 4 13 3 2 203
n -> 1/2880 n + 1/192 n + --- n + 7/32 n + --- n + 2/3
288 360
5 4 13 3 2 1309 53
n -> 1/2880 n + 1/192 n + --- n + 7/32 n + ---- n + ---
288 2880 192
5 4 13 3 2 27
n -> 1/2880 n + 1/192 n + --- n + 7/32 n + -- n + 7/8
288 40
5 4 13 3 2 1309 53
n -> 1/2880 n + 1/192 n + --- n + 7/32 n + ---- n + ---
288 2880 192
5 4 13 3 2 203
n -> 1/2880 n + 1/192 n + --- n + 7/32 n + --- n + 2/3
288 360
5 4 13 3 2 181 39
n -> 1/2880 n + 1/192 n + --- n + 7/32 n + --- n + --
288 320 64
5 4 13 3 2 203 13
n -> 1/2880 n + 1/192 n + --- n + 7/32 n + --- n + --
288 360 24
5 4 13 3 2 1309 53
n -> 1/2880 n + 1/192 n + --- n + 7/32 n + ---- n + ---
288 2880 192
5 4 13 3 2 27
n -> 1/2880 n + 1/192 n + --- n + 7/32 n + -- n + 1
288 40
The Maple code for this goes as follows:
HZ := 1/24*1/(1-z)^6+1/3*1/(1-z^3)^2
+3/8*1/(1-z)^2*1/(1-z^2)^2+1/4*1/(1-z^2)*1/(1-z^4);
PSEQ :=
proc()
option remember;
local n, lambda, l, res, locs, vals, cfs;
res := [];
lambda := lcm(2,3,4);
cfs := series(HZ, z=0, 6*lambda+2);
for l from 1 to lambda do
locs := [seq(l+p*lambda, p=0..5)];
vals := map(loc -> coeff(cfs, z, loc), locs);
res :=
[op(res),
unapply(interp(locs, vals, n), n)];
od;
res;
end;
X :=
proc(n)
local Fseq, lambda, idx;
Fseq := PSEQ();
lambda := nops(Fseq);
idx := n mod lambda;
if idx = 0 then idx := idx + lambda; fi;
op(idx, Fseq)(n);
end;
Addendum. Some additional material which may perhaps inspire
further exploration of these types of problems. Here is a somewhat
optimized routine for computing the cycle indices of the edge
permutation group of the complete graph $K_n.$ It is about twice as
fast as the routine posted at the following MSE
link, which includes
some explanatory material. We obtain e.g. the following cycle index
for $K_6:$
$$Z(G_6) =
{\frac {{a_{{1}}}^{15}}{720}}+1/48\,{a_{{1}}}^{7}{a_{{2}}}^{4}+1
/18\,{a_{{1}}}^{3}{a_{{3}}}^{4}+1/12\,{a_{{1}}}^{3}{a_{{2}}}^{6}
\\+1/4\,a_{{1}}a_{{2}}{a_{{4}}}^{3}+1/6\,a_{{1}}a_{{2}}{a_{{3}}}^{
2}a_{{6}}+1/5\,{a_{{5}}}^{3}\\+1/18\,{a_{{3}}}^{5}+1/6\,a_{{3}}{a_
{{6}}}^{2}.$$
This can of course be used to count multigraphs. Working with ordinary
graphs we obtain the generating function
$$Z(G_6)(1+z) =
{z}^{15}+{z}^{14}+2\,{z}^{13}+5\,{z}^{12}+9\,{z}^{11}+15\,{z}^{
10}+21\,{z}^{9}\\+24\,{z}^{8}+24\,{z}^{7}+21\,{z}^{6}+15\,{z}^{5}+
9\,{z}^{4}+5\,{z}^{3}+2\,{z}^{2}+z+1$$
for a total of $156$ graphs. E.g. the nine non-isomorphic graphs on
six vertices with four edges would appear to be, singletons being
omitted: the path on five vertices, a tree with four leaves, a tree
with three leaves, a square, a triangle with a node attached to it, a
path on four nodes and one on two nodes, a triangle and a detached
connected pair, a tree with three leaves and a connected pair, and two
paths on three nodes.
The code follows. Do consult it to clarify the details of the
technique being used, it is included here in place of introducing new
notation to describe the algorithm.
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_edg :=
proc(n)
option remember;
local all, term, res, cycs, l1, l2, inst1, flat, u, v,
uidx, vidx;
if n=0 then return 1; fi;
all := 0:
for term in pet_cycleind_symm(n) do
flat := pet_flatten_termA(term);
cycs := flat[2]; res := 1;
# edges on different cycles of different sizes
for uidx to nops(cycs) do
u := cycs[uidx];
for vidx from uidx+1 to nops(cycs) do
v := cycs[vidx];
l1 := op(1, u); l2 := op(1, v);
res := res *
a[lcm(l1, l2)]^((l1*l2/lcm(l1, l2))*
op(2, u)*op(2, v));
od;
od;
# edges on different cycles of the same size
for uidx to nops(cycs) do
u := cycs[uidx];
l1 := op(1, u); inst1 := op(2, u);
# a[l1]^(1/2*inst1*(inst1-1)*l1*l1/l1)
res := res *
a[l1]^(1/2*inst1*(inst1-1)*l1);
od;
# edges on identical cycles of some size
for uidx to nops(cycs) do
u := cycs[uidx];
l1 := op(1, u); inst1 := op(2, u);
if type(l1, odd) then
# a[l1]^(1/2*l1*(l1-1)/l1);
res := res *
(a[l1]^(1/2*(l1-1)))^inst1;
else
# a[l1/2]^(l1/2/(l1/2))*a[l1]^(1/2*l1*(l1-2)/l1)
res := res *
(a[l1/2]*a[l1]^(1/2*(l1-2)))^inst1;
fi;
od;
all := all + flat[1]*res;
od;
all;
end;
pet_varinto_cind :=
proc(poly, ind)
local subs1, subs2, polyvars, indvars, v, pot, res;
res := ind;
polyvars := indets(poly);
indvars := indets(ind);
for v in indvars do
pot := op(1, v);
subs1 :=
[seq(polyvars[k]=polyvars[k]^pot,
k=1..nops(polyvars))];
subs2 := [v=subs(subs1, poly)];
res := subs(subs2, res);
od;
res;
end;
VGF :=
proc(n)
option remember;
expand(pet_varinto_cind(1+z, pet_cycleind_edg(n)));
end;
Remark. I only just now noticed that the OP was asking for a hint
rather than a solution. It is hoped that there are enough details to
be filled in here for it to be instructive to assemble a comprehensive
answer.
Addendum Nov 18 2018. It is possible to simplify the routine for
the cycle index of the edge permutation group even further. It is in
fact not necessary to flatten the terms of the cycle index of the
symmetric group. It is perfectly sufficient to use the Maple functions
degree and indets to implement an interface to the monomials
from the cycle index as multisets. This is shown below (duplicate code
omitted). The sequence here is
$$1, 2, 4, 11, 34, 156, 1044, 12346, 274668, 12005168,
\\ 1018997864, 165091172592, 50502031367952,
\\ 29054155657235488, 31426485969804308768,
\\ 64001015704527557894928, \ldots $$
which points us to OEIS A000088.
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_cycleind_edg :=
proc(n)
option remember;
local all, term, termvars, res, l1, l2, inst1, u, v,
uidx, vidx;
if n=0 then return 1; fi;
all := 0:
for term in pet_cycleind_symm(n) do
termvars := indets(term); res := 1;
# edges on different cycles of different sizes
for uidx to nops(termvars) do
u := op(uidx, termvars);
l1 := op(1, u);
for vidx from uidx+1 to nops(termvars) do
v := op(vidx, termvars);
l2 := op(1, v);
res := res *
a[lcm(l1, l2)]
^((l1*l2/lcm(l1, l2))*
degree(term, u)*degree(term, v));
od;
od;
# edges on different cycles of the same size
for u in termvars do
l1 := op(1, u); inst1 := degree(term, u);
# a[l1]^(1/2*inst1*(inst1-1)*l1*l1/l1)
res := res *
a[l1]^(1/2*inst1*(inst1-1)*l1);
od;
# edges on identical cycles of some size
for u in termvars do
l1 := op(1, u); inst1 := degree(term, u);
if type(l1, odd) then
# a[l1]^(1/2*l1*(l1-1)/l1);
res := res *
(a[l1]^(1/2*(l1-1)))^inst1;
else
# a[l1/2]^(l1/2/(l1/2))*a[l1]^(1/2*l1*(l1-2)/l1)
res := res *
(a[l1/2]*a[l1]^(1/2*(l1-2)))^inst1;
fi;
od;
all := all + lcoeff(term)*res;
od;
all;
end;
