고속 푸리에 변환 (Fast Fourier Transform, FFT)

들어가기 전에

서로 수직인 (영벡터가 아닌) $n$차원 벡터 $\mathbf v_1$, $\mathbf v_2$, …, $\mathbf v_n$을 생각해봅시다. 임의의 $n$차원 벡터 $\mathbf a$는 이들의 일차 결합으로 항상 나타낼 수 있습니다.

\[ \mathbf a = c_1 \mathbf v_1 + c_2 \mathbf v_2 + \cdots + c_n \mathbf v_n \]

이때 계수 $c_k$는 양변에 $\mathbf v_k$를 내적하여 구할 수 있습니다.

\[ \mathbf a \boldsymbol\cdot \mathbf v_k = c_1 \mathbf v_1 \boldsymbol\cdot \mathbf v_k +
c_2 \mathbf v_2 \boldsymbol\cdot \mathbf v_k + \cdots + c_n \mathbf v_n \boldsymbol\cdot \mathbf v_k
= c_k \mathbf v_k \boldsymbol\cdot \mathbf v_k \]

\[ \therefore c_k = \frac{ \mathbf a \boldsymbol\cdot \mathbf v_k }{ \mathbf v_k \boldsymbol\cdot \mathbf v_k } \]

역으로, $c_1$, $c_2$, …, $c_n$의 값을 바꿔주면 임의의 $n$차원 벡터를 얻을 수 있습니다. 이 때 $\mathbf v_1$, $\mathbf v_2$, …, $\mathbf v_n$을 $n$차원 공간의 기저라고 합니다.

이산 푸리에 변환

수열의 내적과 직교

일반적으로 내적은 벡터에서만 정의하는 것이 아닙니다. 아주 다양한 대상에 내적을 정의하여 해당 대상을 벡터처럼 다룰 수 있습니다. 주기가 $N$인 복소 수열 (앞으로 별 말이 없으면 수열은 항상 주기가 $N$인 복소 수열을 말합니다.) $\mathbf a$와 $\mathbf b$의 내적을 아래와 같이 정의합시다.

\[ \mathbf a \boldsymbol\cdot \mathbf b = a(0) \bar b(0) + a(1) \bar b(1) + \cdots a(N-1) \bar b(N-1) = \sum_{n=0}^{N-1} a(n) \bar b(n) \]

이때 $\bar b(n)$은 $b(n)$의 켤레 복소수입니다. 벡터의 내적과 아주 유사하죠. 마찬가지로 두 수열의 내적이 0이면 직교한다고 합니다.

이제 정수 $k$에 대해 수열 $\boldsymbol \phi_k$를

\[ \phi_k(n) = \exp \left( \frac{2\pi k}{N} ni \right)
= \cos \left( \frac{2\pi k}{N} n \right) + i \sin \left ( \frac{2\pi k}{N} n \right) \]

이라 하면 이 수열은 주기가 $N$이고 (삼각함수의 성질) $\boldsymbol\phi_0$, $\boldsymbol\phi_1$, …, $\boldsymbol\phi_{N-1}$은 서로 직교합니다.

증명) 서로 다른 정수 $k$와 $l$에 대해

\[ \begin{align}
\boldsymbol \phi_k \boldsymbol\cdot \boldsymbol\phi_l
&= \sum_{n=0}^{N-1} \phi_k(n)\bar\phi_l(n) \\
&= \sum_{n=0}^{N-1} \exp \left( \frac{2\pi k}{N}ni \right) \exp \left( – \frac{2\pi l}{N} ni \right) \\
&= \sum_{n=0}^{N-1} \exp \left[ \frac{2\pi (k-l)}{N} ni \right] \\
&= \frac{\exp\left[ 2\pi (k-l)i \right]-1}{\exp\left[\frac{2\pi}{N}(k-l)i\right]-1} \\
&= 0
\end{align} \]

이므로 $\boldsymbol \phi_k$와 $\boldsymbol\phi_l$은 직교합니다. ■

한편

\[ \boldsymbol\phi_k \boldsymbol\cdot \boldsymbol\phi_k = \sum_{n=0}^{N-1} \exp(0) = N \]

