Up

고속 푸리에 변환

2018년 6월 26일

기저

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

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

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

\begin{gather} \mathbf a \cdot \v_k = c_1 \v_1 \cdot \v_k + c_2 \v_2 \cdot \v_k + \cdots + c_n \v_n \cdot \v_k = c_k \v_k \cdot \v_k \\ \therefore c_k = \frac{ \mathbf a \cdot \v_k }{ \v_k \cdot \v_k } \end{gather}

수열의 내적과 직교

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

\[ \mathbf a \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(i)$는 $\mathbf b$의 $i$번째 항 $b(i)$의 켤레복소수입니다. 벡터의 내적과 매우 유사한 정의인데, 마찬가지로 내적이 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) \]

로 정의하면 $\boldsymbol\phi_k$는 주기가 $N$이고(삼각함수의 성질) $\boldsymbol\phi_0, \boldsymbol\phi_1, \cdots, \boldsymbol\phi_{N−1}$은 서로 직교합니다. 증명하려면 서로 다른 $k$와 $l$에 대해 $\boldsymbol\phi_k\cdot\boldsymbol\phi_l$을 계산해보면 됩니다.

\begin{aligned} \boldsymbol\phi_k \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{aligned}

한편 자기 자신과의 내적은 $N$입니다.

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

크로네커 델타 기호를 사용해 두 경우를 하나로 표현할 수 있습니다.

\[ \boldsymbol\phi_k \cdot \boldsymbol\phi_l = N\delta_{kl} \]

이산 푸리에 변환

$\boldsymbol\phi_0, \boldsymbol\phi_1, \cdots, \boldsymbol\phi_{N−1}$은 영수열(모든 항이 0인 수열)이 아니면서 서로 직교하므로, 임의의 수열 $\mathbf A$는 이들의 일차결합으로 나타낼 수 있습니다.

\[ \mathbf A = c_0\boldsymbol\phi_0 + c_1\boldsymbol\phi_1 + \cdots + c_{N-1}\boldsymbol\phi_{N-1} \]

계수 $c_0,\cdots,c_{N-1}$을 각 항으로 갖는 수열 $\mathbf a$를 생각합시다. 이때 $\mathbf A$는 $\mathbf a$의 이산 푸리에 변환(descrete Fourier transform, DFT)라고 정의하고 기호로는 $\mathcal F\{\mathbf a\}$로 씁니다. 즉,

\[ \mathcal F\{\mathbf a\}(n) = A(n) = \sum_{k=0}^{N-1} a(k) \phi_k(n) \label{eq:dft} \]

역으로 $\mathbf a$는 $\mathbf A$의 역변환이 되고, $\mathbf a = \mathcal F^{-1}\{\mathbf A\}$로 표기합니다. $\boldsymbol\phi_k$가 주기 $N$인 수열이므로 $\mathbf A$ 또한 $\mathbf a$와 마찬가지로 주기 $N$인 수열이 됩니다.

역변환 식은 식 \eqref{eq:dft}의 양변에 $\bar\phi_l(n)$를 곱하여 구할 수 있습니다.

\begin{align} A(n)\bar\phi_l(n) &= \sum_{k=0}^{N-1} a(k)\phi_k(n)\bar\phi_l(n) \nonumber \\ \sum_{n=0}^{N-1} A(n)\bar\phi_l(n) &= \sum_{k=0}^{N-1} \left[ a(k) \sum_{n=0}^{N-1} \phi_k(n)\bar\phi_l(n) \right] = Na(l) \nonumber \\ \therefore a(n) &= \frac1N \sum_{k=0}^{N-1}A(k)\bar\phi_n(k) = \frac1N \sum_{k=0}^{N-1}A(k)\bar\phi_k(n) \label{eq:inv-dft} \end{align}

식 \eqref{eq:dft}과 식 \eqref{eq:inv-dft}을 비교하면 거의 동일한 꼴입니다. 따라서 이산 푸리에 변환은 역변환에 대해 크게 신경쓸 일이 없습니다.

합성곱

