이산 푸리에 변환이란 주기가 \(N\)인 복소수열 \(x_{0},x_{1},\dots,x_{N-1}\)에 대하여 다음과 같이 정의되는 주기가 \(N\)인 복소수열 \(X_{0},X_{1},\dots,X_{N-1}\)로 바꾸는 변환을 말한다.
\[X_{k}=\sum_{n=0}^{N-1}e^{-2\pi ikn/N}x_{n}\]
여기서 \(e^{iz}=\cos z+i\sin z\)이며, 이로부터 \(e^{2\pi i}=e^{-2\pi i}=1\)임을 쉽게 알 수 있다.
~
이산 푸리에 역변환
또한 이렇게 변환된 복소수열 \(X_{k}\)는 다음과 같은 과정을 통하여 복소수열 \(x_{n}\)으로 역변환될 수 있다.
\[x_{n}=\frac{1}{N}\sum_{k=0}^{N-1}e^{2\pi ikn/N}X_{k}\]
~
이산 푸리에 역변환의 증명
\(X_{k}\)의 정의대로 식을 풀어보면 다음과 같다.
\[\sum_{k=0}^{N-1}e^{2\pi ikn/N}X_{k}=\sum_{k=0}^{N-1}\sum_{j=0}^{N-1}e^{2\pi ik(n-j)/N}x_{j}=\sum_{j=0}^{N-1}x_{j}\sum_{k=0}^{N-1}e^{2\pi ik(n-j)/N}\]
\(j=n\)인 경우에는 다음처럼 된다.
\[\sum_{k=0}^{N-1}e^{2\pi ik(n-j)/N}=\sum_{k=0}^{N-1}1=N\]
\(j\neq n\)인 경우에는 다음처럼 된다.
\[\sum_{k=0}^{N-1}e^{2\pi ik(n-j)/N}=\frac{e^{2\pi i(n-j)}-1}{e^{2\pi i(n-j)/N}-1}=0\]
따라서 다음이 성립한다.
\[\frac{1}{N}\sum_{k=0}^{N-1}e^{2\pi ikn/N}X_{k}=\frac{1}{N}Nx_{n}=x_{n}\]
~
이산 푸리에 변환과 Convolution
주기가 \(N\)인 수열 \(a_{n}\)과 \(b_{n}\)을 각각 이산 푸리에 변환한 수열을 \(A_{n}\)과 \(B_{n}\)이라 하자.
또한 새로운 수열 \(C_{n}=A_{n}B_{n}\)을 정의하고, 이를 이산 푸리에 역변환한 수열을 \(c_{n}\)이라 하면 다음과 같은 식이 성립한다.
\[c_{n}=\sum_{k=0}^{N-1}a_{k}b_{n-k}\]
~
이산 푸리에 변환과 Convolution의 증명
\(C_{n}\)의 정의대로 식을 풀어보면 다음과 같다.
\[C_{n}=A_{n}B_{n}=\sum_{j=0}^{N-1}\sum_{k=0}^{N-1}e^{-2\pi i(j+k)n/N}a_{k}b_{j}=\sum_{j=0}^{N-1}e^{-2\pi ijn/N}\sum_{k=0}^{N-1}a_{k}b_{j-k}\]
따라서 이를 역변환한 수열 \(c_{n}\)에 대한 식이 성립한다.
~
고속 푸리에 변환
이산 푸리에 변환을 단순한 반복문으로 계산하면 \(O(N^{2})\)의 시간복잡도를 가지지만, 이를 \(O(N\log N)\)의 시간복잡도에 행할 수 있는 방법이 고속 푸리에 변환이다.
고속 푸리에 변환에는 다양한 알고리즘이 있지만, 여기서는 \(N=2^{M}\) 꼴의 푸리에 변환에 최적화된 Cooley-Tukey 알고리즘에 대해 설명한다.
\(X_{k}\)를 다음처럼 고쳐 적을 수 있다.
\[X_{k}=\sum_{n=0}^{N-1}e^{-2\pi ikn/N}x_{n}=\sum_{n=0}^{N/2-1}e^{-2\pi ik(2n)/N}x_{2n}+\sum_{n=0}^{N/2-1}e^{-2\pi ik(2n+1)/N}x_{2n+1}=\sum_{n=0}^{N/2-1}e^{-2\pi ikn/(N/2)}x_{2n}+e^{-2\pi ik/N}\sum_{n=0}^{N/2-1}e^{-2\pi ikn/(N/2)}x_{2n+1}\]
이제 수열 \(Y_{k}\)와 \(Z_{k}\)를 다음처럼 정의하자.
\[Y_{k}=\sum_{n=0}^{N/2-1}e^{-2\pi ikn/(N/2)}x_{2n}\]
\[Z_{k}=\sum_{n=0}^{N/2-1}e^{-2\pi ikn/(N/2)}x_{2n+1}\]
수열 \(Y_{k}\)와 \(Z_{k}\)를 계산할 수 있다면, 수열 \(Z_{k}\)는 다음과 같은 방법으로 \(O(N)\)의 시간복잡도로 계산할 수 있다.
\[X_{k}=Y_{k}+e^{-2\pi ik/N}Z_{k}\]
\[X_{k+N/2}=Y_{k}-e^{-2\pi ik/N}Z_{k}\]
그런데, 수열 \(Y_{k}\)와 \(Z_{k}\)는 수열 \(x_{2n}\)과 \(x_{2n+1}\)에 대한 이산 푸리에 변환이므로, 분할정복을 통해 총 \(O(N\log N)\)의 시간복잡도로 이산 푸리에 변환을 구할 수 있다.
#include <cstdio>
#include <cmath>
#include <complex>
const int SZ = 20, N = 1 << SZ;
const double PI = acos(-1);
using namespace std;
int Rev(int x) {
int i, r = 0;
for (i = 0; i < SZ; i++) {
r = r << 1 | x & 1;
x >>= 1;
}
return r;
}
void FFT(complex<double> *a, bool f) {
complex<double> x, y, z;
double w;
int i, j, k;
for (i = 0; i < N; i++) {
j = Rev(i);
if (i < j) {
z = a[i];
a[i] = a[j];
a[j] = z;
}
}
for (i = 1; i < N; i <<= 1) {
w = PI / i;
if (f) w = -w;
x = complex<double>(cos(w), sin(w));
for (j = 0; j < N; j += i << 1) {
y = complex<double>(1);
for (k = 0; k < i; k++) {
z = a[i + j + k] * y;
a[i + j + k] = a[j + k] - z;
a[j + k] += z;
y *= x;
}
}
}
if (f) for (i = 0; i < N; i++) a[i] /= N;
}
complex<double> X[N];
int main() {
double x;
int i, n;
scanf("%d", &n);
for (i = 0; i <= n; i++) {
scanf("%lf", &x);
X[i] = complex<double>(x, 0);
}
FFT(X, false);
for (i = 0; i < N; i++) X[i] *= X[i];
FFT(X, true);
for (i = 0; i <= n + n; i++) printf("%.0f ", X[i].real());
}
~
일반화된 이산 푸리에 변환
지금까지는 함수 \(\mathbb{Z}_{N}\rightarrow\mathbb{C}\)에 대한 이산 푸리에 변환에 대해서만 살펴보았다.
한편 이산 푸리에 변환에 필요한 연산은 덧셈과 복소수곱셈뿐이므로 이 둘이 정의된 모든 집합에 대해서 이산 푸리에 변환이 가능하다.
예를 들어, 함수 \(\mathbb{Z}_{M}\rightarrow\mathbb{C}\)는 덧셈과 복소수곱셈이 가능하므로 \(\mathbb{Z}_{N}\rightarrow\mathbb{Z}_{M}\rightarrow\mathbb{C}\)에 대한 이산 푸리에 변환이 가능하다.
이를 통해 두 \(N\times M\) 행렬 \(A_{i,j}\)와 \(B_{i,j}\)에 대하여 다음을 만족하는 행렬 \(C_{i,j}\)를 구할 수 있다.
\[C_{i,j}=\sum_{k=0}^{N-1}\sum_{l=0}^{M-1}A_{k,l}B_{i-k,j-l}\]
template을 사용하여 구현하면 다음과 같다.
#include <cstdio>
#include <cmath>
#include <complex>
const int SZ = 10, N = 1 << SZ;
const double PI = acos(-1);
using namespace std;
int Rev(int x) {
int i, r = 0;
for (i = 0; i < SZ; i++) {
r = r << 1 | x & 1;
x >>= 1;
}
return r;
}
template <typename T>
struct Func {
T a[N];
T &operator [](int x) {
return a[x];
}
Func<T> operator +(Func<T> &x) {
Func<T> r;
int i;
for (i = 0; i < N; i++) r.a[i] = a[i] + x[i];
return r;
}
Func<T> operator -(Func<T> &x) {
Func<T> r;
int i;
for (i = 0; i < N; i++) r.a[i] = a[i] - x[i];
return r;
}
Func<T> operator *(complex<double> &x) {
Func<T> r;
int i;
for (i = 0; i < N; i++) r.a[i] = a[i] * x;
return r;
}
Func<T> operator *=(complex<double> &x) {
int i;
for (i = 0; i < N; i++) a[i] *= x;
return *this;
}
Func<T> operator *=(Func<T> &x) {
int i;
for (i = 0; i < N; i++) a[i] *= x[i];
return *this;
}
void FFT(bool f) {
complex<double> x, y;
T z;
double w;
int i, j, k;
for (i = 0; i < N; i++) {
j = Rev(i);
if (i < j) {
z = a[i];
a[i] = a[j];
a[j] = z;
}
}
for (i = 1; i < N; i <<= 1) {
w = PI / i;
if (f) w = -w;
x = complex<double>(cos(w), sin(w));
for (j = 0; j < N; j += i << 1) {
y = complex<double>(1);
for (k = 0; k < i; k++) {
z = a[i + j + k] * y;
a[i + j + k] = a[j + k] - z;
a[j + k] = a[j + k] + z;
y *= x;
}
}
}
if (f) {
x = complex<double>((double)1 / N);
for (i = 0; i < N; i++) a[i] *= x;
}
}
};
Func<Func<complex<double> > > X;
int main() {
double x;
int i, j, n, m;
scanf("%d%d", &n, &m);
for (i = 0; i <= n; i++) for (j = 0; j <= m; j++) {
scanf("%lf", &x);
X[i][j] = complex<double>(x, 0);
}
for (j = 0; j < N; j++) X[j].FFT(false);
X.FFT(false);
X *= X;
X.FFT(true);
for (j = 0; j < N; j++) X[j].FFT(true);
for (i = 0; i <= n + n; i++) {
for (j = 0; j <= m + m; j++) printf("%.0f ", X[i][j].real());
printf("\n");
}
}
~
XOR Convolution
\(2^{M}\) 크기의 배열 \(A_{i}\)와 \(B_{i}\)에 대하여 다음처럼 정의되는 배열 \(C_{i}\)를 구하는 문제를 생각해 보자. (\(\oplus\)는 Bitwise XOR 연산을 뜻함)
\[C_{i}=\sum_{k=0}^{2^{M}-1}A_{k}B_{i\oplus k}\]
XOR 연산은 각 자리별로 더한 뒤 2로 나눈 나머지를 구하는 것과 같으므로, \(\mathbb{Z}_{2^{M}}\rightarrow\mathbb{C}\)를 다음처럼 펼쳐서 생각해 보자.
\[\mathbb{Z}_{2}\rightarrow\mathbb{Z}_{2}\rightarrow\cdots\rightarrow\mathbb{Z}_{2}\rightarrow\mathbb{C}\]
각각이 이산 푸리에 변환 가능하므로, \(\mathbb{Z}_{2}\)에 대한 이산 푸리에 변환을 \(M\)번 행하는 것임을 알 수 있다.
#include <cstdio>
#include <complex>
const int SZ = 20, N = 1 << SZ;
using namespace std;
int Rev(int x) {
int i, r = 0;
for (i = 0; i < SZ; i++) {
r = r << 1 | x & 1;
x >>= 1;
}
return r;
}
void FFT(int *a, bool f) {
int i, j, k, z;
for (i = 0; i < N; i++) {
j = Rev(i);
if (i < j) {
z = a[i];
a[i] = a[j];
a[j] = z;
}
}
for (i = 1; i < N; i <<= 1) for (j = 0; j < N; j += i << 1) for (k = 0; k < i; k++) {
z = a[i + j + k];
a[i + j + k] = a[j + k] - z;
a[j + k] += z;
}
if (f) for (i = 0; i < N; i++) a[i] /= N;
}
int X[N];
int main() {
int i, n;
scanf("%d", &n);
for (i = 0; i < 1 << n; i++) scanf("%d", &X[i]);
FFT(X, false);
for (i = 0; i < N; i++) X[i] *= X[i];
FFT(X, true);
for (i = 0; i < 1 << n; i++) printf("%d ", X[i]);
}
