Convolve integers in subquadratic time An linear discrete convolution is an operation that turns two vectors of numbers into a third vector of numbers by multiplying elements inside-out. Formally, for two vectors a and b with elements 0 to n - 1, the discrete linear convolution of a and b can be computed with this loop:

for i = 0 to 2*n - 2
    c[i] = 0;
    for j = 0 to n - 1
        if i - j &gt;= 0 and i - j &lt; n
            c[i] = c[i] + a[j] * b[i - j];


As an example, the convolution of a = { 1, 2, 3, 4 } and b = { 5, 6, 7, 8 } is c = { 5, 16, 34, 60, 61, 52, 32 }.
These convolutions appear when doing long multiplication, for example:

1234 * 5678 =
         20 24 28 32
      15 18 21 24
   10 12 14 16
 5  6  7  8
--------------------
 5 16 34 60 61 52 32
--------------------
      32
     52
    61
   60
  34
 16
 5
--------------------
 7006652


Your task is to write a program or function that, given two arrays (or similar) a and b of non-negative integers of equal length n and optionally, n and an output array c, computes the linear discrete convolution of a and b and returns it, assigns it to the parameter c. or prints it out. You may also take input from the user while or before your code is running. The following constraints apply:


Your program must run in subquadratic or o(n&thinsp;2) time. An algorithm like the pseudo-code above that runs in quadratic time &Theta;(n&thinsp;2) is invalid.
You may assume that all integers in in- and output are in the range from 0 to 65535, this also applies to n.
You may not claim that your algorithm runs in subquadratic time because n has an upper bound.
The results must be exact.
This is code golf, the shortest code in octets wins.
You may not use existing library routines or similar to compute a Fourier transform or a number theoretic transform or the respective inverse transformations.
Python 2, 250 bytes, n^log_2(3)

s=lambda a,b,S=1:[x+y*S for x,y in zip(a,b)]
def f(a,A):
 n=len(a);m=n/2
 if m&lt;1:return[a[0]*A[0]]
 if n%2:return f([0]+a,[0]+A)[2:]
 b,c=a[:m],a[m:];B,C=A[:m],A[m:];p=f(b,B);r=f(c,C);z=m*[0];return s(z+z+r,s(z+s(f(s(b,c),s(B,C)),s(p,r),-1)+z,p+z+z))

This isn't well golfed, but more a proof of concept of the algorithm. Let's look at a more readable version.
s=lambda a,b,S=1:[x+y*S for x,y in zip(a,b)]

def f(a,A):
 n=len(a)
 if n==1:return[a[0]*A[0]]
 if n%2:return f([0]+a,[0]+A)[2:]
 m=n/2
 b,c=a[:m],a[m:]
 B,C=A[:m],A[m:]
 p=f(b,B)
 r=f(c,C)
 q=s(f(s(b,c),s(B,C)),s(p,r),-1)
 z=m*[0]
 return s(z+z+r,s(z+q+z,p+z+z))

The idea is to divide and conquer adapting the Karatsuba multiplication algorithm to do convolution. As before, it runs in time n^log_2(3), as it makes 3 recursive calls which shrink the input length in half. This isn't as good as the quasilinear of DFT, but still better than quadratic.
We can think of convolution as multiplication of formal polynomials. Alternatively, if multiplication were done without regrouping, you'd get convolution. So, adapting the algorithm only requires doing vector addition and subtraction instead of the numerical operations. These additions are still linear-time.
Unfortunately for Python, it doesn't have built-in vector addition and subtraction, so we hack together a function s to do both: addition by default, and subtraction when passed an optional argument of -1. Whenever we add, we make sure to pad with zeroes to avoid truncation.
We handle odd list lengths by padding with a zero to an even length, removing the extraneous zeroes that result.