수열 $\mathbf a$와 $\mathbf b$의 이산 푸리에 변환을 항별로 곱하여 나온 수열을 다시 역변환하여 수열 $\mathbf c$를 얻었습니다. $\mathbf a$, $\mathbf b$, $\mathbf c$ 사이의 관계는 어떻게 될까요?

\begin{gather} \begin{aligned} \mathcal F\{\mathbf c\}(n) &= \mathcal F\{\mathbf a\}(n) \mathcal F\{\mathbf b\}(n) \\ &= \sum_{k=0}^{N-1} a(k)\phi_k(n) \sum_{l=0}^{N-1} b(l)\phi_l(n) \\ &= \sum_{k=0}^{N-1}\sum_{l=0}^{N-1} a(k)b(l) \exp\left[ \frac{2\pi(k+l)}N ni \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{aligned} \\ \therefore c(n) = \sum_{k=0}^{N-1} a(k)b(n-k) \label{eq:convolution} \end{gather}

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

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

이산 푸리에 변환에서 합성곱이 갑자기 나오는 이유는 이산 푸리에 변환의 주요한 목적 중 하나가 합성곱을 쉽게 구하기 위해서이기 때문입니다. 합성곱은 의외로 여기저기서 많이 나타나는 수열 연산이므로 이를 빠르게 구하는 건 중요합니다. 그래서 두 수열의 합성곱을 구하고 싶을 때 식 \eqref{eq:convolution}을 직접 계산하는 것이 아니라 각각을 이산 푸리에 변환하고, 항별로 곱한 다음, 역변환을 취하자는 말입니다. 그런데 이게 과연 더 빠를까요?

고속 푸리에 변환

본격적으로 이산 푸리에 변환을 구해봅시다. 식 \eqref{eq:dft}에서 $A(n)$ 하나를 구하는 데 $O(N)$만큼의 시간이 걸리므로 $\mathbf A$ 전체를 구하는 데 $O(N^2)$이면 충분합니다.

당연하게도, 정의 그대로 구현하면 끔찍하게 느립니다! 그래서 아주 옛날부터 이산 푸리에 변환을 $O(N\log N)$에 구하는 다양한 알고리즘들이 제시되었습니다. 이런 알고리즘들을 통틀어 고속 푸리에 변환(fast Fourier transform, FFT)이라 부릅니다. 이 중 가장 일반적이고 흔히 사용되는 것은 제임스 쿨리와 존 튜키가 1965년 발표한 쿨리-튜키 알고리즘입니다.

식 \eqref{eq:dft}을 다시 쓰면,

\[ 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}Nki\right) = \sum_{k=0}^{N-1} a(k)\omega_N^{kn} \]

여기서 $\omega_N = \exp(2\pi i/N)$은 1의 $N$제곱근입니다. $A(n)$은 합으로 정의되니, $k$가 짝수일 때와 홀수일 때로 나눠봅시다. 편의상 $N$은 짝수라고 가정합니다.

\[ A(n) = \sum_{k=0}^{N/2-1} a(2k)\omega_N^{2kn} + \sum_{k=0}^{N/2-1} a(2k+1)\omega_N^{(2k+1)n} = \sum_{k=0}^{N/2-1} a(2k)\omega_N^{2kn} + \omega_N^n \sum_{k=0}^{N/2-1} a(2k+1)\omega_N^{2kn} \label{eq:dft-div} \]

$\mathbf a$의 짝수 번째 항을 모아 만든 수열을 $\mathbf a_\mathrm{even}$, 홀수 번째 항을 모아 만든 수열을 $\mathbf a_\mathrm{odd}$라 합시다. 즉,

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

정의상 $\mathbf a_\mathrm{even}$과 $\mathbf a_\mathrm{odd}$는 주기가 $N/2$인 수열입니다. 따라서 이들의 푸리에 변환은