입니다.

이산 푸리에 변환과 역변환

수학자 조제프 푸리에는 $\boldsymbol\phi_0$, $\boldsymbol\phi_1$, …, $\boldsymbol\phi_{N-1}$을 기저로 사용하여 임의의 수열을 만들어낼 수 있다고 했습니다. (정확히는, 푸리에가 처음 말한 건 수열이 아닌 함수에 대한 것입니다.) 즉, 수열 $\mathbf a$에 대해 수열 $\mathbf A$를

\[ A(n) = a(0)\phi_0(n) + a(1)\phi_1(n) + \cdots + a(N-1)\phi_{N-1}(n) = \sum_{k=0}^{N-1} a(k) \phi_k(n) \tag 1 \]

으로 정의하면 $\mathbf a$를 잘 바꿔서 $\mathbf A$가 어떤 수열이든 되게 할 수 있다는 뜻이죠. 역으로 $\mathbf A$를 주면 위 식을 만족하는 $\mathbf a$를 항상 찾을 수 있다는 소리도 되고요. 여기서 수열 $\mathbf A$를 수열 $\mathbf a$의 이산 푸리에 변환(Discrete Fourier Transform, DFT)이라 하고, $\mathbf A = \mathcal{F} \left\{ \mathbf a \right\}$로 표기합니다. 이때 $\mathbf a$는 $\mathbf A$의 역변환이 되고, $\mathbf a = \mathcal{F}^{-1} \left\{ \mathbf A \right\}$로 표기합니다. 역변환 식은 벡터와 비슷하게 식 (1)의 양변에 $\overline{\phi}_l(n)$을 곱하여 구할 수 있습니다.

\[ A(n) \bar\phi_l(n) = \sum_{k=0}^{N-1} a(k) \phi_k(n) \bar\phi_l(n) \]

\[ \sum_{n=0}^{N-1} A(n) \bar\phi_l(n) = \sum_{n=0}^{N-1} \sum_{k=0}^{N-1} a(k) \phi_k(n) \bar\phi_l(n) = Na(l) \]

\[ \therefore a(n) = \frac{1}{N} \sum_{k=0}^{N-1} A(k) \bar\phi_n(k)= \frac{1}{N} \sum_{k=0}^{N-1} A(k) \bar\phi_k(n) \tag 2\]

식 (2)를 잘 보면 식 (1)에서 $\phi_k(n)$ 대신에 $\bar\phi_k(n)$을 쓰고 앞에 $1/N$이 곱해진 것 말고는 아무 차이가 없습니다. 다시 말해 $\overline{\boldsymbol\phi}_0$, $\overline{\boldsymbol\phi}_1$, …, $\overline{\boldsymbol\phi}_{N-1}$을 기저로 써서 변환한다음 $1/N$을 곱하면 그게 바로 역변환이 됩니다. 따라서 역변환을 굳이 생각하지 않아도 됩니다!

합성곱

수열 $\mathbf a$와 $\mathbf b$의 이산 푸리에 변환의 곱 $\mathcal{F}\{\mathbf a\}\mathcal{F}\{\mathbf b\}$에 역변환을 취하여 나온 수열을 $\mathbf c$라 합시다. 그러면

\[ \begin{align}
\mathcal{F} \{ \mathbf c \} (n) &= \mathcal{F} \{ \mathbf a \} (n) \mathcal{F} \{ \mathbf b \} (n) \\
&= \left[ \sum_{k=0}^{N-1} a(k) \phi_k(n) \right] \left[ \sum_{l=0}^{N-1} b(l) \phi_l(n) \right] \\
&= \sum_{k=0}^{N-1} \sum_{l=0}^{N-1} a(k) b(l) \exp \left[ \frac{2\pi (k+l)}{N} n i \right] \\
&= \sum_{k=0}^{N-1} \underbrace{ \sum_{l=0}^{N-1} a(l) b(k-l)}_{c(k)} \ \underbrace{\exp \left( \frac{2\pi k}{N} ni \right)}_{\phi_k(n)}
\end{align} \]

