고속 푸리에 변환, 이론부터 구현까지

1. 구현하기 전에 필요한 지식

 

 고속 푸리에 변환(Fast Fourier Transform; FFT)은 어떤 벡터의 이산 푸리에 변환을 \( O(n \log n) \) 안에 구할 수 있게 해주는 쉬우면서도 강력한 알고리즘이다. 푸리에 변환이라고 하니 굉장히 생소해 보이고, 신호 처리나 미분방정식 푸는 데 외에는 쓸데가 없어 보이지만, 사실 이걸 갖고 상당히 많은 일을 할 수가 있다. 예를 들자면

 

• 차수가 n인 두 다항식을 \( O(n \log n) \)에 곱한다. (가장 직관적인 방법으로 계산하면 \( O(n^2) \)이 나오게 된다.)
• 크기가 \( 10^n \) 정도인 두 자연수를 \( O(n \log n) \)에 곱한다. (마찬가지로 직관적으로 곱하면 \( O(n^2) \)이다.)
• 자연수 집합의 두 부분집합 A와 B에 대하여 \( A + B := \{ x + y | x \in A, y \in B \} \)을 \( O(n \log n) \)에 구한다. n은 A와 B의 합집합 내부에 속한 원소 중 최댓값이다. (마찬가지로 직관적인 방법으로는 최악의 경우 \( O(|A||B|) \)이 나온다.
• [1, 0, 0, 1, 1, 0, 1, 1] 같은 투명 테이프를 생각해 보자. 이 투명 테이프에서 0은 그 부분이 색칠되어서 불투명함을, 1은 그 부분이 투명함을 나타낸다. 이제 서로 다른 투명 테이프 2개(길이 n, k; n > k)를 겹쳤을 때, 겹쳐졌으면서도 투명한 부분이 가장 많이 나오도록 하는 문제를 생각해보자. 이 문제를  \( O(n \log n) \)에 풀 수 있다. (마찬가지로 직관적인 방법으로는 \( O(nk) \)가 나온다.

 이 문제들을 어떻게 푸는지 대략적으로나마 이해하는 것은 중요하다. 고속 푸리에 변환은 두 벡터의 합성곱(convolution)을 계산하는 데 유용한 알고리즘이며, 이 문제들을 푸는 것은 이 문제와 두 벡터의 convolution을 빠르게 계산하는 문제가 동치라는 사실을 이해해야 하기 때문이다.

 일단 두 벡터 a와 b의 합성곱은 다음과 같이 정의된다. a와 b의 길이는 달라도 된다.

 

         \( c = a * b \) where \( c_k = \sum_{i+j=k} a_i b_j \)

 

 예컨대 a = [1, 0, 1], b = [1, 0, 1]일 경우, a * b를 계산하려면 다음 표를 그린다.

			1
		0		1
	1		0		0
*		1		0		1
	1		0		0
		0		1
			1
세로 방향의 합		1	0	2	0	1

 따라서 a * b = [1, 0, 2, 0, 1]이다. 눈치가 빠른 사람은 알겠지만, a * b의 계산 양상은 다항식을 곱할 때의 계산 양상과 굉장히 비슷하다.

 표가 생긴 것만 봐도 알겠지만, a * b를 계산하는 데는 O(len(a) * len(b))의 시간이 걸릴 것이라는 것을 알 수 있다. convolution을 더 빨리 계산하는 것이 목표이다. convolution의 계산 양상이 다항식 곱셈의 계산 양상과 상당히 비슷하다는 점을 관찰한 뒤, 다음의 발상을 해 보자. 

 

1. a를 그에 대응하는 다항식 f로 바꾸고, b도 다항식 g로 바꾼다. 예를 들어, [1, 0, 1] -> \( 1 + 0x + 1x^2 \).
2. a * b에 대응하는 다항식도 존재할 것이다. 그것을 h라 하자. h가 뭔지는 알 수 없지만, 우리는 최소한 f(0) * g(0) = h(0)이라는 사실을 알 수 있다. f(1) * g(1) = h(1)이라는 사실도 알 수 있다. 이것이 모든 정수, 실수, 심지어 복소수에 대해서도 성립한다.
3. 따라서 우리는 다음 두 문제를 빠르게 풀어야 한다: f가 주어지면 [x0, x1, ..., xk]를 [f(x0), ..., f(xk)]로 빠르게 바꾸는 문제, 그리고 그 반대.

 한 가지 다행인 점은 우리가 x0, x1, ...들을 맘대로 설정할 수 있다는 것이다. 이 값들을 아무렇게나 잡는다면, 대부분의 경우 목적을 달성하는 데 실패한다. 잘 알려진 방법은, x0, x1, ...들을 1의 거듭제곱근으로 잡는 것이다. 1의 거듭제곱근을 넣으면 된다는 사실을 어떻게 발견함? 이라고 물으면, 연속함수에 대한 푸리에 변환이라는 복잡한 주제에서 튀어나왔다는 사실이 잘 알려져 있다. 그러나 이것을 몰라도 FFT를 구현하는 것은 충분히 할 수 있으므로 넘어가도록 하자. 각설하고 1의 n제곱근은 다음과 같이 정의된다.

 

$$ \omega_k = \cos\left({2\pi k \over n}\right) + i \sin \left(2\pi k \over n\right), k = 0 ... n-1 $$

 

 이것은 제곱근의 정의에 따라 자명하다. 주어진 복소수는 항상 크기가 1이고, n제곱하면 편각이 0이 되기 때문이다.

1의 8제곱근의 모습, 당연하지만 8개가 존재한다.

 그런데 고속 푸리에 변환을 할 때 x가 1의 거듭제곱근이라고 해서 항상 \( f(\omega_0), f(\omega_1), ...\)이 빠르게 계산되는 것은 아니고, 항상 n이 2의 거듭제곱이어야 한다. 이유는 이 알고리즘이 분할 정복 알고리즘이기 때문에, 2의 거듭제곱이어야 절반으로 깔끔하게 나누어떨어지기 때문이다. 하지만 n이 2의 거듭제곱이 아닌 값으로 주어졌다고 좌절할 필요는 없으며, 값이 0인 항들을 몇 개 더 추가해 줘서 강제로 사이즈를 2의 거듭제곱으로 맞춰주면 된다.

 

 \( f(\omega_0), f(\omega_1), ...\)들을 빠르게 계산하기 위해서는 이 값들을 하나하나씩 계산하면 안 된다는 것이 자명하고, 한꺼번에 계산해야 할 것이다. 따라서 matrix form을 사용할 것이다. 예시를 바로 보여주기 위해 n = 8일 때 수식을 바로 적도록 하겠다.

 

$$ \begin{bmatrix} \omega_0 & \omega_0 & \omega_0 & \omega_0 & \omega_0 & \omega_0 & \omega_0 & \omega_0 \\ \omega_0 & \omega_1 & \omega_2 & \omega_3 & \omega_4 & \omega_5 & \omega_6 & \omega_7 \\ \omega_0 & \omega_2 & \omega_4 & \omega_6 & \omega_0 & \omega_2 & \omega_4 & \omega_6 \\ \omega_0 & \omega_3 & \omega_6 & \omega_1 & \omega_4 & \omega_7 & \omega_2 & \omega_5 \\ \omega_0 & \omega_4 & \omega_0 & \omega_4 & \omega_0 & \omega_4 & \omega_0 & \omega_4 \\ \omega_0 & \omega_5 & \omega_2 & \omega_7 & \omega_4 & \omega_1 & \omega_6 & \omega_3 \\ \omega_0 & \omega_6 & \omega_4 & \omega_2 & \omega_0 & \omega_6 & \omega_4 & \omega_2 \\ \omega_0 & \omega_7 & \omega_6 & \omega_5 & \omega_4 & \omega_3 & \omega_2 & \omega_1 \end{bmatrix} \begin{bmatrix} a_0 \\ a_1 \\ a_2 \\ a_3 \\ a_4 \\ a_5 \\ a_6 \\ a_7 \end{bmatrix} $$

 

 왼쪽 행렬을 일단 \(W_8\)이라고 이름붙이고, 오른쪽 벡터는 \(A_8\)이라고 하자.

 행렬을 잘 관찰해 보면, 어떤 행은 첨자가 모두 짝수인 반면 어떤 행은 홀수 첨자도 있고 짝수 첨자도 있는 것을 관찰할 수 있다. 이것을 분류 기준으로 삼아서 잘 모아줄 것이다. 행렬의 column의 순서를 바꿔 주는 column operation을 진행하자.

 

$$ \begin{bmatrix} \omega_0 & \omega_0 & \omega_0 & \omega_0 & \omega_0 & \omega_0 & \omega_0 & \omega_0 \\ \omega_0 & \omega_2 & \omega_4 & \omega_6 & \omega_1 & \omega_3 & \omega_5 & \omega_7 \\ \omega_0 & \omega_4 & \omega_0 & \omega_4 & \omega_2 & \omega_6 & \omega_2 & \omega_6 \\ \omega_0 & \omega_6 & \omega_4 & \omega_2 & \omega_3 & \omega_1 & \omega_7 & \omega_5 \\ \omega_0 & \omega_0 & \omega_0 & \omega_0 & \omega_4 & \omega_4 & \omega_4 & \omega_4 \\ \omega_0 & \omega_2 & \omega_4 & \omega_6 & \omega_5 & \omega_7 & \omega_1 & \omega_3 \\ \omega_0 & \omega_4 & \omega_0 & \omega_4 & \omega_6 & \omega_2 & \omega_6 & \omega_2 \\ \omega_0 & \omega_6 & \omega_4 & \omega_2 & \omega_7 & \omega_5 & \omega_3 & \omega_1 \end{bmatrix} \begin{bmatrix} a_0 \\ a_2 \\ a_4 \\ a_6 \\ a_1 \\ a_3 \\ a_5 \\ a_7 \end{bmatrix} $$

 

 이제 행렬을 강제로 가로 세로로 반반 나눌 것이다. 그러면 주어진 행렬식은 Block Matrix들로 이루어진 2 x 2 matrix의 곱으로 해석되어, 다음과 같은 형태로 변한다.

 

$$ F(A) = \begin{bmatrix} ? & ? \\ ? & ? \end{bmatrix} \begin{bmatrix} B \\ C \end{bmatrix} $$

$$ B = \begin{bmatrix} a_0 \\ a_2 \\ a_4 \\ a_6 \end{bmatrix}, C = \begin{bmatrix} a_1 \\ a_3 \\ a_5 \\ a_7 \end{bmatrix} $$

 

 여기서 중요한 것은 왼쪽 절반에 해당하는 행렬이 \(W_4\)라는 사실이다! 따라서 분할 정복을 적용하기에 너무 좋은 환경이 완성되었다. 그렇다면 왼쪽 절반을 다음과 같이 채울 수 있다.

 

$$ F(A) = \begin{bmatrix} W_4 & ? \\ W_4 & ? \end{bmatrix} \begin{bmatrix} B \\ C \end{bmatrix} $$ 

$$ B = \begin{bmatrix} a_0 \\ a_2 \\ a_4 \\ a_6 \end{bmatrix}, C = \begin{bmatrix} a_1 \\ a_3 \\ a_5 \\ a_7 \end{bmatrix} $$

 

 그렇다면 오른쪽 절반은? 오른쪽 절반을 채우는 일은 약간 까다롭지만, 선형대수학적 직관이 있다면 오른쪽 절반에 해당하는 행렬에 약간의 raw operation을 해 줘서 \(W_4\)와 어떤 특수한 행렬의 곱으로 나타낼 수 있다. 직접 계산을 통해 다음 식이 성립함을 확인해 보자.

 

$$ F(A) = \begin{bmatrix} W_4 & DW_4 \\ W_4 & -DW_4 \end{bmatrix} \begin{bmatrix} B \\ C \end{bmatrix} $$

$$ D = \begin{bmatrix} \omega_0 & 0 & 0 & 0 \\ 0 & \omega_1 & 0 & 0 \\ 0 & 0 & \omega_2 & 0 \\ 0 & 0 & 0 & \omega_3  \end{bmatrix} $$

$$ B = \begin{bmatrix} a_0 \\ a_2 \\ a_4 \\ a_6 \end{bmatrix}, C = \begin{bmatrix} a_1 \\ a_3 \\ a_5 \\ a_7 \end{bmatrix} $$

 

 당연하지만, D는 각 행에 상수 배를 취해 주는 elementary row operation들의 결합이다. 위의 식을 정리하면, 최종적으로는 다음과 같은 우아한 결론을 얻는다.

 

$$ F(A) = \begin{bmatrix} W_4 B + DW_4 C \\ W_4 B - DW_4 C \end{bmatrix} $$

 

 그런데 \( W_4(X) \)도 결국에는 \( F(X) \)로 치환할 수 있다는 것을 생각해보면, 결국 다음을 얻는 것이다.

 

$$ F(A) = \begin{bmatrix} F(B) + DF(C) \\ F(B) - DF(C) \end{bmatrix} $$

 

 이 공식은, FFT를 구현하는 데 있어서 필요한 사실상 전부라고 봐도 무방하다.

 

2. 구현하기

 PS를 즐겨 하는 사람들의 경험에 의하면 FFT는 상당히 무거운 알고리즘이라고 한다. 이유는 잘 모르겠다. 머지 소트든 FFT든 시간복잡도는 \( O(n \log n) \)인 건 똑같을 텐데... 하지만 같은 시간복잡도여도 FFT는 내부에서 보다 더 복잡한 연산을 거치기 때문에 더 느린 듯하다. 내가 n = 30만일 때 FFT를 돌려 보니 백준 환경에서 348ms가 나왔다. std::sort는 n = 100만일 때 296ms가 나오니, n = 30만일 때 89ms 정도 나올 것이라고 추측할 수 있는데... 측정이 정확한 건 아니지만 하여튼 FFT는 무겁다는 사실을 알아두자.

 

 여하튼 FFT는 무겁기 때문에 잘 구현하는 것이 중요하다. (만약 PS 목적으로 FFT를 구현하면 잘못하면 시간 초과가 날 수도 있다.) 따라서 FFT를 구현할 때는 보통 재귀를 사용하지 않고 구현하며, 가능한 한 새로운 메모리를 할당하지 않고 in-place로 구현하려고 한다. 이 목적을 모두 달성하면서 FFT를 구현하려고 하다 보니 세간에 떠도는 FFT 구현들은 비슷한 형태로 수렴진화하는 경향성이 있다. 그래서 이 블로그에서 소개하는 FFT도 일단 내가 짠 거긴 한데, 다른 데서 소개하는 FFT와 코드가 비슷할 것이다.

 

 이 코드는 안티 라크소넨의 [알고리즘 트레이닝] (2판) (2022)을 기반으로 내가 공부하여 새롭게 짠 것이다. (단 몇몇 최적화가 빠져 있어서 지금 내가 쓰는 버전이랑은 다르긴 하다...)

 

using cd = std::complex<double>;

std::vector<cd> fft(const std::vector<cd> &a, int pol = 1) {
    int n = 1;
    const int a_sz = static_cast<int>(a.size());
    while (n < a_sz) {
        n <<= 1;
    }
    std::vector<cd> r(n);
    for (int i = 0; i < n; ++i) {
        int s = 0;
        for (int b = 1, d = n / 2; b < n; b <<= 1, d >>= 1) {
            if (b & i) {
                s += d;
            }
        }
        r[s] = i < a_sz ? a[i] : 0;
    }
    
    for (int mte = 2; mte <= n; mte <<= 1) {
        cd we = std::exp(cd(0, pol * 2 * pi / mte));
        for (int i = 0; i < n; i += mte) {
            cd w = 1;
            for (int j = i; j < i + mte / 2; ++j, w *= we) {
                int j2 = j + mte / 2;
                cd temp1 = r[j], temp2 = r[j2];
                r[j] = temp1 + w * temp2;
                r[j2] = temp1 - w * temp2;
            }
        }
    }
    
    if (pol == -1) {
        for (auto &c : r) {
            c /= n;
        }
    }
    
    return r;
}

 

 FFT의 비재귀 구현은 다음의 2페이즈로 나뉜다.

 

1. Divide and Conquer 알고리즘은 보통 Divide -> Conquer -> Merge의 세 단계를 거치는데, FFT의 경우 Merge를 어떤 순서로 해 주느냐가 약간 까다롭다. (Merge sort처럼 0, 1 merge하고 2, 3 merge하고 이런 게 아니다) 그 순서에 맞춰 Merge할 수 있게끔 미리 원소를 옮겨준다.
2. 위에서 썼던 공식을 때려박아준다.

 여기서 2에 나오는 '위에서 썼던 공식'을 다시 한 번 소개한다.

 

$$ F(A) = \begin{bmatrix} F(B) + DF(C) \\ F(B) - DF(C) \end{bmatrix} $$

 

 다음 코드는 이 공식을 사실상 그대로 코드로 옮긴 거라고 봐도 무방하다.

 

    for (int mte = 2; mte <= n; mte <<= 1) {
        cd we = std::exp(cd(0, pol * 2 * pi / mte));
        for (int i = 0; i < n; i += mte) {
            cd w = 1;
            for (int j = i; j < i + mte / 2; ++j, w *= we) {
                int j2 = j + mte / 2;
                cd temp1 = r[j], temp2 = r[j2];
                r[j] = temp1 + w * temp2;
                r[j2] = temp1 - w * temp2;
            }
        }
    }

 일단 mte라는 변수가 있는데, 한 번 merge할 때 몇 개의 수를 묶어서 뭉탱이로 처리할 것이냐를 나타내는 변수다. 당연히 2 -> 4 -> 8 -> 16 ... 순서대로 가야 한다. 그리고 we는 단위근이다. 실제로 we를 mte제곱해 보면 당연히 1이 나오는 것을 확인할 수 있을 것이다. 만약 mte = 8이라면 [0, 1, ..., 7], [8, 9, ..., 15], [16, 17, ..., 23], ... 개의 수를 묶어서 뭉탱이로 처리하게 될 텐데, i는 여기서 0, 8, 16...에 해당한다. 즉 i는 각 뭉탱이의 첫 번째 수다.

 이 코드가 위에 적힌 행렬 연산을 그대로 코드로 옮긴 건 맞지만, 구현할 때 행렬 연산을 그대로 구현한 다음에 그걸 기반으로 FFT를 돌리는 것은 별로 권장되는 방법이 아니다. 예를 들어 다음 알고리즘이 있다고 하자.

 

 [0, 1, 2, 3, 4, 5, 6, 7]을 [0, 1, 2, 3], [4, 5, 6, 7]로 나눔 -> [4, 5, 6, 7]에 D를 곱함 -> [0, 1, 2, 3]+D[4, 5, 6, 7]과 [0, 1, 2, 3]-D[4, 5, 6, 7]에 해당하는 벡터를 각각 만듦 -> 원래 자리에 대입

 

 이렇게 숫자들을 한꺼번에 움직이는 방식은 별로 좋지 않다. 숫자들을 하나하나씩 움직일 수도 있기 때문이다.

 

 [0, 4]를 [0+4w0, 0-4w0]으로 처리 -> [1, 5]를 [1+5w1, 1-5w1]로 처리 -> [2, 6]을 [2+6w2, 2 - 6w2]로 처리 -> [3, 7]을 [3 + 7w3, 3 - 7w3]으로 처리

 

 이 편이 훨씬 공간복잡도가 적게 든다. 추가되는 메모리가 O(1)이다.

 

 여기까지 이해했으면 2페이즈 중 2번째 페이즈를 이해한 것이다. 이제 첫 번째 페이즈를 이해하러 가 보자.

 벡터 [1, 2, 3, 4, 5, 6, 7, 8]은 그 홀짝성에 의해, 묶일 때 짝수 인덱스의 수는 모두 왼쪽으로, 홀수 인덱스의 수는 모조리 오른쪽으로 이동하게 된다. 즉 다음과 같이 바뀐다.

 

 [0, 2, 4, 6, 1, 3, 5, 7]

 

 이 중 [0, 2, 4, 6]은 다시 짝수 인덱스의 수를 모두 왼쪽으로, 홀수 인덱스의 수를 모두 오른쪽으로 옮겨서 다음과 같이 바뀐다.

 

 [0, 4, 2, 6, 1, 5, 3, 7]

 

 (옮기는 김에 1, 3, 5, 7도 같이 처리해 줬다.) 이렇게만 써놓으면 패턴이 뭔지 눈에 잘 안 들어올 수 있는데, 2진법으로 써 보면 다음과 같다.

 

 [000, 100, 010, 110, 001, 101, 011, 111]

 

 각 숫자들을 오른쪽에서 왼쪽으로 읽으면 다음과 같다.

 

 [000, 001, 010, 011, 100, 101, 110, 111]

 

 이것은 놀랍게도 그냥 자연수를 나열한 것과 일치한다. 따라서 우리는 그냥 자연수를 나열한 수열을 하나 만든 뒤, 그 수열의 원소를 하나하나 골라 일일이 비트를 좌우대칭시켜줄 필요가 있다. 그것이 1페이즈가 하는 일이다.

 

    int n = 1;
    const int a_sz = static_cast<int>(a.size());
    while (n < a_sz) {
        n <<= 1;
    }
    std::vector<cd> r(n);
    for (int i = 0; i < n; ++i) {
        int s = 0;
        for (int b = 1, d = n / 2; b < n; b <<= 1, d >>= 1) {
            if (b & i) {
                s += d;
            }
        }
        r[s] = i < a_sz ? a[i] : 0;
    }

 

 1페이즈의 시간복잡도는 2페이즈 코드의 시간복잡도와 마찬가지로 \( O(n \log n) \)이다.

 FFT가 상당히 중요한 알고리즘이기 때문에 많은 학자들과 PS 헤비 유저들이 FFT를 최적화하기 위해 노력했다. 그 결과 이런 비트 좌우반전시켜주기 연산을 최적화하기만 해도, 그 자체가 하나의 논문감이 되는 것으로 보인다. 실제로 이를 최적화하기 위한 학자들의 눈물의 쇼를 구경하고 싶다면, 구글 스칼라에 bit reverse permutation이라고 검색해 보면 무슨 500회, 200회 인용된 논문이 튀어나오고 그런다. 문제 이해 자체는 중학생도 할 만큼 쉬운데 중요도가 이렇게 높다니 참 신기하다는 생각이 든다.

 

 3. 아까 전에 역변환도 해야 한다고 하지 않음?

 

 그렇다. 방금 전에 우리는 \( \omega_0, \omega_1, ... \)를 \( f(\omega_0), f(\omega_1), ... \)로 바꿔 주어야 할 뿐만 아니라 그 반대도 해야 한다고 했었다. 그리고 그것은 이 행렬의 역행렬을 구하면 해결된다.

 

$$ \begin{bmatrix} \omega_0 & \omega_0 & \omega_0 & \omega_0 & \omega_0 & \omega_0 & \omega_0 & \omega_0 \\ \omega_0 & \omega_1 & \omega_2 & \omega_3 & \omega_4 & \omega_5 & \omega_6 & \omega_7 \\ \omega_0 & \omega_2 & \omega_4 & \omega_6 & \omega_0 & \omega_2 & \omega_4 & \omega_6 \\ \omega_0 & \omega_3 & \omega_6 & \omega_1 & \omega_4 & \omega_7 & \omega_2 & \omega_5 \\ \omega_0 & \omega_4 & \omega_0 & \omega_4 & \omega_0 & \omega_4 & \omega_0 & \omega_4 \\ \omega_0 & \omega_5 & \omega_2 & \omega_7 & \omega_4 & \omega_1 & \omega_6 & \omega_3 \\ \omega_0 & \omega_6 & \omega_4 & \omega_2 & \omega_0 & \omega_6 & \omega_4 & \omega_2 \\ \omega_0 & \omega_7 & \omega_6 & \omega_5 & \omega_4 & \omega_3 & \omega_2 & \omega_1 \end{bmatrix} $$

 

 이 글에서는 그냥 남들이 구해놓은 역행렬을 베껴오도록 하자. (이 과정이 그렇게 중요한 것 같지 않다) 결론만 말하면, 위 행렬의 모든 성분을 역수 취해준 다음 1/n을 곱하면 그게 역행렬이 된다. 따라서 위의 FFT 구현에서 \( \omega_1 \)을 \( 1 \over \omega_1 \)로 바꿔준 다음 결과를 n으로 나눠주면 된다. 위의 FFT 구현은 pol이라는 인자를 받고 그 디폴트 값이 1로 지정되어 있는데, 만약 pol 인자의 값을 -1로 설정해 준다면 이 모든 절차가 알아서(...) 진행이 된다. 이는 원래 행렬과 역행렬이 상당히 비슷하게 생겼기 때문이다.

 

 4. 말 나온 김에 convolution까지 구해보자!

 위에서 convolution 구하는 방법도 설명했는데, 이를 그대로 구현하면 다음과 같다.

std::vector<cd> mul(std::vector<cd> a, std::vector<cd> b) {
    const int target_size = static_cast<int>(a.size() + b.size() - 1);
    a.resize(target_size, 0);
    b.resize(target_size, 0);
    
    auto ta = fft(a), tb = fft(b);
    
    std::vector<cd> ret(ta.size());
    for (int i = 0; i < ret.size(); ++i) {
        ret[i] = ta[i] * tb[i];
    }
    
    return fft(ret, -1);
}