Multiply numerical polynomials

Question

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) ==> [1,2]
(1,2) ==> [0,2,3]
(1,3) ==> [0, 0, 3, 4]
(1,4) ==> [0, 0, 0, 4, 5]
(2,2) ==> [0, 1, 6, 6]
(2,3) ==> [0, 0, 3, 12, 10]
(2,4) ==> [0, 0, 0, 6, 20, 15]
(3,4) ==> [0, 0, 0, 4, 30, 60, 35]
(4,4) ==> [0, 0, 0, 1, 20, 90, 140, 70]

Answer 1

APL (Dyalog Extended), ₂₀ 18 bytes

,((×⌿!\)⌹⊢!⍨\⊢)⍳⍤+

Try it online!

A tacit dyadic function that takes m and n as left and right arguments respectively. Mainly uses the matrix division built-in ⌹ to solve the linear equations:

$$ \begin{bmatrix} \binom{1}{1} & \binom{1}{2} & \cdots & \binom{1}{n+m} \\ \binom{2}{1} & \binom{2}{2} & \cdots & \binom{2}{n+m} \\ \vdots & \vdots & \ddots & \vdots \\ \binom{n+m}{1} & \binom{n+m}{2} & \cdots & \binom{n+m}{n+m} \end{bmatrix} \begin{bmatrix} a_1 \\ a_2 \\ \vdots \\ a_{n+m} \end{bmatrix} = \begin{bmatrix} \binom{1}{n} \binom{1}{m} \\ \binom{2}{n} \binom{2}{m} \\ \vdots \\ \binom{n+m}{n} \binom{n+m}{m} \end{bmatrix} \Leftrightarrow Ba = v $$

When \$ B \$ and \$ v \$ are ready, the answer is simply v⌹B. So the main work is to golf the parts to build \$ B \$ and \$ v \$.

How it works

,((×⌿!\)⌹⊢!⍨\⊢)⍳⍤+  ⍝ Left argument: n, Right argument: m
,(            )⍳⍤+  ⍝ Pass [n m] and [1 .. n+m] to the inner function
         ⊢!⍨\⊢      ⍝ Compute B
  (×⌿!\)            ⍝ Compute v
        ⌹           ⍝ Solve the linear equation

Computing \$ B \$

x!y computes \$ \binom{y}{x} \$.

⊢!⍨\⊢  ⍝ Left: [n m], Right: [1 .. n+m]
⊢   ⊢  ⍝ Use right argument for both sides (L, R)
 !⍨\   ⍝ Outer product by flipped binomial:
       ⍝ For each pair of l∊L and r∊R, compute r!l or lCr

Computing \$ v \$

×⌿!\  ⍝ Left: [n m], Right: [1 .. n+m]
  !\  ⍝ Outer product by binomial function
      ⍝ Result is a matrix of two rows [U V]
×⌿    ⍝ Reduce by multiply; compute U×V element-wise

Answer 2

Haskell, 74 bytes

a%b|let a?i|a<1=0^(i-b)^2|d<-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<1=0^(i-b)^2|d<-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.

Answer 3

Haskell, 84 bytes

n#m=[foldr1(-)[i!k*k!n*k!m|k<-[i,i-1..0]]|i<-[1..n+m]]
_!0=1
n!k=n*(n-1)!(k-1)`div`k

Try it online!

Explanation

Given a polynomial \$q\$, the following algorithm (see here) computes the unique coefficients \$a_0,a_1,\ldots\$ such that \$q=a_0p_0+a_1p_1+\ldots\$.

Let $$ \Delta^{(0)}q(x)=q(x)\\ \Delta^{(i+1)}q(x)=\Delta^{(i)}q(x+1)-\Delta^{(i)}q(x). $$

Then \$a_i=\Delta^{(i)}q(0)\$.

Moreover, it can be easily proved that $$ \Delta^{(i)}q(0)=\sum_{k=0}^{i}(-1)^{i-k}{i\choose k}q(k). $$

---

_!0=1
n!k=n*(n-1)!(k-1)`div`k

n!k is the binomial coefficient \$n\choose k\$.

---

n#m=[foldr1(-)[i!k*k!n*k!m|k<-[i,i-1..0]]|i<-[1..n+m]]

The function # computes the answer to the challenge. It simply uses the observation above to compute $$ a_i=\sum_{k=0}^{i}(-1)^{i-k}{i\choose k}{k\choose n}{k\choose m}. $$ The foldr1(-) is the only thing that needs explaining. When applied to a list \$[b_0,b_1,\ldots,b_l]\$, foldr1(-) returns $$ b_0-(b_1-(\ldots-(b_{l-1}-b_l)\ldots)=\sum_{j=0}^{l}(-1)^jb_j $$

---

Haskell, 108 100 bytes

I'm also keeping this one. It is, in my opinion, more elegant, albeit much more verbose.

n#m=take(n+m).tail$head<$>iterate(zipWith(-)=<<tail)[x!n*x!m|x<-[0..]]
_!0=1
n!k=n*(n-1)!(k-1)`div`k