\[ \therefore c(n) = \sum_{k=0}^{N-1} a(k)b(n-k) \tag 3 \]

(세 번째 줄에서 네 번째 줄로 넘어갈 때에는 삼각함수의 주기성을 사용합니다. 조금 유도하기 어렵습니다.) 여기서 $\mathbf c$를 $\mathbf a$와 $\mathbf b$의 합성곱(Convolution)이라 하며 $\mathbf c = \mathbf a * \mathbf b$로 표기합니다. 정의에 의해 $\mathcal{F} \{\mathbf a * \mathbf b \} = \mathcal{F} \{ \mathbf a \} \mathcal{F} \{ \mathbf b \}$입니다.

이 글에서 다루는 수열은 모두 주기가 $N$인 주기 수열이라는 점을 생각하면 합성곱을 이렇게 이해할 수도 있습니다. 위 그림처럼 $\mathbf b$가 오른쪽에서 왼쪽으로 늘어져 있고 그 위에 칸을 맞춰서 $a$가 딱 한 주기만큼만 왼쪽에서 오른쪽으로 늘어져 있습니다. 처음에는 $a(0)$와 $b(0)$가 겹치게 칸을 맞춰줍니다. 마주보는 칸끼리 곱한 다음 곱의 총합을 구하면 그게 $(\mathbf a * \mathbf b)(0)$가 됩니다. 그 다음 $\mathbf a$를 왼쪽으로 한 칸 옮겨 $a(0)$와 $b(1)$이 겹치게 하고 곱의 총합을 구하면 $(\mathbf a * \mathbf b)(1)$입니다. 이렇게 $n = 0, 1, 2, \cdots, N-1$에 대해 $a(0)$와 $b(n)$가 겹치게 놓고서 곱의 총합을 구하면 합성곱의 한 주기가 구해집니다.

푸리에 변환에서 합성곱이 갑자기 나오는 이유는 사실상 이산 푸리에 변환의 목적이 합성곱을 간단하게 구하기 위해서이기 때문입니다. 합성곱은 의외로 여기저기 많이 나오는 수열 연산 중 하나이기 때문에 이를 빠르게 구하는 건 중요합니다. 그래서 두 수열의 합성곱을 구하고 싶을 때 식 (3)을 직접 계산하는 것이 아니라 각각을 이산 푸리에 변환해주고, 변환 결과를 곱한 다음, 역변환을 취하자는 말이죠. 근데, 식 (1)이나 (2)를 보면 식 (3)으로 직접 계산하는 것보다 별로 간편해보이지가 않습니다. 정말 그럴까요?

고속 푸리에 변환

본격적으로 이산 푸리에 변환을 구해봅시다. 식 (1)에서 $A(n)$ 하나를 구하는 데 $O(N)$이 걸리므로 $\mathbf A$ 전체를 구하는 데에는 $O(N^2)$이면 충분하다는 결론이 나옵니다.

… 물론 이렇게 구현하면 끔찍하게 느립니다! 그래서 몇몇 학자들이 고민하던 와중 분할 정복을 사용해 시간 복잡도를 $O(N \log N)$으로 줄이자는 아이디어가 나왔고, 이 기법을 고속 푸리에 변환(Fast Fourier Transform)이라 합니다. 분할 정복만 사용하면 FFT란 이름이 붙기 때문에 아주 다양한 FFT 알고리즘이 있지만 가장 많이 쓰이는 것은 제임스 쿨리와 존 튜키가 1965년 발표한 쿨리-튜키 알고리즘입니다.

식 (1)을 다시 쓰면

\[ A(n) = \sum_{k=0}^{N-1} a(k) \phi_k(n) = \sum_{k=0}^{N-1} a(k) \exp \left( \frac{2\pi n}{N} ki \right)
= \sum_{k=0}^{N-1} a(k) \omega_N^{kn} \]

여기서 $\omega_N = \exp ( 2\pi i / N )$은 1의 $N$제곱근입니다. 분할 정복을 위해 $A(n)$을 두 부분으로 나눠봅시다. 편의상 $N$은 짝수라고 가정합니다.

