# The following is a cut and paste of the code listing appearing in Appendix B of the Andrew Gainer-Dewar reference.
# Γ-species and the enumeration of k-trees, Electron. J. Combin. 19, no. 4, (2012), P45.
# doi: https://doi.org/10.37236/2615

# Minor edits have been made for the code to work on a current version of Sage.

psr = PowerSeriesRing(QQ , 'x')
x = psr.gen()
# Compute the generating function for unlabeled Y-rooted k-trees fixed by permutations of a given cycle type mu.
# Note that mu should partition k
@cached_function
def unlY(mu , n):
  if n <= 0:
    return psr (1)
  else:
    ystretcher = lambda c, part:unlY(Partition(part.power(c)),floor((n-1)/c)).subs({x:x**c})
    descendant_pseries = lambda part: prod(ystretcher(c, part) for c in part)
    return sum(x**i/i * descendant_pseries(mu.power(i)).subs({x:x**i}) for i in range(1, n+1)).exp(n+1)

@cached_function
def unlXY(mu , n):
  if n <= 0:
    return psr (0)
  else:
    ystretcher = lambda c: unlY(Partition(mu.power(c)[: -1]), floor ((n-1)/c)).subs({x:x**c})
    return (x * prod(ystretcher(c) for c in mu)).add_bigoh(n+1)

# Compute the generating functions for unlabeled X-, Y-, and XY-rooted k-trees using quotients
ax = lambda k, n: sum(1/mu.aut() * unlXY(mu , n) for mu in Partitions(k+1))
ay = lambda k, n: sum(1/mu.aut() * unlY(mu , n) for mu in Partitions(k))
axy = lambda k, n: sum(1/mu.aut() * unlXY(Partition(mu + [1]), n) for mu in Partitions(k))

# Compute the generating function for unlabeled un-rooted k-trees using the dissymmetry theorem
a = lambda k, n: ax(k, n) + ay(k, n) - axy(k, n)

# Print the result
# ALERT: User must substitue values for k and n (number  of hedra)
print( a(3 , 15) )