function GetCoeff(n)
  if n le 1 then return 1; end if;
  G:= Sym(n);
  Kanten:= GSet(G, Subsets({1..n}, 2));
  m:= 0;
  for c in Classes(G) do
    O:= Orbits(sub< G | c[3]>, Kanten); 
    // Die Teilmengen von Kanten, welche fix unter c[3] sind, 
    // sind genau die Vereinigungen von Bahnen in O.
    m +:= c[2] * 2^#O;
  end for;
  return m/#G;
end function;

P<x>:= PolynomialRing(Integers());
P ! [ GetCoeff(n): n in [0..10] ];
// Das liefert
// 12005168*x^10 + 274668*x^9 + 12346*x^8 + 1044*x^7 + 156*x^6 + 34*x^5 + 11*x^4 + 4*x^3 + 2*x^2 + x + 1
