19
\$\begingroup\$

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]
\$\endgroup\$
  • 3
    \$\begingroup\$ The specification implies that inputs of length n, m should produce an output of length n + m − 1, but this does not hold for your test case [1,1]*[] = [], and cannot possibly hold for []*[] = ?. Convolution is not well-defined on empty lists. I think you should guarantee that the input lists are nonempty. \$\endgroup\$ – Anders Kaseorg May 16 '16 at 20:55
  • 1
    \$\begingroup\$ @AndersKaseorg You are right, I'll change that. \$\endgroup\$ – flawr May 16 '16 at 21:13
  • \$\begingroup\$ related: codegolf.stackexchange.com/q/67289/29325 \$\endgroup\$ – sergiol Mar 9 '17 at 2:39

14 Answers 14

14
\$\begingroup\$

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.

\$\endgroup\$
  • \$\begingroup\$ Very clever approach! \$\endgroup\$ – Luis Mendo May 16 '16 at 23:46
  • \$\begingroup\$ You don't need to count the parentheses. The expression inside them evaluates to a tacit verb, which can be assigned to a variable. \$\endgroup\$ – Dennis May 17 '16 at 0:11
  • \$\begingroup\$ Great example of adverbs! \$\endgroup\$ – miles May 17 '16 at 0:26
6
\$\begingroup\$

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
\$\endgroup\$
4
\$\begingroup\$

Haskell, 55 49 bytes

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

Defines an operator #.

\$\endgroup\$
  • 1
    \$\begingroup\$ I think the padding [0,0..] can be (0<$b) to give exactly the needed length, allowing the empty base case _#b=0<$b. \$\endgroup\$ – xnor May 17 '16 at 9:58
  • \$\begingroup\$ @xnor Indeed, that saves 6 bytes. \$\endgroup\$ – Anders Kaseorg May 17 '16 at 19:35
  • \$\begingroup\$ Now that I finally understand your answer, I have to say that this is just so damn clever! I'm impressed! \$\endgroup\$ – flawr Jan 1 '17 at 19:09
3
\$\begingroup\$

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.

\$\endgroup\$
3
\$\begingroup\$

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.

\$\endgroup\$
  • \$\begingroup\$ Damnit, I still sould have banned fft, but good job! \$\endgroup\$ – flawr May 16 '16 at 21:33
  • \$\begingroup\$ @flawr Yes, I think we talked about that...? :-P \$\endgroup\$ – Luis Mendo May 16 '16 at 21:34
2
\$\begingroup\$

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!.

\$\endgroup\$
2
\$\begingroup\$

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.
\$\endgroup\$
  • \$\begingroup\$ What‽ Jelly longer than J‽ That's impossible by definition! \$\endgroup\$ – Adám May 17 '16 at 21:53
2
\$\begingroup\$

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]
\$\endgroup\$
2
\$\begingroup\$

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.

\$\endgroup\$
  • \$\begingroup\$ Nice approach!! \$\endgroup\$ – Luis Mendo Dec 21 '17 at 22:02
1
\$\begingroup\$

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:+}}
\$\endgroup\$
1
\$\begingroup\$

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)]
\$\endgroup\$
1
\$\begingroup\$

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.

\$\endgroup\$
1
\$\begingroup\$

Julia, 62 55 bytes

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

Try it online!

\$\endgroup\$
1
\$\begingroup\$

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.

\$\endgroup\$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.