How to get formula of matrix product from formulas of matrices?

I have formulas defining matrices.
The result I want is the formula defining their matrix product.

At no point do I want actual matrices. Those shown below are just for illustration.

The examples I use are Pascal's triangle A and its matrix product with itself B.
(The element-wise product C is also shown, but not what I want.)

I want to calculate sym_b from sym_a:
sym_b = some_magic(sym_a, sym_a)

In the code below I have written sym_b myself.

import sympy as sp


def formula_to_matrix(formula, size):
    matrix = sp.zeros(size, size)
    for _n in range(size):
        for _k in range(_n + 1):
            matrix[_n, _k] = formula.subs(n, _n).subs(k, _k).doit()
    return matrix


n = sp.Symbol('n', integer=True, nonnegative=True)
k = sp.Symbol('k', integer=True, nonnegative=True)
i = sp.Symbol('i', integer=True, nonnegative=True)


# initial matrix A #########################################################

sym_a = sp.binomial(n, k)

mat_a = sp.Matrix([
    [1, 0, 0, 0, 0],
    [1, 1, 0, 0, 0],
    [1, 2, 1, 0, 0],
    [1, 3, 3, 1, 0],
    [1, 4, 6, 4, 1]
])

assert formula_to_matrix(sym_a, 5) == mat_a


# matrix product B (correct) ###############################################

sym_b = sp.Sum(
    sp.binomial(n, i) * sp.binomial(i, k),
    (i, 0, n)
)

mat_b = sp.Matrix([
    [ 1,  0,  0,  0,  0],
    [ 2,  1,  0,  0,  0],
    [ 4,  4,  1,  0,  0],
    [ 8, 12,  6,  1,  0],
    [16, 32, 24,  8,  1]
])

assert formula_to_matrix(sym_b, 5) == mat_b

assert mat_b == mat_a * mat_a  # For matrices the * is matrix multiplication.


# element-wise product C (wrong) ###########################################

sym_c = sym_a * sym_a  # For formulas the * is element-wise multiplication.

mat_c = sp.Matrix([
    [1,  0,  0,  0,  0],
    [1,  1,  0,  0,  0],
    [1,  4,  1,  0,  0],
    [1,  9,  9,  1,  0],
    [1, 16, 36, 16,  1]
])

assert formula_to_matrix(sym_c, 5) == mat_c

I was hoping, that SymPy provides a function for this purpose.

But apparently I have to write it myself, using the rules of matrix multiplication.

def product_of_matrix_formulas(formula_a, formula_b):
    return sp.Sum(
        formula_a.subs(n, n).subs(k, i) * formula_b.subs(n, i).subs(k, k),
        (i, 0, n)
    )

sym_b_calculated = product_of_matrix_formulas(sym_a, sym_a)

assert sym_b_calculated == sym_b