\[ \begin{align}
A(n) &= a(0)\omega_N^0 + a(1)\omega_N^n + a(2)\omega_N^{2n} + \cdots + a(N-1)\omega_N^{(N-1)n} \\
&= \left[ a(0)\omega_N^0 + a(2)\omega_N^{2n} + a(4)\omega_N^{4n} + \cdots + a(N-2)\omega_N^{(N-2)n} \right] \\
& \qquad + \omega_N^n \left[ a(1)\omega_N^0 + a(3)\omega_N^{2n} + a(5)\omega_N^{4n} + \cdots
+ a(N-1)\omega_N^{(N-2)n} \right]
\end{align} \]

주기가 $N/2$인 복소 수열 $\mathbf a_\textrm{even}$과 $\mathbf a_\textrm{odd}$를

\[ \begin{align}
a_\textrm{even} (n) &= a(2n) \\
a_\textrm{odd} (n) &= a(2n+1)
\end{align}
\]

로 정의하면

\[ \begin{align}
A(n) &= \sum_{k=0}^{N/2-1} a_\textrm{even}(k) \omega_N^{2kn}
+ \omega_N^n \sum_{k=0}^{N/2-1} a_\textrm{odd}(k) \omega_N^{2kn} \\
&= \sum_{k=0}^{N/2-1} a_\textrm{even}(k) \omega_{N/2}^{kn}
+ \omega_N^n \sum_{k=0}^{N/2-1} a_\textrm{odd}(k) \omega_{N/2}^{kn} \\
&= \underbrace{\mathcal{F} \{ \mathbf a_\textrm{even} \} (n)}_{A_\textrm{even}(n)}
+ \omega_N^n \underbrace{\mathcal{F} \{ \mathbf a_\textrm{odd} \} (n)}_{A_\textrm{odd}(n)} \\
\end{align} \]

$\mathbf a_\textrm{even}$과 $\mathbf a_\textrm{odd}$를 이산 푸리에 변환한 결과인 $\mathbf A_\textrm{even}$와 $\mathbf A_\textrm{odd}$도 주기가 $N/2$이므로 위 식은 다음과 같이 적으면 더 편합니다.

\[ \begin{align}
A(n) &= A_\textrm{even}(n) + \omega_N^n A_\textrm{odd}(n) \\
A(n + N/2) &= A_\textrm{even}(n) \ – \ \omega_N^n A_\textrm{odd}(n)
\end{align} \tag 4 \]

따라서 재귀적으로 이산 푸리에 변환을 계산하려면

  1. 주어진 수열 $\mathbf a$를 $\mathbf a_\textrm{even}$과 $\mathbf a_\textrm{odd}$으로 나눈다.
  2. 나눠진 각각에 대해 이산 푸리에 변환을 계산한다.
  3. 식 (4)를 이용해 $\mathbf a$의 푸리에 변환을 계산한다.

시간 복잡도를 $T(N)$이라 하면 1번 과정에서 $O(N)$, 2번에서 $2T(N/2)$, 3번에서 $O(n)$만큼 걸리므로 $T(N) = 2T(N/2) + O(N)$이고, 마스터 정리를 써서 계산해보면 $T(N) = O(N \log N)$입니다!

이 방식을 파이썬으로 그대로 구현하면 아래와 같습니다. $N = 1$이면 변환 결과가 자기 자신이라는 데 유의합시다.

주의: 이 코드는 $N$이 2의 거듭제곱일 때만 동작합니다.

재귀 풀기

시간 복잡도는 줄었지만, 아직도 이 속도에 만족하지 못한 사람들이 있었습니다. 느린 속도의 원인은 이걸 재귀로 구현했다는 점이죠. 재귀는 간편하지만, 반복문에 비하면 계산 속도에서 손해를 보기 때문에 재귀를 풀어서 반복문으로 만들어줍시다.

$N=8$에 대해 FFT를 계산해주는 회로를 그리면

