FFT

푸리에변환을 바라보는 관점

1. Domain을 시간에서 주파수로 변경하기. 이것의 장점은 연속/불연속에 상관없이 통합적으로 이해할 수 있으며, 확률론에서의 생성함수가 어떤 의미인지도 여기(도메인의 변경)에 가까운것 같다. 푸리에변환을 검색했을때 나오는 거의 모든 자료가 이러한 관점에서 설명을 한다.

2. 주기가 n인 수열을 생성함수 F(x) = sigma{k>=0}{a_k * x^k}로 나타내고 거기에 x=exp(2pi/n*i*k*x)로 치환한것. NOTE: exp(2pi/n*i)의 거듭제곱은 주기가 n이다. 이 관점의 장점은 받아들이기 쉽다는것이다. 시간이니 주파수니 뜬금없는 개념대신, 수열을 생성함수로 표현하는건 자연스럽고, 그 생성함수에서 주기성을 적용하는것도 자연스럽게 납득할 수 있을것이다. 나무위키왈 라플라스도 라플라스변환이전에 생성함수를 먼저 생각해냈다고. NTT도 주기가 n인 변수를 exp대신 mod의 원시근을 사용한것으로 자연스럽게 받아들여진다. (뇌피셜) 어려운 fft관련 문제들은 이러한 관점에서 생성함수적으로 수식을 변형해 사용할것 같다. https://en.wikipedia.org/wiki/Generating_function#Relation_to_discrete-time_Fourier_transform 2번(생성함수)관점으로 설명하면 대부분을 생성함수에 연동시키고 필요한 내용만 깔끔하게 나올것 같은데, 본문은 아쉽게도 1번(도메인의 변경)관점만 알고있을때 작성되었다.

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

DFT를 빠르게 계산하는 알고리즘. 정의 그대로를 loop로 구하는 Naive한 알고리즘은 $O(N^2)$의 시간이 걸리는데, 이를 $O(NlogN)$에 구할 수 있는 다양한 알고리즘을 통틀어서 FFT라고 부른다.

푸리에 변환의 용도

합성곱 계산의 시간복잡도가 낮아진다. 합성곱이란것 자체가 이미지 프로세싱, 신호처리, 미분방정식 등 다양한 분야에서 활용되는 중요한 연산이다. ps에서의 구체적인 사용예시는 (big-integer/polynomial)의 (곱셈/나눗셈)정도?

(이산)(circular)합성곱

수열 $a, b$에 대해
$c_i=\sum\limits_{0\leq{}j<n}a_jb_{i-j}$로 정의되는 새로운 수열 $c$를 만드는 연산. 이 때 $b$에 음수 인덱스가 생기는데, $b$를 circular하게 확장하여 사용한다. 즉 양수인덱스 $i$에 대해 $b_{-i}=b_{n-i}$이다. 합성곱을 뜻하는 표기법은 다음과 같다 $c=a*b$
나이브한 반복문으로 구하면 $O(n^2)$이다. 이 합성곱을 $O(n^2)$보다 빠르게 구하는게 사실상 불가능해 보이지만 푸리에변환을 통해 빠르게 구할 수 있다는 것이 핵심이다.
푸리에변환은 다음과 같은 성질이 있다.
$F^{-1}(F(a)F(b))=a*b$
즉, 수열 $a, b$를 FFT로 $O(n\ log\ n)$에 푸리에변환하고, 단순한 $O(n)$의 함수곱셈 후 FFT로 $O(n\ log\ n)$에 역변환하면 총 시간복잡도 $O(n\ log\ n)$에 $a*b$를 구할 수 있다.

Cooley-Tukey FFT

길이가 2의 거듭제곱꼴인 수열 $a$에 대해 그 푸리에변환은 다음과 같다.
$F(a)_k=\sum\limits_{0\leq{}j<n}e^{-2\pi{}ijk/n}a_j$
$=\sum\limits_{0\leq{}j<n/2}e^{-2\pi{}i(2j)k/n}a_{2j}+\sum\limits_{0\leq{}j<n/2}e^{-2\pi{}i(2j+1)k/n}a_{2j+1}$
$=\sum\limits_{0\leq{}j<n/2}e^{-2\pi{}ijk/(n/2)}a_{2j}+e^{-2\pi{}ik/n}\sum\limits_{0\leq{}j<n/2}e^{-2\pi{}ijk/(n/2)}a_{2j+1}$
여기서 $a$의 짝수번째 수들로 만든 수열을 $a.even$, 홀수번째 수들로 만든 수열을 $a.odd$라고 하자. 각 수열은 길이가 $n/2$임을 상기하며 위 식을 다음과 같이 쓸 수 있다.
$F(a)_k=F(a.even)_k+e^{-2\pi{}ik/n}F(a.odd)_k$
단계마다 n이 2배 줄어들고, 각 단계마다 $a.even$과 $a.odd$를 나누는데 $O(n)$의 시간이 걸린다. 마스터정리에 의해 $O(n\ log\ n)$이다.

푸리에 역변환 알고리즘