\begin{align} \mathcal F\{\mathbf a_\mathrm{even}\}(n) = \mathbf A_\mathrm{even}(n) &= \sum_{k=0}^{N/2-1} a_\mathrm{even}(k)\omega_{N/2}^{kn} = \sum_{k=0}^{N/2-1} a(2k)\omega_N^{2kn} \\ \mathcal F\{\mathbf a_\mathrm{odd}\}(n) = \mathbf A_\mathrm{odd}(n) &= \sum_{k=0}^{N/2-1} a_\mathrm{odd}(k)\omega_{N/2}^{kn} = \sum_{k=0}^{N/2-1} a(2k+1)\omega_N^{2kn} \end{align}

이를 식 \eqref{eq:dft-div}에 대입하면 $\mathbf A$를 $\mathbf A_\mathrm{even}$과 $\mathbf A_\mathrm{odd}$로 표현할 수 있습니다.

\[ A(n) = A_\mathrm{even}(n) + \omega_N^n A_\mathrm{odd}(n) \label{eq:fft-dnc-1} \]

한편 $\mathbf A_\mathrm{even}$과 $\mathbf A_\mathrm{odd}$도 주기가 $N/2$이므로

\[ A(n+N/2) = A_\mathrm{even}(n+N/2) + \omega_N^{n+N/2} A_\mathrm{odd}(n+N/2) = A_\mathrm{even}(n) - \omega_N^n A_\mathrm{odd}(n) \label{eq:fft-dnc-2} \]

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

  1. 주어진 수열 $\mathbf a$를 $\mathbf a_\mathrm{even}$과 $\mathbf a_\mathrm{odd}$로 나눈다.
  2. 나눈 각각에 대해 이산 푸리에 변환을 계산한다.
  3. 식 \eqref{eq:fft-dnc-1} 및 \eqref{eq:fft-dnc-2}에 따라 $\mathbf a$의 이산 푸리에 변환을 계산한다.

다만 $N$이 짝수라고 가정하고 만든 알고리즘이므로 $\mathbf a_\mathrm{even}$과 $\mathbf a_\mathrm{odd}$의 이산 푸리에 변환을 같은 방법으로 구하기 위해서는 $N/2$, $N/4$, $...$ 또한 짝수여야 합니다. 고로 이 알고리즘은 $N$이 2의 거듭제곱일 때만 작동하는데, 이 정도만 해도 당장은 충분합니다.

이 방식을 파이썬으로 구현하면 아래와 같습니다. $N=1$인 경우에는 이산 푸리에 변환이 자기 자신입니다.

from cmath import exp, pi