이때 식 (4)는 왼쪽과 오른쪽 회로가 동등하다는 걸 말해줍니다. 여기서 화살표에 적힌 수는 곱셈을 의미합니다. 마찬가지로 $N=4$에 대해 FFT를 계산해주는 회로는 $N=2$에 대해 FFT를 계산해주는 회로 두 개를 써서 그릴 수 있고, 이렇게 재귀적으로 그리면 다음 그림이 나옵니다. 이 계산법을 나비처럼 생겼다고 해서 버터플라이 알고리즘이라고 부릅니다.

재배열 두 단계를 거치면서 $\mathbf a$의 순서가 완전히 뒤죽박죽이 되었습니다. 여하튼 $\mathbf a$를 위와 같은 순서로 처음에 한 번만 재배열하고, 그다음부터는 이 회로를 그대로 반복문을 써서 시뮬레이션하면 됩니다. 관건은 $\mathbf a$의 순서가 대체 어떻게 되는가인데 순서의 이진수 표현에 주목합시다. 재배열 결과를 $\mathbf a_\textrm{rearr}$이라 하면

\[ \begin{matrix}
n & a_\textrm{rearr}(n) \\ \hline
000 & a(000) \\
001 & a(100) \\
010 & a(010) \\
011 & a(110) \\
100 & a(001) \\
101 & a(101) \\
110 & a(011) \\
111 & a(111)
\end{matrix} \]

따라서 $a_\textrm{rearr}(n) = a(n$의 마지막 $\log N$비트를 반전한 수$)$라고 유추할 수 있습니다.

증명) 먼저 $N=1$일 때 성립합니다. 그리고 $N=2^k$일 때 성립한다고 하다고 가정합시다. ($k \ge 0$)

$2N$에 대한 FFT를 하려면 먼저 $\mathbf a_\textrm{even}$과 $\mathbf a_\textrm{odd}$로 나눠야죠. 정의에 의해서 $n=0, 1, 2, \cdots, N-1$에 대해

\[ a_\textrm{even} ( \underbrace{\times\times\cdots\times}
_{\substack{\textrm{bit representation of } n \\ \textrm{total } k \textrm{ bits} }} )
= a [ 2 ( \underbrace{\times\times\cdots\times}_n ) ]
= a ( \underbrace{\times\times\cdots\times}_n 0 ) \]

\[ a_\textrm{odd} ( \underbrace{\times\times\cdots\times}
_{\substack{\textrm{bit representation of } n \\ \textrm{total } k \textrm{ bits} }} )
= a [ 2 ( \underbrace{\times\times\cdots\times}_n ) + 1 ]
= a ( \underbrace{\times\times\cdots\times}_n 1 ) \]

$\mathbf a_\textrm{even}$의 FFT를 계산하기 위해 재배열 한 것을 $\mathbf a_\textrm{even, rearr}$이라 하면 가정에 의해

\[ a_\textrm{even, rearr} ( \underbrace{\times\times\cdots\times}_{\textrm{bit repr. of } n} )
= a_\textrm{even} ( \underbrace{\times\times\cdots\times}_{\substack{\textrm{reversed} \\ \textrm{bit repr. of } n}} )
= a ( \underbrace{\times\times\cdots\times}_{n \textrm{ reversed}} 0 ) \]

같은 방법으로

\[ a_\textrm{odd, rearr} ( \underbrace{\times\times\cdots\times}_{\textrm{bit repr. of } n} )
= a_\textrm{odd} ( \underbrace{\times\times\cdots\times}_{\substack{\textrm{reversed} \\ \textrm{bit repr. of } n}} )
= a ( \underbrace{\times\times\cdots\times}_{n \textrm{ reversed}} 1 ) \]

$\mathbf a_\mathrm{rearr}$은 $\mathbf a_\textrm{even, rearr}$와 $\mathbf a_\textrm{odd, rearr}$을 이어붙인 것이므로

\[ a_\mathrm{rearr} (0 \underbrace{\times\times\cdots\times}_{\textrm{bit repr. of } n} )
= a_\textrm{even, rearr} ( \underbrace{\times\times\cdots\times}_{\textrm{bit repr. of } n} )
= a ( \underbrace{\times\times\cdots\times}_{n \textrm{ reversed}} 0 ) \]

