# 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) ==> [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]

• Just an FYI we have MathJax now so you can use $ to write in-line equations like$ p_{m} $instead of using images. – FryAmTheEggman Oct 13 '18 at 22:55 • Looks like every test is missing the$a_0=0\$entry; so are the inputs actually going to be strictly positive integers or should (0,0) yield an empty list? – Jonathan Allan Oct 13 '18 at 23:08 • @JonathanAllan Strictly positive -- none of the polynomials have a constant coefficient so it would be kind of boring. – Hood Oct 13 '18 at 23:15 • @FryAmTheEggman That's good to know because making the images is quite annoying. – Hood Oct 13 '18 at 23:16 • @JonathanAllan I guess I shouldn't have given p_0 as an example of a basis element... – Hood Oct 13 '18 at 23:16 ## 6 Answers # 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. # APL (Dyalog Extended), 20 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  # 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)divk  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)divk  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)divk  Try it online! _!0=1 n!k=n*(n-1)!(k-1)divk  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.
• Cool! I guess this is a sort of discrete Taylor expansion. – Hood Oct 13 '18 at 23:50
• @Hood Indeed ;) – Delfad0r Oct 13 '18 at 23:58
• This description also explains why p_m*p_n never involves binomial coefficients lower than min(m,n): the discrete derivative has a product rule: D(f*g)=D(f(x))*g(x+1) + f(x)*D(g(x)). Applying this k times, you can see that D^k(f*g)=D^n(f(x))*g(x+k) + D^{k-1}(f(x))*D(g(x+k-1)) + ... + f(x)*D^n(g(x)). If k<n then all of the derivatives of p_n have value zero at zero, so D^k(p_n*g)(0)=0. – Hood Oct 14 '18 at 15:51

# Wolfram Language (Mathematica), 57 bytes

Saved one byte thanks to Misha Lavrov.

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


Try it online!

• We can shave off a byte by replacing +##{1,1} with 0{,}+##, as in this golfing tip. – Misha Lavrov Oct 15 '18 at 1:57

# Pari/GP, 65 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~)[2..#b]


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

Try it online!

• Nice answer. As a completely minor golf, curryinng gives -1 bytes (m->n->). – Mr. Xcoder Oct 14 '18 at 18:49

# 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\$$.