Revision history [back]

First, Sage can compute cycle index out of the box, although it's expressed in terms of symmetric polynomials and so some conversion is needed - for example,

def Z(n):
    ind = SymmetricGroup(n).cycle_index()
    R.<a> = InfinitePolynomialRing(QQ)
    return sum(prod(a[i] for i in t)*c for t,c in ind)

Second, if one wants to implement a recurrence formula, it's worth to notice that we deal with polynomials here and employ the corresponding algebraic structure. Also, to avoid recomputing the same result over and over again, it's worth to use @cached_function decorator:

@cached_function
def Z(n):
    R.<a> = InfinitePolynomialRing(QQ)
    if n==0:
        return R.one()
    return sum( a[k] * Z(n-k) for k in (1..n) ) / n