FFT(Faster Fourier Transformation)에 대해서

https://codeforces.com/blog/entry/43499

https://codeforces.com/blog/entry/43499

^영어가 된다면 이것부터 보는 것을 과감히 추천^

 

<22.09.25 수정 : 오밤중에 써서 그런가 매끄럽지 않은 부분이 많아 수정함.>

 

 이번에 다수의 FFT를 이용한 문제 풀이를 하면서 이 알고리즘을 너무 어렵지 않은 선에서 간략히 정리해보고자 한다.

우선, "Fast"하지 않은, 그냥 기본적인 푸리에 변환을 어떤 경우에 이용할 수 있는지를 설명하는 것이 좋을 것 같다.

 

 위키피디아에 따르면 푸리에 변환은 "시간이나 공간에 대한 함수를 시간 또는 공간 주파수 성분으로 분해하는 변환"을 말한다. 이렇게 말하면 감이 잘 잡히지 않을 수 있는데, 어떤 함수를 그 함수를 정확히 나타낼 수 있을 만큼의 "점"들로 변환하는 것을 말한다.

 

 다항함수에 푸리에 변환을 하는 방법은 상당히 간단한데, 그냥 해당 함수에 아무 x값을 넣고 y값을 기록하는 일을 해당 다항함수의 차수 + 1만큼 해주면 되기 때문이다. 이유는 간단하다. 해당 함수의 계수가 n차 방정식에 대해서 n + 1개 만큼 있기 때문이다.

3차방정식을 예로 들면, $ a x^3 + b x^2 + c x + d $ 의 모든 계수를 알아내기 위해서는 서로 다른 방정식 위의 점 4개 (x0, y0),(x1, y1),(x2, y2),(x3, y3) 가 있으면 되고, 해당 점을 이용해서 4개의 식을 만들어내서 4원1차방정식을 만들어 풀면 원래의 함수를 복원할 수 있다.

 

 여기서 "대체 왜 이걸 점들로 바꾸냐?"고 생각이 들 수도 있다. 이를 알기 위해서는 여러 다항식을 푸리에 변환해야 한다. 2개의 다항식을 푸리에 변환한 후, 나온 점들을 각각 다항식 A에 대해서는

A(x) = {(x0, y0), (x1, y1), (x2, y2), ..., (xn - 1, yn - 1)}

이라고 하고, B에 대해서는

B(x) = {(a0, b0), (a1, b1), (a2, b2), ..., (an - 1, bn - 1)}

라고 하자. 여기서

$$ C = A * B $$

, 즉 C는 A와 B의 합성곱이라고 하자. C다항식의 푸리에 변환은 간단히 :

C(x) = {(a0 * x0, b0 * y0), (a1 * x1, b1 * y1), (a2 * x2, b2 * y2), ..., (an-1 * xn - 1, bn-1 * yn - 1)}

로 구할 수 있다! 이말은, 단순히 점들을 선형적을 곱함으로써 합성곱을 계산할 수 있다는 의미이다!

실제로 x값을 대충 자연수들로 지정하여 푸리에 변환을 하는 것을 DFT(Discrete Fourier Transformation)이라 부른다.

다만... 이것만으로는 $O(n^2)$의 한계를 넘기에는 부족하다. 크게 2가지 문제가 있다.

1. 점들을 대입해서 y값을 알아내는데 시간이 걸림

2. 합성곱이 된 함수를 복원하는 것도 시간이 걸림

이 2개의 작업이 푸리에변환을 이용한 합성곱의 시간복잡도가 $O(n)$보다 복잡하기 때문에 배보다 배꼽이 더 큰 상황이 발생한다...

 

이러한 단점을 보완한 것이 FFT이다.

"아무" 수나 넣는 대신, 우리는 1의 거듭제곱근(영어로 root of unity)를 이용한다. 

https://en.wikipedia.org/wiki/Root_of_unity

Root of unity - Wikipedia

위의 링크에 붙어있는 이미지는 보면, 1의 5거듭제곱근들을 볼 수 있다. 즉, 저 복소평면 위의 점들은 5제곱을 하면 전부 "1"이 된다! 물론 1도 포함이다 ㅎㅎ 그냥 원을 5등분 한거 아니냐고 할 수 있는데 실제로 굉장히 직관적이고 감사하게도 중심이 원점이고 반지름이 1인 원 위에 5개의 거듭제곱근들이 (1, 0i)을 원점에 대해 각각 108도 만큼 회전한 곳에 놓여 있다.

 

여기서 포인트는 그냥 1의 거듭제곱근을 선택하는게 아니라, 2^n승에 해당하는 1의 거듭제곱근을 선택하는 것이다.

이러면 2가지 이점이 있는데,

1. 그 거듭제곱근 중에서 1, -1, $i$,$-i$ 등이 포함되는 되어 연산이 간단해짐.

2. 홀수와 짝수 차수의 식들 끼지 모으면 같은 형태의 식이 계속 등장하기 때문에 Divide and Conquer를 이용해서 연산을 더욱 줄이면서 얻는 정보는 이전과 같은 상태를 유지할 수 있음.

이 그것이다.

<이게 사실 중요한 파트이고 어려운 파트이지만 여기까지 이해하기는 어렵기도 하고 필자도 100% 자신은 없기에 나중에 따로 글을 쓰던가, 아님 위의 블로그에 맡길 생각이다.>

 

아래는 해당 FFT를 c++로 구현한 것이다.

 

<FFT 노트(분할 정복을 재귀적으로 구현하지 않음으로써 성능을 높인 버전에 해당함.)>

#define _USE_MATH_DEFINES
#include <iostream>
#include <vector>
#include <cmath>
#include <algorithm>
#include <complex>

using namespace std;
typedef complex<double> base;

#define SZ(x) ((int)(x).size())
#define ALL(x) (x).begin(), (x).end()

void fft(vector<base> &a, bool invert){
    int n = SZ(a);
    for(int i = 1, j = 0; i < n; i++){
        int bit = n >> 1;
        while(j >= bit){
            j -= bit;
            bit >>= 1;
        }
        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;
        }
    }
}

void multiply(const vector<int> &a, const vector<int> &b, vector<int> &res){
    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);
    res.resize(n);
    for(int i = 0; i < n; i++){
        res[i] = (int)abs(round(fa[i].real()));
    }
}

 

FFT를 큰 수의 빠른 곱셈에도 이용하는 이유는 FFT는 "합성곱을 빠르게 하는 것에 특화"되어 있다고 한단계 추상화해서 생각하면 알 수 있다. 자연수 계수가 0이상 9이하의 정수이며 x = 10인 다항식으로 생각할 수 있기 때문이다.