Multiply numerical polynomials
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.