Module karatsuba_convolution

Calculate the convolution of two vectors using a Karatsuba-style decomposition and the CRT.

This is not a new idea, but we did have the pleasure of reinventing it independently. Some references:

• https://cr.yp.to/lineartime/multapps-20080515.pdf
• https://2π.com/23/convolution/

Given a vector v \in F^N, let v(x) \in F[X] denote the polynomial v_0 + v_1 x + … + v_{N - 1} x^{N - 1}. Then w is equal to the convolution v * u if and only if w(x) = v(x)u(x) mod x^N - 1. Additionally, define the negacyclic convolution by w(x) = v(x)u(x) mod x^N + 1. Using the Chinese remainder theorem we can compute w(x) as w(x) = 1/2 (w_0(x) + w_1(x)) + x^{N/2}/2 (w_0(x) - w_1(x)) where w_0 = v(x)u(x) mod x^{N/2} - 1 w_1 = v(x)u(x) mod x^{N/2} + 1

To compute w_0 and w_1 we first compute v_0(x) = v(x) mod x^{N/2} - 1 v_1(x) = v(x) mod x^{N/2} + 1 u_0(x) = u(x) mod x^{N/2} - 1 u_1(x) = u(x) mod x^{N/2} + 1

Now w_0 is the convolution of v_0 and u_0 which we can compute recursively. For w_1 we compute the negacyclic convolution v_1(x)u_1(x) mod x^{N/2} + 1 using Karatsuba.

There are 2 possible approaches to applying Karatsuba which mirror the DIT vs DIF approaches to FFT’s, the left/right decomposition or the even/odd decomposition. The latter seems to have fewer operations and so it is the one implemented below, though it does require a bit more data manipulation. It works as follows:

Define the even v_e and odd v_o parts so that v(x) = (v_e(x^2) + x v_o(x^2)). Then v(x)u(x) = (v_e(x^2)u_e(x^2) + x^2 v_o(x^2)u_o(x^2)) + x ((v_e(x^2) + v_o(x^2))(u_e(x^2) + u_o(x^2)) - (v_e(x^2)u_e(x^2) + v_o(x^2)u_o(x^2))) This reduces the problem to 3 negacyclic convolutions of size N/2 which can be computed recursively.

Of course, for small sizes we just explicitly write out the O(n^2) approach.

Traits

• Convolve
  Template function to perform convolution of vectors.
• RngElt
  This trait collects the operations needed by Convolve below.