Discrete Convolution or Polynomial Multiplication

Question

Given two non empty lists of integers, your submission should calculate and return the discrete convolution of the two. Interestingly, if you consider the list elements as coefficients of polynomials, the convolution of the two lists represents the coefficients of the product of the two polynomials.

Definition

Given the lists A=[a(0),a(1),a(2),...,a(n)] and B=[b(0),b(1),b(2),...,b(m)] (setting a(k)=0 for k<0 and k>n and b(k)=0 for k<0 and k>m) then the convolution of the two is defined as A*B=[c(0),c(1),...,c(m+n)] where c(k) = sum [ a(x)*b(y) for all integers x y such that x+y=k]

Rules

• Any convenient input and output formatting for your language is allowed.
• Built-ins for convolution, creating convolution matrices, correlation and polynomial multiplication are not allowed.

Examples

[1,1]*[1] = [1,1]
[1,1]*[1,1] = [1,2,1]
[1,1]*[1,2,1] = [1,3,3,1]
[1,1]*[1,3,3,1] = [1,4,6,4,1]
[1,1]*[1,4,6,4,1] = [1,5,10,10,5,1]

[1,-1]*[1,1,1,1,1] = [1,0,0,0,0,-1]
[80085,1337]*[-24319,406] = [-1947587115,7,542822]

Answer 1

J, 10 8 bytes

[:+//.*/

Usage:

ppc =: [:+//.*/    NB. polynomial product coefficients 
80085 1337 ppc _24319 406
_1947587115 7 542822

Description: The program takes takes two lists, makes a multiplication table, then adds the numbers on the positive diagonals.

Answer 2

MATL, 19 bytes

PiYdt"TF2&YStpsw]xx

Try it online!

Explanation

This builds a block-diagonal matrix with the two inputs, reversing the first. For example, with inputs [1 4 3 5], [1 3 2] the matrix is

[ 5 3 4 1 0 0 0
  0 0 0 0 1 3 2 ]

Each entry of the convolution is obtained by shifting the first row one position to the right, computing the product of each column, and summing all results.

In principle, the shifting should be done padding with zeros from the left. Equivalently, circular shifting can be used, because the matrix contains zeros at the appropriate entries.

For example, the first result is obtained from the shifted matrix

[ 0 5 3 4 1 0 0
  0 0 0 0 1 3 2 ]

and is thus 1*1 == 1. The second is obtained from

[ 0 0 5 3 4 1 0
  0 0 0 0 1 3 2 ]

and is thus 4*1+1*3 == 7, etc. This must be done m+n-1 times, where m and n are the input lengths. The code uses a loop with m+n iterations (which saves some bytes) and discards the last result.

P          % Take first input (numeric vactor) implicitly and reverse it
i          % Take second input (numeric vactor) 
Yd         % Build diagonal matrix with the two vectors
t          % Duplicate
"          % For each column of the matrix
  TF2&YS   %   Circularly shift first row 1 step to the right
  t        %   Duplicate
  p        %   Product of each column
  s        %   Sum all those products
  w        %   Swap top two elements in stack. The shifted matrix is left on top
]          % End for
xx         % Delete matrix and last result. Implicitly display

Answer 3

Haskell, 55 49 bytes

(a:b)#c=zipWith(+)(0:b#c)$map(a*)c++[]#b
_#c=0<$c

Defines an operator #.

Answer 4

Matlab / Octave, 41 bytes

@(p,q)poly([roots(p);roots(q)])*p(1)*q(1)

This defines an anonymous function. To call it, assign it to a variable or use ans.

Try it here.

Explanation

This exploits the facts that

• The (possibly repeated) roots characterize a polynomial up to its leading coefficient.
• The product of two polynomials has the roots of both.

The code computes the roots of the two polynomials (function roots) and concatenates them into a column array. From that it obtains the coefficients of the product polynomial with a leading 1 (function poly). Finally the result is multiplied by the leading coefficients of the two polynomials.

Answer 5

05AB1E, 18 17 bytes

Code

0Ev²¹g<Å0«y*NFÁ}+

Explanation

The theory behind:

To find the convolution, let's take the example [1, 2, 3], [3, 4, 5]. We position the values of the first array upside down and vertically, like this:

3
2
1

Now, we place the second array like a ladder and multiply it by :

3 ×       [3  4  5]
2 ×    [3  4  5]
1 × [3  4  5]

Resulting into:

        9   12   15
    6   8   10
3   4   5

Then, we add them up, resulting into:

        9   12   15
    6   8   10
3   4   5       

3   10  22  22   15