\[ a_\mathrm{rearr} (1 \underbrace{\times\times\cdots\times}_{\textrm{bit repr. of } n} )
= a_\textrm{odd, rearr} ( \underbrace{\times\times\cdots\times}_{\textrm{bit repr. of } n} )
= a ( \underbrace{\times\times\cdots\times}_{n \textrm{ reversed}} 1 ) \]

가 되어 $2N=2^{k+1}$일 때에도 성립합니다. ■

정수 $n$의 마지막 $k$자리를 뒤집는 연산은 많이 쓰이기 때문에 여러 구현법이 있습니다. 여기서는 다음 구현을 쓰겠습니다. #출처

남은 건 정말로 구현뿐입니다. 그 전에 잠시 역변환을 어떻게 구현할 지 생각해봅시다. 역변환은 기저로  $\overline{\boldsymbol\phi}_0$, $\overline{\boldsymbol\phi}_1$, …, $\overline{\boldsymbol\phi}_{N-1}$을 쓴다고 했으므로 위 과정을 그대로 해주되 $\omega_N$ 대신 $\overline\omega_N = 1/\omega_N$을 넣어주면 됩니다. 그러면 FFT를 그대로 쓸 수 있죠. 아, 그리고 최종 결과를 $N$으로 나눠야 한다는 점 까먹지 마세요.

C++로 짠 FFT 코드입니다. is_reverse가 참이면 역변환을 계산해줍니다.

주의: 여전히 이 코드도 $N$이 2의 거듭제곱일 때에만 동작합니다.

N이 2의 거듭제곱이 아닌 경우

분할 정복의 특성상 지금까지는 $N$이 2의 거듭제곱이라 가정하고 알고리즘을 세웠습니다. $N$이 2의 거듭제곱이 아닌 경우를 풀려면 두 가지 선택지가 있습니다.

  1. $N-1$에 대해 FFT한 결과를 가지고 $N$에 대한 FFT를 해주는 알고리즘 만들기. 이렇게 하면 홀수일 경우 하나 작은 짝수로 바꾸어 풀게 됩니다.
  2. 2의 거듭제곱이 되도록 수열을 늘립니다. 물론 FFT 결과는 다를 겁니다.

1번: 직접 해봅시다.

2번: FFT의 가장 주요한 목표가 합성곱 빠르게 구하기라는 점에서 착안하여, FFT는 못 구해도 합성곱은 구할 수 있게 수열을 바꿔버립니다. 주어진 수열 $\mathbf a$와 $\mathbf b$에 대해 $2N$보다 작지 않은 2의 거듭제곱 중에서 제일 작은 것을 $N’$이라 합시다. 그리고 주기가 $N’$인 복소 수열 $\mathbf a’$과 $\mathbf b’$을 다음과 같이 정의합니다.

\[ \mathbf a’ = \{ a(0), a(1), \cdots, a(N-1), \underbrace{0, 0, \cdots, 0}_{N’-N} \} \]

\[ \mathbf b’ = \{ b(0), b(1), \cdots, b(N-1), \underbrace{0, 0, \cdots, 0}_{N’-2N}, b(0), b(1), \cdots, b(N-1) \} \]

이때 $\mathbf a’$과 $\mathbf b’$의 합성곱의 첫 $N$개 항은 $\mathbf a$와 $\mathbf b$의 합성곱과 똑같습니다.

증명) 합성곱 설명할 때 썼던 그림을 그려보면 자명합니다. (??)

따라서 임의의 $N$에 대해 합성곱을 구해주는 함수를 구현할 수 있습니다.

응용

BOJ 1067:: 이동

https://www.acmicpc.net/problem/1067

딱 보면 합성곱의 냄새가 납니다. 정확히는, $\mathbf b$를 뒤집은 다음 합성곱을 구해야 합니다.

원래 수열이 정수 수열이므로 합성곱도 정수로 나와야 하지만 계산 결과는 복소수입니다. 실수부를 잘 반올림해서 정수로 만들어줍시다. 생각 외로 오차가 크지 않아서 아주 잘 계산됩니다.