Try it online!

---

_!0=1
n!k=n*(n-1)!(k-1)`div`k

n!k is the binomial coefficient \$n\choose k\$.

---

n#m=take(n+m).tail$head<$>iterate(zipWith(-)=<<tail)[x!n*x!m|x<-[0..]]

The function # computes the answer to the challenge. Let \$q=p_n\cdot p_m\$.

• [x!n*x!m|x<-[0..]] is the infinite list $$ \left[{0\choose n}{0\choose m},{1\choose n}{1\choose m},{2\choose n}{2\choose m},\ldots\right] $$ i.e. $$ [q(0),q(1),q(2),\ldots]. $$
• Given a list \$[b_0,b_1,b_2,\ldots]\$, the function zipWith(-)=<<tail returns \$[b_1-b_0,b_2-b_1,\ldots]\$. Thus the code head<$>iterate(zipWith(-)=<<tail)[x!n*x!m|x<-[0..]] returns the infinite list $$ [\Delta^{(0)}q(0),\Delta^{(1)}q(0),\Delta^{(2)}q(0),\ldots]=[a_0,a_1,a_2,\ldots]. $$
• Finally with take(n+m).tail we drop \$a_0\$ and we take the first \$n+m\$ coefficients, thereby solving the challenge.

Answer 4

Pari/GP, 46 bytes

m->n->b=binomial;[b(m,i-n)*b(i,m)|i<-[1..m+n]]

Try it online!

I did not expect the formula to be so simple: the \$i\$-th element of the output for \$(m,n)\$ is \${m \choose i-n} {i \choose m}\$.

In fact, it can be proven that \$\sum_{i=0}^{m+n} {m \choose i-n} {i \choose m} {x \choose i} = {x \choose m} {x \choose n}\$. The RHS counts the number of ways to choose two subsets of \$\{1,2,\dots,x\}\$ of size \$m\$ and \$n\$ respectively; while the LHS counts, for each \$i\$, the number of ways to choose a subset of size \$i\$, and two subsets of size \$m\$ and \$n\$ respectively whose union is the chosen subset of size \$i\$.

---

Pari/GP, 62 bytes

Saved one byte thanks to Mr. Xcoder.

m->n->b=matpascal(m+n);([b[i,m+1]*b[i,n+1]|i<-[1..#b]]/b~)[^1]

Sadly, GP does not have a built-in for pointwise product.

Try it online!

Answer 5

Wolfram Language (Mathematica), 37 bytes

b[r=Range@+##,#]b[#,r-#2]&
b=Binomial

Try it online!

-5 bytes thanks to att.

---

Wolfram Language (Mathematica), 42 bytes

Table[b=Binomial;b[#,i-#2]b[i,#],{i,+##}]&

Try it online!

---

Wolfram Language (Mathematica), 57 bytes

Saved one byte thanks to Misha Lavrov.

Inverse[b=Array[Binomial,0{,}+##]].(b[[;;,#]]b[[;;,#2]])&

Try it online!

Answer 6

JavaScript (ES6), 109 bytes

This is really just a port of Delfad0r's great answer, implemented as nested recursive functions.

Takes input as (n)(m).

n=>m=>(g=i=>i?[...g(i-1),(h=k=>k&&h(k-1)-(b=(y,x=k)=>y?x*b(y-1,x-1)/y:i-k&1||-1)(k,i)*b(n)*b(m))(i)]:[])(n+m)

Try it online!

How?

The function \$b(y,x)\$ computes:

$$-(-1)^{i-k}{x\choose y}$$

It's slightly shorter to return \$-(-1)^{i-k}\$ as the last iteration of this function than processing it separately. Because we compute the product of three binomial coefficients, the signs of the first two of them are simply cancelling each other out.

The function \$h()\$ is used to compute:

$$\sum_{k=1}^{i}-b(k,i)\times b(n,k)\times b(m,k)\\ =\sum_{k=1}^{i}{(-1)^{i-k}{i\choose k}{k\choose n}{k\choose m}}$$

And the function \$g()\$ is used to compute each term \$a_i\$ for \$1\le i \le n+m\$.

Answer 7

Maxima, 55 bytes

Golfed version. Try it online!

f(m,n):=makelist(binomial(m,i-n)*binomial(i,m),i,1,m+n)

Ungolfed version. Try it online!

f(m,n) := block(
    local(i),
    makelist(binomial(m,i-n)*binomial(i,m), i, 1, m+n)
)$

for x in [[1, 1], [1, 2], [1, 3], [1, 4], [2, 2], [2, 4], [3, 3], [3, 4], [4, 4]] do (
    print(x[1], " -> ", f(x[1], x[2]))
)$

Answer 8

Haskell, 61 bytes

n#m=[n!(m+n-i)*i!n|i<-[1..n+m]]
_!0=1
n!k=n*(n-1)!(k-1)`div`k

Try it online!

Based on @Delfad0r's answer, but using the formula in my PARI/GP answer.