So, the convolution is [3, 10, 22, 22, 15].

The code itself:

We are going to do this step by step using the [1, 2, 3], [3, 4, 5] as the test case.

0Ev²¹g<Å0«y*NFÁ}+

We first push 0 and then we Evaluate the first input array. We map over each element using v.

So, for each element, we push the second array with ² and then the length of the first array using ¹g and decrease this by 1 (with <). We convert this into a list of zeros with (length 1st array - 1) zeros, using Å0 and append this to our list. Our stack now looks like this for the first item in the input list:

[3, 4, 5, 0, 0]

We multiply this array with the current item, done with y*. After that, we push N, which indicates the index of the current item (zero-indexed) and rotate the array that many times to the right using FÁ}. Finally, we add this to our initial value (0). So, what basically is done is the following:

[0, 0, 9, 12, 15] +
[0, 6, 8, 10, 0] +
[3, 4, 5, 0, 0] =

[3, 10, 22, 22, 15]

Which is then implicitly printed. Uses CP-1252 encoding. Try it online!.

Answer 6

Jelly, 9 bytes

0;+
×'Ṛç/

Try it online! or verify all test cases.

How it works

×'Ṛç/  Main link. Arguments: p, q (lists)

×'     Spawned multiplication; multiply each item of p with each item of q.
  Ṛ    Reverse the rows of the result.
   ç/  Reduce the rows by the helper link.


0;+    Helper link. Arguments: p, q (lists)

0;     Prepend a 0 to p.
  +    Perform vectorized addition of the result and q.

Answer 7

Octave, 48 bytes

@(p,q)ifft(fft([p q*0]).*fft([q p*0]))(1:end-1)

Try it here.

Explanation

Discrete convolution corresponds to multiplication of the (discrete-time) Fourier transforms. So one way to multiply the polynomials would be transform them, multiply the transformed sequences, and transform back.

If the discrete Fourier transform (DFT) is used instead of the Fourier transform, the result is the circular convolution of the original sequences of polynomial coefficients. The code follows this route. To make circular convolution equal to standard convolution, the sequences are zero-padded and the result is trimmed.

Answer 8

Husk, 5 bytes

mΣ∂Ṫ*

Try it online!

Note: When supplying the zero-polynomial/empty list, you need to specify its type (ie. []:LN)!

Explanation

mΣ∂Ṫ*  -- implicit inputs xs ys, for example: [1,-1] [1,1]
   Ṫ*  -- compute the outer product xsᵀ·ys: [[1,1],[-1,-1]]
  ∂    -- diagonals: [[1],[1,-1],[-1]]
mΣ     -- map sum: [1,0,1]

Answer 9

Matlab, 33 bytes

@(x,y)sum(spdiags(flip(x').*y),1)

Try it online!

Creates a matrix of all element-wise products of the inputs, then sums along the diagonals. The ,1 at the end forces matlab to sum along the correct direction when one of the input vectors has length 1.

In Octave spdiags does not work for vectors, resulting in an error when one of the inputs has length 1. Matlab 2016b or newer is required for the explicit expansion of the element-wise product.

Answer 10

Itr, 6 bytes

##C*µS

reads two lists from standard input, returns a list

online interpreter

Explanation

##     ; read two lists from standard input
  C    ; go through all diagonals of the Cartesian product
   *   ; multiplying the elements from the two lists
    µS ; sum up the diagonals
       ; implicit output

Answer 11

Ruby, 83 bytes

Almost directly pulled out of an answer I did earlier regarding expanding roots into a polynomial.

->a,b{q=a.map{|c|r=b.map{|e|e*c};b=[0]+b;r}
q.pop.zip(*q).map{|e|(e-[p]).reduce:+}}

Answer 12

Python, 90 bytes

lambda p,q:[sum((p+k*[0])[i]*(q+k*[0])[k-i]for i in range(k+1))for k in range(len(p+q)-1)]

Answer 13

JavaScript (ES6), 64 bytes

(a,b)=>a.map((n,i)=>b.map((m,j)=>r[j+=i]=m*n+(r[j]||0)),r=[])&&r

Returns the empty array if either input is empty. Based on my answer to Polynomialception.

Answer 14

Julia, 62 55 bytes

~=endof
p%q=[sum(diag(p*rot180(q'),-k))for k=1-~q:~p-1]

Try it online!

Answer 15

Clojure, 104 bytes

#(vals(apply merge-with +(sorted-map)(for[i(range(count %))j(range(count %2))]{(+ i j)(*(% i)(%2 j))})))

Merging to sorted-map guarantees that values are returned in the correct order. I wish there were a few more test cases.