기저
서로 수직인 (영벡터가 아닌) $n$차원 벡터 $\v_1, \v_2, \cdots, \v_n$을 생각해봅시다. 임의의 $n$차원 벡터 $\mathbf a$는 항상 이들의 일차 결합으로 나타낼 수 있습니다.
이때 계수 $c_k$는 양변에 $\v_k$를 내적하여 구할 수 있습니다.
수열의 내적과 직교
내적은 벡터에서만 정의하는 것이 아닙니다. 아주 다양한 대상에 내적을 정의하여 해당 대상을 벡터처럼 다룰 수 있습니다. 주기가 $N$인 복소 수열(앞으로 별 말이 없으면 수열은 항상 주기가 $N$인 복소 수열을 말합니다.) $\mathbf a$와 $\mathbf b$의 내적을 아래와 같이 정의합시다.
여기서 $\bar{b}(i)$는 $\mathbf b$의 $i$번째 항 $b(i)$의 켤레복소수입니다. 벡터의 내적과 매우 유사한 정의인데, 마찬가지로 내적이 0인 두 수열은 직교한다고 말합니다.
이제 정수 $k$에 대해 수열 $\boldsymbol\phi_k$를
로 정의하면 $\boldsymbol\phi_k$는 주기가 $N$이고(삼각함수의 성질) $\boldsymbol\phi_0, \boldsymbol\phi_1, \cdots, \boldsymbol\phi_{N−1}$은 서로 직교합니다. 증명하려면 서로 다른 $k$와 $l$에 대해 $\boldsymbol\phi_k\cdot\boldsymbol\phi_l$을 계산해보면 됩니다.
한편 자기 자신과의 내적은 $N$입니다.
크로네커 델타 기호를 사용해 두 경우를 하나로 표현할 수 있습니다.
이산 푸리에 변환
$\boldsymbol\phi_0, \boldsymbol\phi_1, \cdots, \boldsymbol\phi_{N−1}$은 영수열(모든 항이 0인 수열)이 아니면서 서로 직교하므로, 임의의 수열 $\mathbf A$는 이들의 일차결합으로 나타낼 수 있습니다.
계수 $c_0,\cdots,c_{N-1}$을 각 항으로 갖는 수열 $\mathbf a$를 생각합시다. 이때 $\mathbf A$는 $\mathbf a$의 이산 푸리에 변환(descrete Fourier transform, DFT)라고 정의하고 기호로는 $\mathcal F{\mathbf a}$로 씁니다. 즉,
역으로 $\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)$를 곱하여 구할 수 있습니다.
식 \eqref{eq:dft}과 식 \eqref{eq:inv-dft}을 비교하면 거의 동일한 꼴입니다. 따라서 이산 푸리에 변환은 역변환에 대해 크게 신경쓸 일이 없습니다.
합성곱
수열 $\mathbf a$와 $\mathbf b$의 이산 푸리에 변환을 항별로 곱하여 나온 수열을 다시 역변환하여 수열 $\mathbf c$를 얻었습니다. $\mathbf a$, $\mathbf b$, $\mathbf c$ 사이의 관계는 어떻게 될까요?
(세 번째 줄에서 네 번째 줄로 넘어갈 때에는 삼각함수의 주기성을 사용합니다. 조금 유도하기 어렵습니다.) 여기서 $\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}을 다시 쓰면,
여기서 $\omega_N = \exp(2\pi i/N)$은 1의 $N$제곱근입니다. $A(n)$은 합으로 정의되니, $k$가 짝수일 때와 홀수일 때로 나눠봅시다. 편의상 $N$은 짝수라고 가정합니다.
$\mathbf a$의 짝수 번째 항을 모아 만든 수열을 $\mathbf a_\mathrm{even}$, 홀수 번째 항을 모아 만든 수열을 $\mathbf a_\mathrm{odd}$라 합시다. 즉,
정의상 $\mathbf a_\mathrm{even}$과 $\mathbf a_\mathrm{odd}$는 주기가 $N/2$인 수열입니다. 따라서 이들의 푸리에 변환은
이를 식 \eqref{eq:dft-div}에 대입하면 $\mathbf A$를 $\mathbf A_\mathrm{even}$과 $\mathbf A_\mathrm{odd}$로 표현할 수 있습니다.
한편 $\mathbf A_\mathrm{even}$과 $\mathbf A_\mathrm{odd}$도 주기가 $N/2$이므로
따라서 이산 푸리에 변환을 재귀적으로 계산하려면,
- 주어진 수열 $\mathbf a$를 $\mathbf a_\mathrm{even}$과 $\mathbf a_\mathrm{odd}$로 나눈다.
- 나눈 각각에 대해 이산 푸리에 변환을 계산한다.
- 식 \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}$이라 하면
이로부터 다음을 유추할 수 있습니다.
$N$이 2의 거듭제곱이면
먼저 $N=1$일 때 자명하게 성립합니다. 그리고 $N=2^k$일 때 성립한다고 가정합시다. ($k\ge0$) $2N$에 대한 FFT를 하려면 먼저 $\mathbf a_\mathrm{even}$과 $\mathbf a_\text{odd}$로 나눠야죠. 정의에 의하여 $0\le n\le N-1$에 대해
$\mathbf a_\text{even}$의 FFT를 계산하기 위해 재배열 한 것을 $\mathbf a_\mathrm{even,rearr}$이라 하면 가정에 따라
같은 방법으로
$\mathbf a_\text{rearr}$은 $\mathbf a_\text{even, rearr}$와 $\mathbf a_\text{odd, rearr}$을 이어붙인 것이므로
가 되어 $2N=2^{k+1}$일 때에도 성립합니다.
정수 $n$을 $k$비트로 나타낸 뒤 뒤집는 연산은 여러 구현법이 있습니다. 여기서는 다음 구현법을 쓰겠습니다. 출처
unsigned bitreverse(unsigned n, 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의 거듭제곱인 경우만 다루었습니다. 아닌 경우를 해결하는 방법은 두 가지가 있습니다.
- $N$이 홀수인 경우도 분할 정복 가능하게 알고리즘을 잘 바꿉니다.
- 수열을 잘 늘려서 주기가 2의 거듭제곱이 되게 만듭니다. FFT 결과는 다르겠지만, FFT의 목표가 합성곱 구하기라는 점에 착안하여 합성곱은 잘 구할 수 있게 수열을 늘리면 됩니다.
2번을 구체적으로 설명하면, 주기가 $N$인 두 수열 $\mathbf a$와 $\mathbf b$의 합성곱을 구한다고 해봅시다. $2N$ 이상인 2의 거듭제곱 중 제일 작은 것을 $N’$이라 하고, 주기가 $N’$인 수열 $\mathbf a’$과 $\mathbf b’$을 다음과 같이 정의합니다.
이때 $\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$차 다항식
의 곱은 다음과 같습니다.
이걸 좀 알아보기 쉽게 길이 $N’$인 새로운 수열 $\mathbf a’$과 $\mathbf b’$을 정의합시다. 여기서 $N’$은 $2N-1$ 이상인 적당한 2의 거듭제곱입니다.
그러면 식 \eqref{eq:polynomial-product}은
따라서 다항식의 곱의 계수는 $\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입니다. 이때
라 하면 골프 로봇을 두 번 이하로 써서 거리 $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;
}