빠른 다항식 곱셈 (BOJ 11385:: 씽크스몰)

https://www.acmicpc.net/problem/11385

수열 $\mathbf a$와 $\mathbf b$를 계수로 하는 두 $N-1$차 다항식

\[ \begin{align}
f(x) &= a(0) + a(1) x + a(2) x^2 + \cdots + a(N-1) x^{N-1} \\
g(x) &= b(0) + b(1) x + b(2) x^2 + \cdots + b(N-1) x^{N-1}
\end{align} \]

의 곱은

\[ h(x) = \sum_{k=0}^{N-1} \sum_{l=0}^{k} a(l) b(k-l) x^k + \sum_{k=N}^{2N-2} \sum_{l=k-N+1}^{N-1} a(l) b(k-l) x^k \]

입니다. 이걸 좀 보기 쉽도록 $a(N), a(N+1), \cdots, a(2N-2)$와 $b(N), b(N+1), \cdots, b(2N-2)$를 0으로 정의하면 (즉, $\mathbf a$와 수열 $\mathbf b$에 0을 이어붙여 길이를 늘립니다.)

\[ h(x) = \sum_{k=0}^{2N-1} \sum_{l=0}^{2N-1} a(l) b(k-l) x^k \]

가 됩니다. 합성곱이네요. 핵심은 $\mathbf a$와 $\mathbf b$에 0을 이어붙인다는 건데 길이가 $2N-1$ 이상이 되도록 충분히 많이 이어붙이면 됩니다. 그렇게 해서 나온 길이가 2의 거듭제곱이면 되겠죠. $f$와 $g$의 차수가 달라도 상관 없습니다. 0을 적당히 붙여서 차수가 작은 쪽을 큰 쪽에 맞춰버리면 됩니다.

씽크스몰 문제는 곱의 계수가 $10^{18}$정도 나오기 때문에 실수 오차에 좀 주의해야 합니다. 25-29번 줄을 보면 미리 $\omega_s$의 값을 구하는 부분이 있습니다. 원래 방식대로 1에다가 $\omega$를 계속 곱해나가다간 실수 오차로 틀렸습니다를 받게 됩니다. long long 쓰는 것도 잊지 맙시다.

빠른 정수 곱셈 (BOJ 11576:: 큰 수 곱셈 (2))

https://www.acmicpc.net/problem/15576

이를테면 정수 123456은 다항식 $123x+456$에 $x=1000$을 대입한 것으로 볼 수 있습니다. 이 아이디어를 쓰면 정수를 몇 자리씩 끊어서 다항식으로 취급한 다음 다항식 곱셈을 하고 $x$에 값을 대입해주는 방식으로 정수 곱셈을 할 수 있습니다. 이때 $x$는 10의 거듭제곱이니 그냥 받아올림 하듯이 처리해주면 되죠. 저는 6자리씩 끊었습니다.

BOJ 10531:: Golf Bot

https://www.acmicpc.net/problem/10531

골프 로봇이 거리 $n$만큼 공을 날려보낼 수 있으면 $a(n)=1$, 그렇지 않으면 $a(n)=0$으로 정의합시다. (단, $a(0) = 1$입니다.) 이때

\[ c(n) = \sum_{k=0}^n a(k) a(n-k) \]

라 하면 골프 로봇을 두 번 이하로 써서 $n$에 도달할 수 있으면 $c(n) \neq 0$, 그렇지 않으면 $c(n) = 0$이 됩니다. 식을 변형하기 위해 골프 로봇이 날릴 수 있는 최대 거리를 $D$이라 하고 $a(D+1) = a(D+2) = \cdots = a(2D+1) = 0$라 정의하면

\[ c(n) = \sum_{k=0}^{2D+1} a(k) a(n-k) \]

로 식을 바꿀 수 있습니다. 여기서 수열 $\mathbf a$의 주기는 $2D+2$입니다. 이제 합성곱 꼴이 나왔으므로 FFT를 쓸 수 있습니다.