ps:고속_푸리에_변환
목차
고속 푸리에 변환 (Fast Fourier Transform, FFT)
개념
- 푸리에 변환은 신호처리등의 용도로 각종 공학분야에서 너무도 많이 쓰이는 것이지만, PS와 관련돼서는 합성곱의 빠른 계산 그리고 다항식의 빠른 곱셈에 사용된다.
- 그래서, 관련 내용을 공부하기 위해 자료를 찾을 때에는 PS쪽의 컨텍스트로 작성된 글인지 확인해서 봐야 한다.. 시간에 대한 함수를 주파수 성분으로 변환한다는 개념으로 설명하는 자료로 공부를 시작하면, 그렇게 이해한 내용을 PS쪽의 컨텍스트까지 연결시키는 데까지 시간이 너무 오래 걸린다.
- PS에서의 응용에 포커스를 맞춰서 이산 푸리에 변환을 설명하는 자료들도, '합성곱'의 계산으로 설명하는 방식과 '다항식 곱셈'으로 설명하는 방식이 있다. 물론 결국은 그 두 개가 같은 거긴 한데 (수열의 합성곱이 생성함수의 곱이다), 처음에 개념부터 간단히 잡아보려는 단계에서는 섞어서 보면 헷갈릴 여지가 꽤 있다.
- 경험상으로는 '다항식 곱셈'의 관점에서 푸리에 변환을 coefficient representation을 point-value representation으로 바꾸는 과정이라는 관점으로 설명하는 방식이 처음에 이해하기는 더 쉬웠다. 그리고 이쪽이 신호처리쪽에서 쓰는 정통 설명에서도 출발했을 때에도 이쪽으로의 연결이 좀 더 쉬운 느낌.
- 가장 이해하기 쉬운 자료는 이것
- 하지만 이쪽으로 시작하고 나면 '합성곱'의 개념으로 연결해 나가기가 조금 어려웠고, 구현과도 잘 연결이 안됐다.
- 반면에 '합성곱'의 개념으로 접근하는 자료인 The Casterian - 고속 푸리에 변환 은, 처음에는 조금 덜 와닫는 편이다. 하지만 처음부터 인풋을 함수가 아닌, 주기를 갖고 순환하는 수열로 정의하고, 합성곱 역시 거기에 맞춰서 계산되는 cyclic convolution 로 정의해놓고 설명하는 부분이, 실제 구현과도 잘 맞아들어간다. 그래서 '다항식 곱셈'과 연결하는 과정에서도 차이가 무엇인지가 좀더 명확하게 들어오게 되었다.
- 그 외에 두 개념을 섞어서 설명하는 곳들 (예) 의 경우는 오히려 이해하기가 더 어려웠다; 뭔가 명확하게 차이를 짚어주는 것이 도움이 될만한 부분을 얼렁뚱땅 넘어가는 느낌인데.. 이부분이 설명이 더 필요했다는 것도 나중에 다 알고보니 보이는 것이지, 처음에는 뭔가 연결이 안되는 느낌은 있는데 정확히 거기가 어딘지는 모르겠고 매우 답답했다;
- 헷갈렸던 포인트를 정리해보자..
- a = [a1,a2,a3], b = [b1,b2,b3] 의 두 수열 또는 이것을 계수로 하는 두 다항식이 있을 때, 이것의 합성곱을, 두 다항식을 곱한 다항식의 계수 [a1*b1, a1*b2+a2*b1, a1*b3+a2*b2+a3*b1, a2*b3+a3*b2, a3*b3] 로 설명하는 경우가 있다. 합성곱이 인풋보다 길이가 길어지는데, 그냥 그런가보다 하고 넘어갔다가 그 뒤에서 계속 혼란을 겪었다. (사실 수식을 제대로 쳐다보고 갔으면 헷갈릴 일도 없는데, 그냥 말로 써진 설명만 적당히 읽으며 넘어갔더니 이런 일이..ㅜ)
- 우선 명확하게 해야 할 것은, fft와 ifft는 수열의 길이를 바꾸지 않는다. 길이 n짜리 배열에 fft/ifft를 돌리면 결과 배열도 길이 n이다. 물론 정의역은 다르다
- a와 b가 주기를 갖는 무한 수열이라고 생각하면 합성곱도 주기를 갖고 [a1*b1+a2*b3+a3*a2, a1*b2+a2*b1+a3*b3, a1*b3+a2*b2+a3*b1] 가 나타나는 무한 수열이 된다. 이렇게 생각하면 결과도 인풋과 길이가 같아진다. 이것이 실제 fft알고리즘을 돌려서 나오는 결과이다.
- 그러면 이렇게 사이클릭 컨볼루션을 계산해주는 FFT 방법으로, 다항식 곱셈에서 필요로 하는 [a1*b1, a1*b2+a2*b1, a1*b3+a2*b2+a3*b1, a2*b3+a3*b2, a3*b3] 과 같은 그냥 컨볼루션 결과를 얻으려면, 수열을 [a1,a2,a3,0,0], [b1,b2,b3,0,0] 와 같이 변경해서 돌려야 한다.
- 이렇게 인풋에 0을 추가해서 변경하는 과정이 굉장히 킬포임에도, 설명이나 코드를 볼때 강조가 잘 안되고 너무도 자연스럽게 넘어가는 것은, 우리가 쓰는 FFT 알고리즘인 Cooley–Tukey FFT algorithm은 수열을 2^n 크기로 변환을 해서 돌리는 것을 요구하기 때문이다. 그래서 2^n 크기로 변환시키려고 어차피 0을 붙이다 보니, 저 과정이 티가 안난다. (실제로 [a1,a2,a3,0,0,0,0,0], [b1,b2,b3,0,0,0,0,0]의 형태로 변환해서 돌린다.) 오히려 2^n에 맞춰 리사이징을 할때, 사이클릭 컨볼루션을 계속 정확히 얻을 수 있도록 하려면, 단순히 0을 붙이는 것이 아니라 수열을 한번 더 붙여줘야 한다 ([a1,a2,a3,0,0,0,0,0], [b1,b2,b3,0,0,b1,b2,b3]와 같이 변환해야 한다). 이렇게 보면 사이클릭 컨볼루션이 더 손이 많이 가는 작업같은 느낌은 들지만, 입력수열이 애초부터 2^n 크기인 경우에는 사이클릭 컨볼루션을 구하기 위해서는 리사이징 없이 그냥 넘겨도 되지만, 다항식곱셈을 하려면 여전히 0을 더 붙여줘야 한다.
- 관련 자료들
응용
- 위에서 열심히 설명했 듯, FFT를 써서 효율적으로 풀 수 있는 문제는 기본적으로는 합성곱과 다항식의 곱셈이다.
- 사실 이것만으로는 용도가 제한적인것 같은데, 합성곱과 다항식의 곱셈을 염두에 두고 다른 문제들을 바라보면, 전혀 관계 없는 듯한 문제가 합성곱이나 다항식의 곱셈 문제로 변환되는 경우가 있다.
큰 수의 곱셈
- 다항식의 곱셈을 응용하면 큰 수의 곱셈을 효율적으로 계산할 수 있다
- 어떤 큰 수를, n가 특정 값일때 어떤 다항식의 값으로 표현할 수 있다. 예를 들어 4321이라는 수는, f(x)=4*x^3 + 3*x^2 + 2*x + 1 이란 다항식에 대해서 f(10) = 4321 이 된다. 그러면 두개의 큰 수 n과 m을 f(10)과 g(10)으로 표현하고 h(x)=f(x)*g(x)을 fft로 계산하면, n*m = h(10)이 된다.
- 위키에 의하면, 빠른 곱셈을 하는 알고리즘들로 카라츠바, 톰쿡, 그리고 FFT 기반으로 Schönhage–Strassen과 Fürer's algorithm이 소개되는데, 여기에서 FFT 기반 알고리즘의 시간복잡도는 O(nlogn)보다 더 큰 것으로 기록되어있다. 대체 무슨 차이인지 궁금한데 나중에 찾아보자..
All possible sums
- 두 자연수의 배열 a[]와 b[]가 있을때 a[i]+b[j]로 나올수 있는 모든 가능한 값과, 그 값을 만드는 (i,j)의 갯수를 빠르게 구할 수 있다.
- 대충 a+b=x 꼴을 주고서 만족하는 a,b의 갯수를 물으면 일단 FFT로 풀리는지 생각해보자. BOJ에도 이것을 살짝만 변형한 문제들이 꽤 많다. P_i+P_j=x, P_i+2*P_j=x, a^2+b^2=c^2, 등등
- FFT를 쓰니까 O(NlogN) 인데 N이 a,b의 크기가 아니라, a,b에서 가장 큰 원소의 값이라는 것에 유의.
- 예를 들어 설명하면, a=[1, 2, 3] and b=[2, 4] 라고 하면, f(x)=x^1 + x^2 + x^3, g(x) = x^2 + x^4 이라고 정의한다. 그리고 f(x)*g(x)를 계산하면 x^3 + x^4 + 2*x^5 + x^6 + x^7 인데, 이것은 set(a[i]+b[j]) = {3,4,5,6,7} 이고, a[i]+b[j]==5인 (i,j)는 2개라는 뜻이다.
- 처음 들었을 때에는 아이디어에 굉장히 감탄했다. 이게 FFT를 공부하고 바로 이어서 공부하게 되다보니, FFT로 다항식 곱셈을 처리하는 발상 - 계수를 fft로 변환하고 그것들을 곱한뒤 ifft해서 처리한다는 발상 - 이 워낙 기상천외해서, 이 발상의 임팩트가 상대적으로 약하게 느껴질 수도 있는데.. 이것도 내게는 굉장히 기발했다.
다항식의 빠른 나눗셈
- [To be filled]
- 이것은 다시 Multi-point Evaluation 이나, 키타마사법 등에 사용된다.
Multi-point Evaluation
구현
- FFT는 FT를 fast하게 (O(nlogn)에) 푸는 방법을 총칭하는 것이지 특정 알고리즘의 이름은 아니다. 내가 본 모든 자료에서는 FFT의 구현을 위해서 Cooley–Tukey FFT algorithm을 사용했다. 그 외에 Good-Thomas, Bruun's, Rader's, Bluestein's FFT 등이 있긴 하다고 하는데 본적은 없다.
- 사실, 이 알고리즘은 Cooley와 Tukey가 1960년대에 발표해서 이런 이름이 붙었지만, 이미 160년 전에 가우스가 이 알고리즘을 발견했다는 것이 이후에 뒤늦게 알려졌다.
- 기본적으로는 재귀형태이지만, 재귀를 풀어서 쓰려면 약간의 테크닉이 필요하다.
- Python3에서는 complex(math.cos(ang), math.sin(ang)) 대신 cmath.exp(complex(0,ang)) 으로 더 간단하게 쓰는 방법도 있다. 속도는 비슷한것 같다.
- 비슷하게, complex(1)을 1+0j 로 쓰는 것도 가능하다. 이 편이 좀 더 빠른듯 했다.
- 그리고, 사실상 PS에서는 FFT만을 단독으로 사용하는 일은 없다. 거의 모든 경우에서, FFT를 하고 계수를 곱하고 IFFT를 하는 것이 하나의 작업으로 묶인다. 그것을 고려하면 조금 더 효율적인 코드 작성도 가능하다.
- 그러나, 이런 최적화는 사소한 차이일뿐, Python에서 FFT를 직접 구현해서 쓰는 것은 굉장히 느리다. BOJ에 있는 몇가지 대표적인 FFT문제에 대해서 속도를 재보면 c++과 너무도 심하게 차이가 난다. Pypy를 동원해도 아슬아슬하게 컷에 걸리는 정도.
코드
def fft(a: Iterable[complex], inverse: Optional[bool] = False) -> List[complex]:
"""Returns the (inverse) discrete Fourier transform of the given sequence.
The length of the sequence should be power of 2.
"""
a = list(a)
n = len(a)
if n & (n - 1):
raise ValueError(f'len(a) ({len(a)}) is not power of 2.')
j = 0
for i in range(1, n):
bit = n >> 1
while j >= bit:
j -= bit
bit >>= 1
j += bit
if i < j:
a[i], a[j] = a[j], a[i]
p = (-2 if inverse else 2) * math.pi
s = 2
while s <= n:
ws = cmath.exp(complex(0, p / s))
for i in range(0, n, s):
w = 1 + 0j
for j in range(i, i + (s >> 1)):
tmp = a[j + (s >> 1)] * w
a[j + (s >> 1)] = a[j] - tmp
a[j] += tmp
w *= ws
s <<= 1
return [x / n for x in a] if inverse else a- 그러나 실제로 이 코드는 사용하지 않는다. FFT를 사용하지 않는 FFT 참고.
FFT를 사용하지 않는 FFT
- 제목만으로는 무슨 장난같은 소리인가 싶은데.. FFT를 안 쓰고도 FFT를 필요로 하는 문제를 풀기 위한 꼼수이다.
- 정석적인 FFT 구현이 파이썬에서 너무 느리다는 것을 고민하다 보니 이런 변태적인 아이디어가 떠올랐다. 내가 검색한 결과로는 이전에 이런 짓을 했던 사람을 발견하지 못했다..
- PS에서 FFT의 활용이 (합성곱 = 다항식 곱셈 ⇒ 큰 수의 곱셈) 으로 사용된다는 점을 생각하자. 그런데 Python에서는 이미 내부적으로 정수형이 큰수를 지원하고, 큰수의 곱셈이 빠르게 처리되도록 구현되어 있다. Python의 큰 정수 표현법 2 - 사칙연산 을 참고하면, 내부적으로는 카라츠바 알고리즘을 사용해서 곱셈을 수행하도록 구현되어 있는 듯 하다. 아무튼 그게 FFT가 아니라 카라츠바 알고리즘이더라도 내부적으로는 C 모듈로 구현되어 있으므로 순수 파이썬 코드보다는 훨씬 빠를 것을 기대할 수 있다. 그렇다면 역으로 이러한 “빠른 큰수의 곱셈”을 이용해서 “빠른 다항식 곱셈”을 계산할 수도 있지 않을까하는 생각에 도달하게 된다.
- 간단하게 생각해서 계수가 [3,4,5]인 다항식과 [4,5,6]인 다항식을 곱하는 것을. 300040005와 '400050006라는 두개의 큰 수를 만들어서 이걸 곱하는 것이다. 그러면 120031005800490030 라는 결과가 나오는데, 이것을 끊어서 나오는 [12,31,58,49,30]이 실제 두 다항식을 곱한 결과의 계수인 것.
- 이걸 좀 더 최적화 하기 위해 테스트 해본 바로는, 십진 표현보다는 이진 표현이 빠르고, 이진수를 만들때 비트연산으로 큰 수를 만드는 것보다.. 문자열 연산으로 이진수를 만들고, 숫자로 변환하는 것이 더 빨랐다. (사실 좀 어이없다;)
- 곱셈 결과에서 계수들을 정확히 추출하기 위해서는 carry가 생기지 않아야 하기 때문에, 처음에 큰 수를 만들때 각 숫자 앞에 0을 붙이게 되는데, 당연히 이 0의 갯수가 작아야 속도가 빠르다. 0의 갯수는 결과로 나올수 있는 계수의 최댓값을 고려해서 붙여야 하는데, min(len(a), len(b)) * max(a) * max(b) 뭐 이런식으로 대충 상한을 잡을수는 있을 것이다. 그러나 귀찮으니 그냥 64를 디폴트로 줬다.. 문제에서 다른 언급이 없다면 결과값이 64비트 정수형 안데 들어가긴 할테니..
- 사실 조금 궁리해보면 carry가 있는 계산 값에서도 계수들을 추출하는 방법이 있을것 같기도 하다. 대충 0의 갯수를 두 종류로 줘서 결과값을 각각 계산하고 그것을 비교하면 복원이 가능하지 않을까 싶긴 한데.. 생각하기 귀찮다.
- 짜보면서도, fft에 비해 n값이 커지고, 알고리즘 시간복잡도도 O(nlogn) 대신 O(n^log3)이 되고, 문자열 변환하는 오버로드도 생기는 등등의 요인으로 이게 실제로 빨라질지는 의구심이 있었다. 그러나 이 방법을 사용하자, Golf Bot 문제의 런타임이 기존 5876ms가 1996ms로 줄어들었다
- 그 외의 특징은
- 이 방법을 쓰면 배열의 길이를 2^n 형태로 늘릴 필요가 없다. 그것도 속도 향상의 요인.
- 정수만 사용하므로 실수 오차도 없다.
- 코드도 짧아진다.
Teflib FFT
- 이렇게 기존보다 빠른 구현을 완성하고, 관련 문제들을 풀다가 이상한 결과를 발견했다. 큰 수 곱셈 문제는 파이썬의 경우는 그냥 정수 곱셈을 써서 풀면 풀리고 그게 가장 빨라야 할 문제인데, 대부분의 답이 5000ms 정도인 가운데 몇개만 100ms 정도의 솔루션이 있었다. 차이는 계산을 decimal.Decimal 로 변환해서 한 것뿐. 너무나 이상한 결과..
- 검색해도 여기에 대한 내용을 쉽게 못찾고 있다가, 공식 매뉴얼에서 근거를 찾았다.
Q. CPython 구현은 커다란 수에서 빠릅니까?
A. 예. CPython 과 PyPy3 구현에서, decimal 모듈의 C/CFFI 버전은 임의의 정밀도로 올바르게 자리 올림 되는 십진 부동 소수점 산술을 위한 고속 libmpdec 라이브러리를 통합합니다 1. libmpdec는 중간 크기의 숫자에는 카라추바 곱셈(Karatsuba multiplication)을 사용하고 매우 큰 숫자에는 수론적 변환(Number Theoretic Transform)을 사용합니다.
- 그래서, 레알 큰 수의 곱셈을 할때는 그냥 기본 정수형을 갖고 하는 것보다, decimal 로 바꿔서 하는게 더 빠른것 같다.
- 관련 자료가 잘 안나오는 것은, 진짜로 레알 큰 수의 곱셈을 할 사람들은 이미 numpy나 다른 라이브러리를 쓰고 있어서가 아닐까..
- 암튼.. 그래서 나도 Decimal을 쓰도록 구현을 수정해서 multiply v2.0 을 만들었다.
- 이제 진짜로 FFT를 하긴 한다! 그래서 제목을 'FFT를 사용하지 않는 FFT' 대신에 다른 이름을 붙여야 할것 같은데 뭐라고 불러야 할지가 어렵다.. 그냥 구현체가 teflib에 있으므로, 텦립FFT 라고 일단 부르기로 했다;
- 속도는 다시 한번 향상되었다
- 이동 (N≤60,000) : 1736ms ⇒ 844ms ⇒ 444ms (cpp: 28ms)
- Golf Bot (N≤200,000) : 5876ms ⇒ 1384ms ⇒ 724ms (cpp: 88ms)
- 르모앙의 추측 (N≤1,000,000) : 시간초과 ⇒ 시간초과 ⇒ 2900ms (PyPy:7436ms ⇒ 6244ms ⇒ 시간초과, cpp: 252ms)
- 드디어 PyPy를 사용하지 않고도 통과하는 데에 성공. 사실 오히려 PyPy에서는 시간초과가 나기 시작했다. (원래 대부분의 실행시간이 c-extension인 경우에는 pypy가 더 느리다)
코드
한계점
- 속도도 빠르고, 코드도 짧고 장점이 많긴 한데, 정석이 아니다보니 한계도 있다.
- 모듈러를 취한 값을 계산하는 문제
- 정석 방식으로도 그냥은 안되고 NTT를 써야 하는 문제이기는 하다. 내 구현으로는 안되는게 맞긴 하지만 굳이 해보자면 그냥 모듈러를 안취한 값을 다 구하고서, 마지막에 그 값에 모듈러를 하는 방법이 있기는 하다. 숫자가 말도 안되게 커질것 같아보이지만, a*b%m 를 64bit 정수형으로 계산할 수 있게 하려면 m은 2^32 이내로 주어져야 한다. 그러면 처음 주어지는 값들은 모듈러를 먼저 하고 시작할 수 있으니까 대충 10^9 이내. FFT의 계산 속도를 감안해서 n의 크기는 보통 많아야 10^6 이내로 주어진다. 그러면 결과값에 들어올 수 있는 가장 큰 값은 10^9 * 10^9 * 10^6 해서 10^24 정도이다. 수 한개당 24자리를 쓰고 그게 백만개니까 2400만 자리의 숫자를 곱하는 건데.. 씽크스몰 문제에서 이미 수 한개당 18자리를 써서 1800만 자리의 숫자의 곱셈을 5~6초 정도에 통과한 적이 있다. 문제 기본 제한시간을 2초정도로만 준다면 해볼만하지 않을까..
- 인풋에 음수가 포함된 경우
- 내 구현에서 다항식을 큰수로 변환한것은, 인풋에 0 또는 양의 정수만 있는 경우를 가정했기에 가능했고, 음수가 있으면 이 방법으로는 안된다. 사실 다행히도 대부분의 문제에서는 음수가 들어오는 경우가 별로 없긴 하다. 그러나 몇몇 경우, 예를 들면 다항식 나눗셈을 할 때, 에는 음수를 포함한 수열에 FFT를 돌려야 한다.
- 처음 떠오른 방법은 위에서 한것처럼 모듈러 처리를 하는 것이다. -x 대신에 C-x 를 인풋으로 넣고 계산한다음에, 나온 결과를 C로 나눈 나머지를 구해서 그 값이 C/2보다 크면 다시 거기에서 C를 뺀 값을 쓰고, 아니면 그 값을 그냥 쓰는 것이다. 가능하긴 한데.. 이렇게 하려면 C/2가 정상적으로 나올수 있는 최댓값보다 커야 한다. 그러면.. 입력에 들어오는 수가 절대값 10^a 이하이고, 갯수가 10^b라고 하면 최댓값은 10^(2a+b).. 그래서 C를 이 값으로 잡게 된다면, 변환된 입력인 C-x의 크기가 10^(2a+b)니까, 최댓값은 10^(4a+3b). 그러면 수 하나를 표기하기 위해 써야 하는 자리수가 2~3배 더 늘어난다. 음.. 못할거는 아니긴 하다만;
- 다르게 떠오른 방법은 인풋 수열을 양수와 음수로 쪼개서 계산하는 것. 다항식의 곱셈이라고 생각하면 f(x)를 모든 계수가 양수인 f1(x)와 모든 계수가 음수인 f2(x)의 합으로 표현하는 것. 그러면 f(x)*g(x) = f1(x)*g1(x) + f1(x)*g2(x) + f2(x)*g1(x) + f2(x)*g2(x) 로 나눠서 계산이 가능하다. 곱셈을 4번 하고 더해주면 된다. 식은 좀 더 깔끔한거 같은데. 속도는 위의 방법보다 더 느릴것 같기도 하고..
- 인풋에 실수가 포함된 경우
- PS 범위에서 설마 이렇게 지저분한 문제가 나올까..ㅜㅜ 나온다면 뭐 답에서 요구하는 유효숫자만큼 곱해서 정수로 바꿔놓고 계산하면 불가능하지는 않을것 같긴 하다.
NTT
ps/고속_푸리에_변환.txt · 마지막으로 수정됨: 2023/08/04 15:09 저자 teferi
토론