def fft(a):
    N = len(a)
    if N == 1: return a
    A_even = fft(a[0::2])
    A_odd = fft(a[1::2])
    w_N = [exp(2j*pi*n/N) for n in range(N//2)]
    return [A_even[n] + w_N[n] * A_odd[n] for n in range(N//2)] + \
        [A_even[n] - w_N[n] * A_odd[n] for n in range(N//2)]

재귀 풀기

$N=8$에 대해 FFT를 다음과 같이 회로처럼 도식화하여 그려봅시다. 화살표에 쓰인 수는 곱셈을 의미합니다. 식 \eqref{eq:fft-dnc-1}과 \eqref{eq:fft-dnc-2}은 좌우 두 회로가 동등하다는 것을 말해줍니다. 다만 $\mathbf a$의 짝수 번째 항은 위쪽, 홀수 번째 항은 아래쪽으로 재배열된 것만 주의합시다.

마찬가지로 $N=4$에 대해 FFT를 해주는 회로는 $N=2$에 대한 FFT 회로 두 개를 써서 만들 수 있고, 이렇게 재귀적으로 그리면 다음 회로가 나옵니다. 교차하는 부분이 나비처럼 생겨서 흔히 버터플라이 다이어그램이라고 부릅니다.

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

\begin{align*} a_\text{rearr}(000) &= a(000) \\ a_\text{rearr}(001) &= a(100) \\ a_\text{rearr}(010) &= a(010) \\ a_\text{rearr}(011) &= a(110) \\ a_\text{rearr}(100) &= a(001) \\ a_\text{rearr}(101) &= a(101) \\ a_\text{rearr}(110) &= a(011) \\ a_\text{rearr}(111) &= a(111) \end{align*}

이로부터 다음을 유추할 수 있습니다.

$N$이 2의 거듭제곱이면 \[ a_\text{rearr}(n) = a(n\text{을 } \log N \text{비트로 나타낸 뒤 반전한 수}) \]
먼저 $N=1$일 때 자명하게 성립합니다. 그리고 $N=2^k$일 때 성립한다고 가정합시다. ($k\ge0$) $2N$에 대한 FFT를 하려면 먼저 $\mathbf a_\mathrm{even}$과 $\mathbf a_\text{odd}$로 나눠야죠. 정의에 의하여 $0\le n\le N-1$에 대해
\begin{align} a_\text{even} ( \underbrace{\bits} _{\substack{\text{bit representation of } n \\ (\text{total } k \text{ bits}) }} ) &= a [ 2 ( \underbrace{\bits}_{\text{bit repr. of }n} ) ] = a ( \underbrace{\bits}_{\text{bit repr. of }n} 0 ) \\ a_\text{odd} ( \underbrace{\bits} _{\substack{\text{bit representation of } n \\ (\text{total } k \text{ bits}) }} ) &= a [ 2 ( \underbrace{\bits}_{\text{bit repr. of }n} ) + 1 ] = a ( \underbrace{\bits}_{\text{bit repr. of }n} 1 ) \end{align}

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

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

같은 방법으로

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

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

\begin{align} a_\text{rearr}(0\underbrace{\bits}_{\text{bit repr. of }n}) &= a_\text{even, rearr}(\underbrace{\bits}_{\text{bit repr. of }n}) = a(\underbrace{\bits}_{\substack{\text{reversed} \\ \text{bit repr. of }n}}0) \\ a_\text{rearr}(1\underbrace{\bits}_{\text{bit repr. of }n}) &= a_\text{odd, rearr}(\underbrace{\bits}_{\text{bit repr. of }n}) = a(\underbrace{\bits}_{\substack{\text{reversed} \\ \text{bit repr. of }n}}1) \end{align}

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

정수 $n$을 $k$비트로 나타낸 뒤 뒤집는 연산은 여러 구현법이 있습니다. 여기서는 다음 구현법을 쓰겠습니다. (출처)

inline unsigned bitreverse(const unsigned n, const unsigned k) {
    unsigned r, i;
    for (r = 0, i = 0; i < k; ++i)
        r |= ((n >> i) & 1) << (k - i - 1);
    return r;
}

남은 건 정말로 구현 뿐입니다. 그 전에 잠시 역변환을 어떻게 구현할 지 생각해봅시다. 역변환은 기저로 $\bar{\boldsymbol\phi}_0$, $\cdots$, $\bar{\boldsymbol\phi}_{N-1}$을 쓴다고 했으므로 위 과정을 그대로 해주되 $\omega_N$ 대신 $\bar\omega_n=1/\omega_N$을 넣어주면 됩니다. 그리고 최종 결과를 $N$으로 나눠주면 끝이죠.

다음은 C++로 구현한 FFT 코드입니다. is_reverse가 참이면 역변환을 계산합니다.

using namespace std;
using cdbl = complex<double>;

const double PI = acos(-1);

void fft(vector<cdbl> &a, bool is_reverse=false) {
    const unsigned n = a.size(), k = __builtin_ctz(n);
    unsigned s, i, j;
    for (i = 0; i < n; i++) {
        j = bitreverse(i, k);
        if (i < j)
            swap(a[i], a[j]);
    }
    for (s = 2; s <= n; s *= 2) {
        double t = 2*PI/s * (is_reverse? -1 : 1);
        cdbl ws(cos(t), sin(t));
        for (i = 0; i < n; i += s) {
            cdbl w(1);
            for (j = 0; j < s/2; j++) {
                cdbl tmp = a[i + j + s/2] * w;
                a[i + j + s/2] = a[i + j] - tmp;
                a[i + j] += tmp;
                w *= ws;
            }
        }
    }
    if (is_reverse)
        for (i = 0; i < n; i++)
            a[i] /= n;
}

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

분할 정복의 특성상 지금까지는 $N$이 2의 거듭제곱인 경우만 다루었습니다. 아닌 경우를 해결하는 방법은 두 가지가 있습니다.

  1. $N$이 홀수인 경우도 분할 정복 가능하게 알고리즘을 잘 바꿉니다.
  2. 수열을 잘 늘려서 주기가 2의 거듭제곱이 되게 만듭니다. FFT 결과는 다르겠지만, FFT의 목표가 합성곱 구하기라는 점에 착안하여 합성곱은 잘 구할 수 있게 수열을 늘리면 됩니다.

2번을 구체적으로 설명하면, 주기가 $N$인 두 수열 $\mathbf a$와 $\mathbf b$의 합성곱을 구한다고 해봅시다. $2N$ 이상인 2의 거듭제곱 중 제일 작은 것을 $N'$이라 하고, 주기가 $N'$인 수열 $\mathbf a'$과 $\mathbf b'$을 다음과 같이 정의합니다.

\begin{align} \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) \} \end{align}

이때 $\mathbf a' * \mathbf b'$의 첫 $N$개 항은 $\mathbf a * \mathbf b$와 같습니다. 합성곱 설명할 때 썼던 그림을 그려보면 자명합니다. 따라서 임의의 $N$에 대해 합성곱을 구해주는 함수를 구현할 수 있습니다.

void convolution(vector<cdbl> &a, vector<cdbl> &b) {
    const unsigned n = a.size();
    unsigned np = n, i;

    if (n & (n-1)) {
        for (np = 1; np < 2*n; np *= 2);
        a.resize(np);
        b.resize(np);
        for (i = 0; i < n; i++)
            b[np-n+i] = b[i];
    }

    fft(a);
    fft(b);
    for (i = 0; i < np; i++)
        a[i] *= b[i];
    fft(a, true);

    if (n & (n-1)) {
        a.resize(n);
        b.resize(n);
    }
}

합성곱의 결과는 a에 저장됩니다.

응용

빠른 다항식 곱셈

수열 $\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 \label{eq:polynomial-product} \]

이걸 좀 알아보기 쉽게 길이 $N'$인 새로운 수열 $\mathbf a'$과 $\mathbf b'$을 정의합시다. 여기서 $N'$은 $2N-1$ 이상인 적당한 2의 거듭제곱입니다.

\begin{align} a'(n) &= \begin{cases} a(n), & 0 \le n \le N-1 \\ 0, & N \le n \le N'-1 \end{cases} \\ b'(n) &= \begin{cases} b(n), & 0 \le n \le N-1 \\ 0, & N \le n \le N'-1 \end{cases} \end{align}

그러면 식 \eqref{eq:polynomial-product}은

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

따라서 다항식의 곱의 계수는 $\mathbf a'$과 $\mathbf b'$의 합성곱으로 구할 수 있습니다. 두 다항식의 차수가 다른 경우엔 작은 쪽을 큰 쪽에 맞춰주면 됩니다.

BOJ 11385:: 싱크스몰

하라는 대로 다항식 곱셈을 하면 됩니다. 단, 계수가 1018정도까지 가능하므로 실수오차에 좀 주의해야 합니다. 위 FFT 코드에서는 $\omega$를 여러 번 곱해서 $\omega$의 거듭제곱을 구했는데 여기선 25-29번째 줄에서 $\omega$의 거듭제곱들을 미리 계산하고 있습니다. 그리고 당연히 long long형을 사용해야 한다는 점도 잊지 맙시다.

#include<bits/stdc++.h>

using namespace std;
using cdbl = complex<long double>;

const long double PI = acos(-1.L);

inline unsigned bitreverse(const unsigned n, const unsigned k) {
    unsigned r, i;
    for (r = 0, i = 0; i < k; ++i)
        r |= ((n >> i) & 1) << (k - i - 1);
    return r;
}

void fft(vector<cdbl> &a, bool is_reverse=false) {
    const unsigned n = a.size(), k = __builtin_ctz(n);
    unsigned s, i, j;
    for (i = 0; i < n; i++) {
        j = bitreverse(i, k);
        if (i < j)
            swap(a[i], a[j]);
    }
    for (s = 2; s <= n; s *= 2) {
        vector<cdbl> w(s/2);
        for (i = 0; i < s/2; i++) {
            long double t = 2*PI*i/s * (is_reverse? -1 : 1);
            w[i] = cdbl(cos(t), sin(t));
        }
        for (i = 0; i < n; i += s) {
            for (j = 0; j < s/2; j++) {
                cdbl tmp = a[i + j + s/2] * w[j];
                a[i + j + s/2] = a[i + j] - tmp;
                a[i + j] += tmp;
            }
        }
    }
    if (is_reverse)
        for (i = 0; i < n; i++)
            a[i] /= n;
}

void convolution(vector<cdbl> &a, vector<cdbl> &b) {
    const unsigned n = a.size();
    fft(a);
    fft(b);
    for (unsigned i = 0; i < n; i++)
        a[i] *= b[i];
    fft(a, true);
}


int main(void) {
    ios_base::sync_with_stdio(false);
    cin.tie(nullptr);
    cout.tie(nullptr);

    unsigned n, m, np = 1;
    cin >> n >> m;
    n++;
    m++;
    for (; np < 2*max(n, m)-1; np *= 2);

    vector<cdbl> a(np, 0), b(np, 0);
    for (unsigned i = 0; i < n; i++)
        cin >> a[i];
    for (unsigned i = 0; i < m; i++)
        cin >> b[i];

    convolution(a, b);

    unsigned long long res = 0;
    for (unsigned i = 0; i < n+m-1; i++)
        res ^= (unsigned long long)(a[i].real() + (a[i].real() > 0? 0.5L : -0.5L));

    cout << res;

    return 0;
}

빠른 정수 곱셈

예를 들어 정수 123456은 다항식 $123x+456$에 $x=1000$을 대입한 것으로 볼 수 있습니다. 따라서 정수를 몇 자리씩 끊어서 다항식으로 만든 다음 곱하고 $x$에 값을 대입하면 정수 곱셈을 빠르게 할 수 있습니다.

BOJ 11576:: 큰 수 곱셈 (2)

6자리씩 끊어서 계산했습니다. 이러면 다항식의 곱의 계수가 대략 3×1017까지이므로 long long형을 써서 충분히 계산 가능합니다.

BOJ 1067:: 이동

$Y$를 뒤집은 다음 $X$와 합성곱을 하면 됩니다. 정답은 합성곱의 항 중 최댓값입니다.

BOJ 10531:: Golf Bot

골프 로봇이 거리 $n$만큼 공을 날려보낼 수 있으면 $a(n)=1$, 아니면 $a(n)=0$으로 정의합시다. 단, $a(0)=1$로 정의하고 골프 로봇이 보낼 수 있는 최대 거리보다 큰 $n$에 대해서는 $a(n)$이 0입니다. 이때

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

라 하면 골프 로봇을 두 번 이하로 써서 거리 $n$만큼 날려보낼 수 있는 경우 $c(n)\neq0$, 불가능한 경우 $c(n)=0$입니다. 이제 $\mathbf a$를 주기 $2D+2$인 수열로 생각하고 자기 자신과의 합성곱을 계산합니다.

int main(void) {
    ios_base::sync_with_stdio(false);
    cin.tie(nullptr);
    cout.tie(nullptr);

    int n, m, k, d, cnt = 0;
    vector<cdbl> a(524288, 0);

    cin >> n;
    for (int i = 0; i < n; i++) {
        cin >> k;
        a[k] = 1;
    }
    a[0] = 1;

    fft(a);
    for (int i = 0; i < 524288; i++)
        a[i] *= a[i];
    fft(a, true);

    cin >> m;
    for (int i = 0; i < m; i++) {
        cin >> d;
        if (int(a[d].real() + (a[d].real() > 0? 0.5 : -0.5)) != 0)
            cnt++;
    }

    cout << cnt;

    return 0;
}