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

ConvolutionElt

Bound alias for the wide operand type (used for both lhs and output).

ConvolutionRhs

Bound alias for the narrow operand type (rhs only).

Convolve

Trait for computing cyclic and negacyclic convolutions.

Functions

mds_circulant_karatsuba_16

Circulant matrix-vector multiply for width 16 via Karatsuba convolution.

mds_circulant_karatsuba_24

Circulant matrix-vector multiply for width 24 via Karatsuba convolution.