如何从矩阵的元素公式推导矩阵乘积的元素公式？

我有定义矩阵元素的公式，想要直接推导它们矩阵乘积的元素公式，全程不需要实际生成矩阵，下方展示的矩阵仅作示例。以帕斯卡三角矩阵A及其自乘的矩阵乘积B为例（同时展示了逐元素乘积C，但这不是我需要的）。我希望能实现一个some_magic函数，让sym_b = some_magic(sym_a, sym_a)，下方代码中我已手动写出了sym_b的表达式。

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)


# 初始矩阵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


# 矩阵乘积B（正确结果） ###############################################

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  # 矩阵的*运算符代表矩阵乘法


# 逐元素乘积C（非需求结果） ###########################################

sym_c = sym_a * sym_a  # 公式的*运算符代表逐元素相乘

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

内容的提问来源于stack exchange，提问作者Watchduck