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

고속 푸리에 변환은 이산 푸리에 변환(Discrete Fourier Transform, DFT)과 그 역변환을 빠르게 수행하는 효율적인 알고리즘으로 알려져 있다. 이에 대한 이론적인 설명을 자세하게 하려는 목적은 아니다. FFT가 PS에 어떤 식으로 응용될 수 있는지에 대해서 아주 간단하게 소개하고자 한다.

 

PS 문제들을 풀다보면 FFT를 사용해야하는 경우를 대면하게 되는데, 이러한 문제들은 백날 생각해도 모르면 못 풀고, 조금을 생각하더라도 알면 푼다고 생각하기 때문에 꼭 알아둬야 되는 부분이라고 생각한다. 우선, PS에서 FFT가 사용되는 대표적인 예는 convolution를 $O(n \log_2 n)$만에 계산할 때이다.

 

Convolution이란, 여기에서 설명된 것 처럼 크기가 $n$인 배열 $a$, $b$가 있다고 했을 때, 아래와 같은 식으로 계산된 $c$를 convolution이라고 부른다. 이 배열들은 0-indexed 이다.

 

$c_i = \sum\limits^{i}_{j=0}{a_j \times b_{i-j}}$

 

일반적으로 $c$를 계산하는 시간복잡도는 $O(n^2)$이지만, FFT를 이용하면 계산하는 시간복잡도는 $O(n \log_2 n)$로 매우 빠르다. FFT에 대한 자세한 설명은 생략하고, $a$, $b$가 주어졌을 때 $c$를 계산하는 코드와 응용 문제들을 설명하겠다.

 

이러한 방법을 이해하기 위해서는 이론에 대해서 깊이 있는 공부를 해야하기 때문에 이해하지 않고 팀노트에 등록하는 것도 좋아보인다.

 

#define _USE_MATH_DEFINES #include <math.h> #include <complex> #include <vector> #include <algorithm> using namespace std; #define sz(v) ((int)(v).size()) #define all(v) (v).begin(), (v).end() typedef complex<double> base; void fft(vector<base> &a, bool invert) { int n = sz(a); for (int i=1,j=0;i<n;i++){ int bit = n >> 1; for (;j>=bit;bit>>=1) j -= bit; j += bit; if (i < j) swap(a[i],a[j]); } for (int len=2;len<=n;len<<=1){ double ang = 2*M_PI/len*(invert?-1:1); base wlen(cos(ang), sin(ang)); for (int i=0;i<n;i+=len){ base w(1); for (int j=0;j<len/2;j++){ base u = a[i+j], v = a[i+j+len/2]*w; a[i+j] = u+v; a[i+j+len/2] = u-v; w *= wlen; } } } if (invert){ for (int i=0;i<n;i++) a[i] /= n; } } template <typename T> vector<T> multiply(const vector<T>& a, const vector<T>& b) { vector<base> fa(all(a)), fb(all(b)); int n = 1, m = sz(a)+sz(b)-1; while (n < m) n <<= 1; fa.resize(n); fb.resize(n); fft(fa, false); fft(fb, false); for (int i=0;i<n;i++) fa[i] *= fb[i]; fft(fa, true); vector<T> ret(m); for (int i=0;i<m;i++) ret[i] = static_cast<T>(fa[i].real()+(fa[i].real()>0?0.5:-0.5)); return ret; }

#include <bits/stdc++.h> using namespace std; #define sz(v) ((int)(v).size()) #define all(v) (v).begin(), (v).end() using complex_t = complex<double>; void fft(vector<complex_t>& a) { int n = sz(a), L = 31 - __builtin_clz(n); static vector<complex<long double>> R(2, 1); static vector<complex_t> rt(2, 1); // (^ 10% faster if double) for (static int k = 2; k < n; k *= 2) { R.resize(n); rt.resize(n); auto x = polar(1.0L, acos(-1.0L) / k); for (int i=k;i<k+k;i++) rt[i] = R[i] = i&1 ? R[i/2] * x : R[i/2]; } vector<int> rev(n); for (int i=0;i<n;i++) rev[i] = (rev[i / 2] | (i & 1) << L) / 2; for (int i=0;i<n;i++) if (i < rev[i]) swap(a[i], a[rev[i]]); for (int k = 1; k < n; k *= 2){ for (int i = 0; i < n; i += 2 * k) for (int j=0;j<k;j++) { // complex_t z = rt[j+k] * a[i+j+k]; // (25% faster if hand-rolled) /// include-line auto x = (double *)&rt[j+k], y = (double *)&a[i+j+k]; /// exclude-line complex_t z(x[0]*y[0] - x[1]*y[1], x[0]*y[1] + x[1]*y[0]); /// exclude-line a[i + j + k] = a[i + j] - z; a[i + j] += z; } } } template <typename T> vector<T> conv(const vector<T>& a, const vector<T>& b) { if (a.empty() || b.empty()) return {}; vector<T> res(sz(a) + sz(b) - 1); int L = 32 - __builtin_clz(sz(res)), n = 1 << L; vector<complex_t> in(n), out(n); copy(all(a), begin(in)); for (int i=0;i<sz(b);i++) in[i].imag(b[i]); fft(in); for (complex_t& x: in) x *= x; for (int i=0;i<n;i++) out[i] = in[-i & (n - 1)] - conj(in[i]); fft(out); for (int i=0;i<sz(res);i++){ res[i] = static_cast<T>(imag(out[i]) / (4 * n) + (is_integral_v<T> ? (imag(out[i]) > 0 ? 0.5 : -0.5) : 0)); } return res; }

 

응용 문제 1) 부분 배열 곱

 

문제) 크기가 $N$인 정수 배열 $A$와 크기가 $M$인 정수 배열 $B$가 있다. ($M \leq N \leq 500,000$). 크기가 $M$인 $A$의 (연속한) 부분 배열 $C$가 있을 때, 함수 $f$의 정의는 다음과 같다. $f(C) = \sum\limits^{M}_{i=0}{B[i] \times C[i]}$. 이 때, $f(C)$가 최대가 되는 $C$를 찾아 $f(C)$ 값을 구하시오.

 

응용 문제 2) HackerRank w6 - Best spot

 

응용 문제 3) BOJ 1067번 이동