I want to multiply matrices defined by formulas, and get a result defined by a formula. E.g. Pascal's triangle is defined by the binomial coefficient. Its matrix product with itself is correct_square_matrix. (Sloane's A038207, which is related to hypercubes.) I can calculate correct_square_matrix for a particular size. But what I want is correct_square_formula. The function I am looking for is called some_magic in the code below. 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) return matrix n = sp.Symbol('n', integer=True, nonnegative=True) k = sp.Symbol('k', integer=True, nonnegative=True) # initial matrix ############################################################# pascal_formula = sp.binomial(n, k) pascal_matrix = formula_to_matrix(pascal_formula, 5) assert pascal_matrix == 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] ]) # matrix product ############################################################# correct_square_matrix = 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 pascal_matrix * pascal_matrix == correct_square_matrix # correct_square_formula = some_magic(pascal_formula, pascal_formula) # assert correct_square_matrix == formula_to_matrix(correct_square_formula, 5) # element-wise product ####################################################### wrong_square_matrix = 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] ]) wrong_square_formula = pascal_formula * pascal_formula assert wrong_square_matrix == formula_to_matrix(wrong_square_formula, 5) Continue reading...