Introducing the imaginary cyclic convolution, speeding up convolution by a factor of 2
Разница между en7 и en8, 7 символ(ов) изменены
I recently had a very interesting idea for how to greatly speed up convolution (a.k.a. polynomial multiplication). ↵

```py↵
def convolution(A, B):↵
  C = [0] * (len(A) + len(B) - 1)↵
  for i in range(len(A)):↵
    for j in range(len(B)):↵
      C[i + j] += A[i] * B[j]↵
  return C↵
```↵

The standard technique for how to do convolution fast is to make use of cyclic convolution (polynomial mult mod $x^n - 1$). ↵

```py↵
def cyclic_convolution(A, B):↵
  n = len(A) # A and B needs to have the same size↵
  C = [0] * n↵
  for i in range(n):↵
    for j in range(n):↵
      C[(i + j) % n] += A[i] * B[j]↵
  return C↵
```↵

Cyclic convolution can be calculated in $O(n \log n)$ using FFT, which is really fast. The issue here is that in order to do convolution using cyclic convolution, we need to pad with a lot of 0s to not be affected by the wrap around. All this 0-padding feels very inefficient.↵

So here is my idea. What if we do polynomial multiplication mod $x^n - i$ instead of mod $x^n - 1$? Then when we get wrap around, it will be multiplied by the imaginary unit, so it wont interfere with the real part! I call this the _imaginary cyclic convolution_.↵

```py↵
def imaginary_cyclic_convolution(A, B):↵
  n = len(A) # A and B needs to have the same size↵
  C = [0] * n↵
  for i in range(n):↵
    for j in range(n):↵
      C[(i + j) % n] += (1 if i + j < n else 1j) * A[i] * B[j]↵
  return C↵
```↵

Imaginary cyclic convolution is the perfect algorithm to use for implementing convolution. Using it, we no longer need to do copious amount of 0 padding, since the imaginary unit will take care of the wrap around. In fact, the size (the value of $n$) required is exactly half of what we would need if we had used cyclic convolution.↵

One question 
still remains, how do we implement imaginary cyclic convolution efficiently? ↵

The trick is rather simple. Let $\omega = i^{\frac{1}{n}}$. Now note that if $C(\omega x) = A(\omega x) B(\omega x) \mod x^n - 1$ then $C(x) = A(x) B(x) \mod x^n - i$. So here is the algorithm ↵

```py↵
def imaginary_cyclic_convolution(A, B):↵
  n = len(A) # A and B needs to have the same size↵
  w = (1j)**(1/n) # n-th root of imaginary unit↵
  ↵
  # Transform the polynomial
s A(x) -> A(wx) and B(x) -> B(wx)↵
  A = [A[i] * w**i for i in range(n)]↵
  B = [B[i] * w**i for i in range(n)]↵

  C = cyclic_convolution(A, B)↵
  ↵
  # Transform the polynomial C(wx) -> C(x)↵
  C = [C[i] / w**i for i in range(n)]↵
  return C↵
```

История

 
 
 
 
Правки
 
 
  Rev. Язык Кто Когда Δ Комментарий
en9 Английский pajenegod 2022-09-16 16:34:17 0 (published)
en8 Английский pajenegod 2022-09-16 16:33:59 7
en7 Английский pajenegod 2022-09-16 15:12:46 2 Tiny change: '< n else 1i) * A[i] *' -> '< n else 1j) * A[i] *'
en6 Английский pajenegod 2022-09-16 14:57:46 126
en5 Английский pajenegod 2022-09-16 14:36:07 32
en4 Английский pajenegod 2022-09-16 14:32:50 8 Tiny change: 'n\n C = circular_convoluti' -> 'n\n C = cyclic_convoluti'
en3 Английский pajenegod 2022-09-16 14:29:11 105
en2 Английский pajenegod 2022-09-16 14:23:27 6 Tiny change: 'len(A)\n omega = (1i)**(' -> 'len(A)\n w = (1i)**('
en1 Английский pajenegod 2022-09-16 14:17:53 2316 Initial revision (saved to drafts)