Concatenation of series via Bell polynomials

bell_circ.py

import math

from sympy import bell  # https://en.wikipedia.org/wiki/Bell_polynomials
"""
# Pretty print bell polynomials
from sympy import Symbol, symbols
print(bell(4, Symbol('t')))
print(bell(6, 2, symbols('x:6')[1:]))
print()
"""


########## Utility functions
def circ(f, g):
    return lambda x: f(g(x))

def eval(f, x):
    # f is understood here as polynomial represented by its coefficient list
    return sum(f_k * x**k for k, f_k in enumerate(f))

def to_exp(f):
    return [f_k / math.factorial(k) for k, f_k in enumerate(f)]

def from_exp(f):
    return [f_k * math.factorial(k) for k, f_k in enumerate(f)]


def _from_exp_jointly_padded(f, g):
    # Map to exp form and pad with zeros
    n = len(g) * len(f)  # Should be at least as big as order of concatenated functions
    def from_exp_padded(fun):
        return from_exp(fun) + (-len(fun) + n) * [0]
    return map(from_exp_padded, (f, g))


def bell_circ(g, f):
    # https://en.wikipedia.org/wiki/Fa%C3%A0_di_Bruno%27s_formula#Formal_power_series_version
    f, g = _from_exp_jointly_padded(f, g)

    def coeff(n):
        def t(k):
            idx = n - k + 1
            arg = f[1:idx+1]
            return g[k] * bell(n, k, arg)
        return sum(map(t, range(1, n+1)))

    res = to_exp(map(coeff, range(len(f))))
    res[0] = g[0]
    return res


if __name__=="__main__":
    g = [2, 0, 5, -6]  # 2 + 5 x^2 - 6 x^3 + x^4
    f = [0, 3, -4, 1]  # 3 x - 4 x^2 + x^3
    gf = bell_circ(g, f)
    print(gf)

    # Test
    def test_diff(x, eps):
        lhs = eval(g, eval(f, x))
        rhs = eval(gf, x)
        diff = abs(lhs-rhs)
        assert diff <= eps, f"diff={diff}"

    EPS_INTS = 0  # Note: Assumes also polynomials have float coeffs
    for x in [n for n in range(100)]:
        test_diff(x, EPS_INTS)

    EPS_FLOATS = 0.05  # Bigger epsilon needed for longer polynomials
    for x in [n/math.pi for n in range(100)]:
        test_diff(x, EPS_FLOATS)