Multiply numerical polynomials A numerical polynomial is a polynomial \$p\$ in one variable with rational coefficients such that for every integer \$i\$, \$p(i)\$ is also an integer. The numerical polynomials have a basis given by the binomial coefficients:
$$p_n = {x \choose n} = \frac{x(x-1)\cdots(x-n+1)}{n!}$$
For instance:
\$p_0 = 1\$
\$p_1 = x\$
\$p_2 = \frac{x(x-1)}{2} = \frac{1}{2}x^2 - \frac{1}{2}x\$
\$p_3 = \frac{x(x-1)(x-2)}{6} = \frac{1}{6}x^3 - \frac{1}{2}x^2 + \frac{1}{3}x\$
The product of any two numerical polynomials is a numerical polynomial, so there are formulas expressing \$p_m\times p_n\$ as a linear combination of \$p_0, p_1, ..., p_{m+n}\$.
Your job is to produce these formulas.
Goal:
Input: A pair of positive integers \$m\$ and \$n\$
Output: The list of integers \$[a_1,...,a_{m+n}]\$ of length \$m+n\$ such that
$$p_m\times p_n = \sum_{i=1}^{m+n} a_ip_i$$
This is code golf, so shortest code wins.
Examples:
Input: (1,1)
We have \$p_1 = x\$, so \$p_1\times p_1 = x^2\$. The leading term is \$1x^2\$, and the leading term of \$p_2\$ is \$\frac{1}{2!}x^2\$, so we set \$a_2 = \frac{2!}{1} = 2\$. Subtracting off \$2p_2\$ we have \$p_1\times p_1-2p_2 = x^2 - (x^2 - x) = x\$. Thus, we see that \$p_1\times p_1 = p_1 + 2p_2\$, so the output should be \$[1,2]\$.
Input: (1,2)
\$p_2 = \frac{1}{2}x(x-1)\$, so \$p_1\times p_2 = \frac{1}{2}x^2(x-1)\$, which has leading term \$\frac{1}{2}x^3\$. The leading term of \$p_3\$ is \$\frac{1}{3!}x^3\$, so we set \$a_3 = \frac{3!}{2} = 3\$. \$p_1\times p_2 - 3p_3 = x^2-x = 2p_2\$, so we deduce that \$p_1\times p_2=0p_1 + 2p_2 + 3p_3\$, so the output should be \$[0,2,3]\$.
Input (2,2)
The leading term of \$p_2^2\$ is \$\frac{1}{4}x^4\$, so we start with \$p_2^2-\frac{4!}{4}p_4\$. This has leading term \$x^3\$, so we subtract off \$\frac{3!}{1}p_3\$ to get \$p_2^2-\frac{4!}{4}p_4-\frac{3!}{1}p_3\$. This expression turns out to be equal to \$p_2\$, so rearranging we get that \$p_2^2 = 0p_1+p_2+6p_3+6p_4\$, so the output should be \$[0,1,6,6]\$.
Test Cases:
(1,1) ==&gt; [1,2]
(1,2) ==&gt; [0,2,3]
(1,3) ==&gt; [0, 0, 3, 4]
(1,4) ==&gt; [0, 0, 0, 4, 5]
(2,2) ==&gt; [0, 1, 6, 6]
(2,3) ==&gt; [0, 0, 3, 12, 10]
(2,4) ==&gt; [0, 0, 0, 6, 20, 15]
(3,4) ==&gt; [0, 0, 0, 4, 30, 60, 35]
(4,4) ==&gt; [0, 0, 0, 1, 20, 90, 140, 70]
Haskell, 74 bytes



a%b|let a?i|a&lt;1=0^(i-b)^2|d&lt;-a-1=div(i*d?(i-1)-(d-i)*d?i)a=map(a?)[1..a+b]


Try it online!

Uses a recursive formulation I found. Let \$g_b(a,i)\$ be the coefficient of \$p_i\$ in \$p_a p_b\$, which is implemented in the code by a?i (with b fixed). We can express \$g_b\$ recursively as

$$g_b(a,i)=\frac{i\cdot g_b(a-1,i-1)-(a-1-i)\cdot g_b(a-1,i-1)}{a}$$

with base case
$$g_b(0,i)=
\begin{cases}
1,\text{ if }b=i\\
0,\text{ otherwise}\\
\end{cases}
$$

The base case corresponds to the expansion for \$p_b\$, which has a single nonzero coefficient of \$p_i\$ at \$i=b\$. The code implements this indicator function as 0^(i-b)^2.

The code has guards belonging to the definition of ? in the let binding. It's perhaps easier to read as the following (non-working) code where the definition of a?i expects b as a global variable. Perhaps confusingly, the a in the definition of ? is not always the same as the input to %, since it recurses down.

a?i|a&lt;1=0^(i-b)^2|d&lt;-a-1=div(i*d?(i-1)-(d-i)*d?i)a
a%b=map(a?)[1..a+b]


The main function % lists the coefficients for each i from 1 to a+b. It's possible to code the recursion in % directly by zipping shifted lists, but I did not find a short way to do it especially with the dependence on the position i.