Concatenation of analytic functions using Bell polynomials

bell_concat.py

"""
This is the script discussed in the video

https://youtu.be/JCXZw1qC1ZA

$g(x) := 2 + 5 x^2 - 6 x^3 + x^4$

$g(x^2) =$
$= 2 + 5 x^4 - 6 x^6 + x^8$

$g(\pi \cdot x^3) =$
$= 2 + 5 \pi^2 \cdot x^6 - 6 \pi^3 \cdot x^9 + \pi^4 \cdot x^{12}$

$g(3 x - 4 x^2 + x^3) =$
$=\,?$
$=2 + 5\cdot 3^2 \cdot x^2 + \dots ?$

----
Six element set:
$\{A,B,C,D,E,F\}$
Paritions into two sets:
$\{\{A,B,C,D,E\},\{F\}\}, \{\{A,B,C,D,F\},\{E\}\}$, ... six such partitions
$\{\{A,B,C,D\},\{E,F\}\}, \{\{A,B,C,E\},\{D,F\}\}$, ... fifteen such partitions
$\{\{A,B,C\},\{D,E,F\}\}, \{\{A,B,D\},\{C,E,F\}\}$, ... ten such partitions

This gets encoded alla
$B_{6,2}(x_1,x_2,x_3,x_4,x_5)=6x_5x_1+15x_4x_2+10x_3^2$

Then
$B_{n,k}(x_1,x_2,\dots,x_{n-k+1}) = \sum^*_{(j_1,\dots,j_{n-k+1})} n! \cdot\prod_{m=1}^{n-k+1}\large\frac{1}{j_m!}\left(\large\frac{1}{m!}x_m\right)^{j_m}$

where $\sum^*_{(\cdots)}$ is a sum over sequences of $j$'s bound to some arithmetic conditions.
https://en.wikipedia.org/wiki/Bell_polynomials
https://docs.sympy.org/latest/modules/functions/combinatorial.html

----

Given series
$g(x)=\sum_{n>0}g_n\frac{1}{n!}x^n$

$f(x)=\sum_{n>0}f_n\frac{1}{n!}x^n$

Then
$g(f(x)) = \sum_{n>0}c_n\frac{1}{n!}x^n$

with
$c_n = \sum_{k=1}^n g_k\cdot B_{n,k}(f_1,\dots,f_{n-k+1})$

(This is trivial to extend to the case with non-zero $g_0$, but not so much for non-zero $f_0$.)

See also
https://en.wikipedia.org/wiki/Exponential_formula
https://en.wikipedia.org/wiki/Fa%C3%A0_di_Bruno%27s_formula#Formal_power_series_version
"""

from fractions import Fraction
import math
import os
from sympy import bell

# Pretty print bell some polynomials
from sympy import Symbol, symbols
#print(f"bell_6_2 = {bell(6, 2, symbols('x:6')[1:])}\n")
#print(f"bell_4 = {bell(4, Symbol('t'))}\n")


def div(x, y):
    # This implementations is one for f's with integer coefficients for the
    # sake of returning exact ratios. Otherwise, for mere floats, return just x/y.
    assert int(x)==x and int(y)==y
    return Fraction(int(x), (y))

def naturals():
    n = 0
    while True:
        yield n
        n += 1

def monomials(x):
    return (x**n for n in naturals())

def c_times(v, w):
    return (x * y for x, y in zip(v, w))

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

def c_div_fact(f):
    return (div(f_k, math.factorial(k)) for k, f_k in enumerate(f))

def pad_zeros(f, n):
    return f + [0] * (n - len(f))

def inner(v, w):
    return sum(c_times(v, w))

def eval(f, x):
    return inner(f, monomials(x))

def pprint(f):
    def show(c):
        return int(c) if int(c)==c else float(round(c, 3))
    print(f"{list(map(show, f))}\n")


def _bell_circ_exp_coeffs(g, f):
    """
    Inputs g and f are expected in "exp_coeffs" format, i.e. the numbers a_k for the series
    sum_{k>=0} a_k x^k / k!
    (including the 0'th coefficient)
    The return value of this function is also such a sequence.
    """

    # For this implementation, expect f and g are finite polynomials for the sake
    # of just passing finite lists in the examples.
    # Benefit: Can also compute order upfront and yield in finite for-loop.
    g = list(g)
    f = list(f)
    yield g[0]
    assert f[0]==0

    order = 1 + len(g[1:]) * len(f[1:])
    f = pad_zeros(f, order)
    g = pad_zeros(g, order)

    for n in range(1, order):
        bells = (bell(n, k, f[1:(n - k + 1) + 1]) for k in range(1, n + 1))
        yield inner(g[1:], bells)


def bell_circ_coeffs(g, f):
    """
    As opposed to _bell_circ_exp_coeffs, the inputs and return value for this function
    are the numbers a_k for
    sum_{k>=0} a_k x^k
    I.e. without the factorials.
    """
    return c_div_fact(
        _bell_circ_exp_coeffs(
            c_times_fact(g),
            c_times_fact(f)
        )
    )


if __name__=="__main__":
    os.system("clear")

    g = [2, 0, 5, -6]  # 2 + 5 x^2 - 6 x^3
    f = [0, 3, -4, 1]  # 3 x - 4 x^2 + x^3
    gf = list(bell_circ_coeffs(g, f))
    pprint(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_RATIOS = 0
    for x in [div(n, 3) for n in range(100)]:
        test_diff(x, EPS_RATIOS)

    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)

    log = [0] + [(-1)**(n+1) * div(1, n) for n in range(1, 5)]  # log(1+x)
    exp = [0] + [div(1, math.factorial(n)) for n in range(1, 5)]  # exp(x)-1
    exp_log = bell_circ_coeffs(exp, log)
    log_exp = bell_circ_coeffs(log, exp)
    pprint(exp_log)
    pprint(log_exp)