푸리에 역변환은 다음과 같다. 증명은 생략한다.
$a_k=\frac{1}{n}\sum\limits_{0\leq{}j<n}e^{2\pi{}ijk/n}F(a)_j$
푸리에 변환과 거의 비슷한데 $e$의 지수에 마이너스가 없고 $\frac{1}{n}$이 붙었다. FFT를 활용할 수 있도록 살짝 변형해보자. 우선 $e^{2\pi{}ijk/n}$과 $F(a)_j$는 주기가 $n$인 $j$의 함수들이므로 그 곱 또한 주기가 $n$인 함수가 된다. 즉, 아무곳에서나 연속된 $n$개의 항을 더해도 값이 같다. 그러므로 다음과 같이 $n$부터 왼쪽 $n$개항의 합으로 바꿔보자.
$=\frac{1}{n}\sum\limits_{0\leq{}j<n}e^{2\pi{}i(n-j)k/n}F(a)_{n-j}$
$=\frac{1}{n}\sum\limits_{0\leq{}j<n}e^{-2\pi{}ijk/n}F(a)_{n-j}$
이 식이 의미하는 바는 다음과 같다. $F(a)$ 배열의 $[1,n)$구간 원소들을 뒤집고 FFT한 다음 $n$으로 나눠주면 푸리에 역변환이 된다!

임의의 n으로 확장하는 법

수열 $a$에 $n*2$ 보다 큰 최소의 2의 거듭제곱수 크기만큼 0을 채워준다.
수열 $b$는 뒤집어서 circular하게 확장시켜주고 다시 뒤집는다(음수인덱스에 대해 circular하게 확장시켜주는 의미).
이렇게 확장한 수열 $a, b$에 대해 $c$를 구하면 $c$의 끝부분 $n$개의 수가 원래 수열의 circular $a*b$가 된다. 직접 손으로 해보면 왜 되는지 알 수 있으므로 증명은 생략한다.

다항식 곱셈 구하는법

다항식 곱셈은 다음과 같이 쓸 수 있다.
수열 $a, b$에 대해
$c_i=\sum\limits_{0\leq{}j\leq{}i}a_jb_{i-j}$
자세히보면 $\sum{}$의 범위가 합성곱과 다르다. 그렇지만 다음과 같이 구할 수 있다. 수열 $a, b$를 n*2 보다 큰 최소의 2의 거듭제곱수 크기만큼 0을 채워서 구한 c의 앞부분 $n$개가 원래 수열의 다항식곱셈이 된다. 같은 이유로 증명은 생략한다.

확장

(+,*)에서 알 수 있듯이, Ring에만 적용할 수 있다. min같은건 여기서 나가리. (min,+) Convolution이란게 실제로 있는거같은데(https://codeforces.com/blog/entry/98663), FFT가 아니라 볼록성으로 가지고 뭐 하는 일종의 슬로프트릭이라는듯? mod 2에서는 (xor, and)이려나? 아직 이해못했는데 Walsh Hadamard Transform 이게 비슷한것 같다. xor컨볼루션이란게 사실 그냥 비트별로 컨볼루션 독립적으로 돌리면 되는거 아닌가 안풀어봐서 잘 모르겠음. 자료1 자료2 관련문제1

(+,xor)컨볼루션은 xor을 곱셈2개로 쪼개서 (+,*)컨볼루션 2개로 쪼개 계산하는게 가능하다는거같다. 어느정도 웰논인건지 어째선지 Atcoder Beginner에도 나왔다(...) XOR을 곱의합으로 쪼개는 트릭이 종종 보이는거같다. 외워두는게 좋을듯 ABC196F

구현

Python3
C++
구현 특이사항으로, $e^{-2\pi{}ik/n}$을 매번 계산하는게 상당한 시간이 걸린다. $e^{-2\pi{}i/n}$를 계산하고 제곱하여 $e^{-2\pi{}ik/n}$를 구하는건 복소수 곱셈만으로 가능하므로 훨씬 빠르다. $e^{-2\pi{}i/n}$ 또한 처음 한번만 계산하고, 이후에는 재귀인자로 제곱해서 넘겨주면 훨씬 빠르다.

순환합성곱으로 다항식곱셈을 구할 수 있듯이, 반대로 다항식곱셈으로 순환합성곱 구하는걸 생각해봤다. 같은 아이디어로 a를 circular expansion하고 b를 0 expansion하여 곱한 결과의 특정부분이 circular convolution이 된다. 그리고 다항식곱셈은 bigint곱셈에서 올림을 막기 위한 적당한 크기의 padding bit들을 추가해주면 구할 수 있다. 따라서 python이나 java등의 bigint로 합성곱을 간단하게 구할 수 있다!
199바이트 합성곱
가독성 향상버전

연습문제

입력배열B를 뒤집으면 $A*B$의 최댓값을 구하는 문제가 된다.

available[x] && available[y] => available[x+y]이므로 큰수곱셈이 된다.

.

.

.

여담

특정 소수모듈러 수체계에 대해, 정수연산만으로 FFT를 구할 수 있는 NTT(Number Theoretic Transform?)라는것도 있다. 푸리에변환이 지수함수e의 주기성을 사용하는걸 NTT에서는 소수P의 원시근(계속 제곱하면 1~P-1이 모두 나오는 수?)으로 대체하여 계산하는 알고리즘이라고 한다. 현재까지 NTT를 강요하는 문제는 만난적이 없기 때문에 굳이 NTT를 공부하지는 않았다.
추가메모리O(1)과 재귀없이 fft를 하는 최적화가 가능하지만 시간복잡도상의 변화는 없으니 굳이 다루지는 않겠다. 다만 상수차이가 꽤 커서 팀노트에는 최적화버젼이 같이 들어가면 좋을것이다. NTT또한 FFT를 이해하였다면 팀노트에 적어서 사용하는데 문제는 없을 것이다.

참고자료

https://en.wikipedia.org/wiki/Convolution
https://en.wikipedia.org/wiki/Fast_Fourier_transform
https://ko.wikipedia.org/wiki/%ED%91%B8%EB%A6%AC%EC%97%90_%EB%B3%80%ED%99%98
https://topology-blog.tistory.com/6
https://koosaga.com/139
https://cubelover.tistory.com/19
https://cubelover.tistory.com/12