cp's OEIS Frontend

Sage

Pol. = InfinitePolynomialRing(QQ)
@cached_function
def Z(n):
    if n == 0: return Pol.one()
    return sum(x[k]*Z(n-k) for k in (1..n))/n
@cached_function
def monprod(M):
    p = Pol.one()
    V = [m.variables() for m in M]
    T = cartesian_product(V)
    for t in T:
        r = [Pol.varname_key(str(u))[1] for u in t]
        j = [Pol(M[u[0]]).degree(u[1]) for u in enumerate(t)]
        lcm_r = lcm(r)
        p *= x[lcm_r]^(prod(r)/lcm_r*prod(j))
    return p
@cached_function
def pol_isotop(n,k):
    P = Z(n)
    p = Pol.zero()
    coeffs = P.coefficients()
    mons = P.monomials()
    C = cartesian_product(k*[mons])
    Csorted = [tuple(sorted(u)) for u in C]
    Cset = set(Csorted)
    for c in Cset:
        p += Csorted.count(c)*prod([coeffs[mons.index(u)] for u in c])*monprod(c)
    return p
@cached_function
def rule_sub(r,m):
    D = 0
    for d in divisors(r):
        try: D += d*m.degrees()[-d-1]
        except: break
    return D
def a(n,k=3):
    P = Z(n)
    coeffs = P.coefficients()
    Q = pol_isotop(n,k)
    inds = [Pol.varname_key(str(u))[1] for u in Q.variables()]
    p = 0
    for mon in enumerate(P.monomials()):
        m = Pol(mon[1])
        p += coeffs[mon[0]]*Q.subs({x[i]:rule_sub(i,m) for i in inds})
    return p
# Philip Turecek, Jun 17 2023