고속 푸리에 변환¶
이 글에서는 길이가 $n$인 두 다항식을 $O(n \log n)$ 시간에 곱할 수 있는 알고리즘을 다룹니다. 이는 $O(n^2)$ 시간이 걸리는 단순한 곱셈보다 빠릅니다. 두 개의 큰 수를 곱하는 것도 당연히 다항식의 곱셈으로 환원할 수 있으므로, 두 개의 큰 수 역시 $O(n \log n)$ 시간에 곱할 수 있습니다(여기서 $n$은 수의 자릿수입니다).
고속 푸리에 변환(FFT)의 발견은 1965년에 알고리즘을 발표한 Cooley와 Tukey의 업적으로 여겨집니다. 그러나 사실 FFT는 그 이전에도 여러 번 발견되었으나, 현대적인 컴퓨터가 발명되기 전까지는 그 중요성이 이해되지 않았습니다. 일부 연구자들은 FFT의 발견을 1924년 Runge와 König의 업적으로 봅니다. 하지만 실제로는 Gauss가 이미 1805년에 이러한 방법을 개발하였으나, 이를 발표하지는 않았습니다.
여기서 소개하는 FFT 알고리즘은 $O(n \log n)$ 시간에 동작하지만, 임의로 큰 계수를 가진 임의로 큰 다항식의 곱셈이나 임의로 큰 정수의 곱셈에는 사용할 수 없다는 점에 유의하십시오. 이 알고리즘은 계수가 작은 크기 $10^5$의 다항식이나 크기 $10^6$인 두 수의 곱셈을 쉽게 처리할 수 있으며, 이는 보통 경쟁 프로그래밍 문제를 푸는 데 충분합니다. $10^6$비트 크기의 수를 곱하는 규모를 넘어서면 계산 중에 사용하는 부동소수점 수의 범위와 정밀도가 정확한 최종 결과를 내기에 충분하지 않습니다. 다만 임의로 큰 다항식/정수 곱셈을 수행할 수 있는 더 복잡한 변형들이 존재합니다. 예를 들어 1971년에 Schönhage와 Strasser는 환(ring) 구조에서 FFT를 재귀적으로 적용하여 임의로 큰 수를 곱하는 변형을 개발하였으며, 이는 $O(n \log n \log \log n)$에 동작합니다. 그리고 최근(2019년)에는 Harvey와 van der Hoeven이 진정한 $O(n \log n)$에 동작하는 알고리즘을 발표하였습니다.
이산 푸리에 변환¶
차수가 $n - 1$인 다항식이 있다고 합시다:
일반성을 잃지 않고, 계수의 개수인 $n$이 $2$의 거듭제곱이라고 가정합니다. $n$이 $2$의 거듭제곱이 아니라면, 빠진 항 $a_i x^i$를 추가하고 계수 $a_i$를 $0$으로 두면 됩니다.
복소수 이론에 따르면 방정식 $x^n = 1$은 $n$개의 복소수 해(이를 $n$차 단위근이라고 부릅니다)를 가지며, 그 해는 $k = 0 \dots n-1$에 대해 $w_{n, k} = e^{\frac{2 k \pi i}{n}}$의 형태입니다. 또한 이 복소수들은 매우 흥미로운 성질을 몇 가지 가지고 있습니다. 예를 들어 주 $n$차 단위근 $w_n = w_{n, 1} = e^{\frac{2 \pi i}{n}}$을 사용하여 다른 모든 $n$차 단위근을 나타낼 수 있습니다: $w_{n, k} = (w_n)^k$.
다항식 $A(x)$(또는 이와 동등하게 계수 벡터 $(a_0, a_1, \dots, a_{n-1})$)의 이산 푸리에 변환(DFT)은 점 $x = w_{n, k}$에서의 다항식 값으로 정의됩니다. 즉, 다음 벡터입니다:
마찬가지로 역 이산 푸리에 변환을 정의합니다: 다항식의 값 $(y_0, y_1, \dots, y_{n-1})$의 역 DFT는 다항식의 계수 $(a_0, a_1, \dots, a_{n-1})$입니다.
따라서 정방향 DFT가 $n$차 단위근인 점들에서 다항식의 값을 계산한다면, 역 DFT는 그 값들을 이용하여 다항식의 계수를 복원할 수 있습니다.
DFT의 응용: 다항식의 빠른 곱셈¶
두 다항식 $A$와 $B$가 있다고 합시다. 각각에 대해 DFT를 계산합니다: $\text{DFT}(A)$와 $\text{DFT}(B)$.
이 다항식들을 곱하면 어떻게 됩니까? 당연히 각 점에서 값들이 단순히 곱해집니다. 즉,
이는 한 벡터의 각 원소에 다른 벡터의 대응하는 원소를 곱하는 방식으로 벡터 $\text{DFT}(A)$와 $\text{DFT}(B)$를 곱하면, 다름 아닌 다항식의 DFT인 $\text{DFT}(A \cdot B)$를 얻는다는 뜻입니다:
마지막으로 역 DFT를 적용하면 다음을 얻습니다:
우변에서 두 DFT의 곱은 벡터 원소끼리의 쌍별 곱을 의미합니다. 이는 $O(n)$ 시간에 계산할 수 있습니다. DFT와 역 DFT를 $O(n \log n)$에 계산할 수 있다면, 두 다항식의 곱(그리고 결과적으로 두 개의 큰 수의 곱)도 같은 시간 복잡도로 계산할 수 있습니다.
두 다항식의 차수가 같아야 한다는 점에 유의해야 합니다. 그렇지 않으면 DFT의 두 결과 벡터의 길이가 달라집니다. 값이 $0$인 계수를 추가하여 차수를 맞출 수 있습니다.
또한 두 다항식을 곱한 결과는 차수가 $2 (n - 1)$인 다항식이므로, 각 다항식의 차수를 두 배로 늘려야 합니다(이 역시 $0$을 채워 넣어 해결합니다). 값이 $n$개인 벡터로부터는 계수가 $2n - 1$개인 원하는 다항식을 복원할 수 없습니다.
고속 푸리에 변환¶
고속 푸리에 변환은 DFT를 $O(n \log n)$ 시간에 계산할 수 있게 해 주는 방법입니다. FFT의 기본 아이디어는 분할 정복을 적용하는 것입니다. 다항식의 계수 벡터를 두 벡터로 나누고, 각각에 대해 DFT를 재귀적으로 계산한 다음, 그 결과를 결합하여 전체 다항식의 DFT를 계산합니다.
그러면 차수가 $n - 1$인 다항식 $A(x)$가 있다고 합시다. 여기서 $n$은 $2$의 거듭제곱이고 $n > 1$입니다:
이 다항식을 더 작은 두 다항식으로 나눕니다. 하나는 짝수 위치의 계수만 포함하고, 다른 하나는 홀수 위치의 계수를 포함합니다:
다음이 성립함을 쉽게 알 수 있습니다
다항식 $A_0$과 $A_1$은 다항식 $A$의 절반만큼의 계수만 가집니다. $\text{DFT}(A_0)$과 $\text{DFT}(A_1)$을 이용해 $\text{DFT}(A)$를 선형 시간에 계산할 수 있다면, 시간 복잡도에 대해 점화식 $T_{\text{DFT}}(n) = 2 T_{\text{DFT}}\left(\frac{n}{2}\right) + O(n)$을 얻으며, 이는 마스터 정리에 의해 $T_{\text{DFT}}(n) = O(n \log n)$이 됩니다.
이를 어떻게 해낼 수 있는지 알아봅시다.
벡터 $\left(y_k^0\right)_{k=0}^{n/2-1} = \text{DFT}(A_0)$과 $\left(y_k^1\right)_{k=0}^{n/2-1} = \text{DFT}(A_1)$을 계산했다고 가정합시다. $\left(y_k\right)_{k=0}^{n-1} = \text{DFT}(A)$에 대한 식을 찾아봅시다.
처음 $\frac{n}{2}$개의 값에 대해서는 앞서 언급한 방정식 $A(x) = A_0(x^2) + x A_1(x^2)$를 그대로 사용할 수 있습니다:
그러나 나머지 $\frac{n}{2}$개의 값에 대해서는 약간 다른 식을 찾아야 합니다:
여기서 다시 $A(x) = A_0(x^2) + x A_1(x^2)$과 두 항등식 $w_n^n = 1$, $w_n^{n/2} = -1$을 사용했습니다.
따라서 전체 벡터 $(y_k)$를 계산하는 원하는 공식을 얻습니다:
(이 $a + b$와 $a - b$ 패턴을 때때로 나비라고 부릅니다.)
이렇게 해서 DFT를 $O(n \log n)$ 시간에 계산하는 방법을 배웠습니다.
역 FFT¶
차수가 $n - 1$인 다항식 $A$의 점 $x = w_n^k$에서의 값들인 벡터 $(y_0, y_1, \dots y_{n-1})$이 주어졌다고 합시다. 이 다항식의 계수 $(a_0, a_1, \dots, a_{n-1})$을 복원하려고 합니다. 이 잘 알려진 문제를 보간이라고 하며, 이를 푸는 일반적인 알고리즘들이 있습니다. 하지만 이 특수한 경우에는(단위근에서의 점의 값을 알고 있으므로) 훨씬 더 간단한 알고리즘(실질적으로 정방향 FFT와 동일한 알고리즘)을 얻을 수 있습니다.
DFT는 그 정의에 따라 행렬 형태로 쓸 수 있습니다:
이 행렬을 방데르몽드 행렬이라고 합니다.
따라서 벡터 $(y_0, y_1, \dots y_{n-1})$의 왼쪽에 이 행렬의 역행렬을 곱하여 벡터 $(a_0, a_1, \dots, a_{n-1})$을 계산할 수 있습니다:
간단히 확인해 보면 이 행렬의 역행렬이 다음과 같은 형태임을 알 수 있습니다:
따라서 다음 공식을 얻습니다:
이를 $y_k$에 대한 공식과 비교하면
두 문제가 거의 같다는 것을 알 수 있습니다. 따라서 계수 $a_k$도 순방향 FFT와 똑같은 분할 정복 알고리즘으로 구할 수 있습니다. 다만 $w_n^k$ 대신 $w_n^{-k}$를 사용해야 하고, 마지막에 결과 계수들을 $n$으로 나누어야 합니다.
따라서 역 DFT의 계산은 순방향 DFT의 계산과 거의 같으며, 역시 $O(n \log n)$ 시간에 수행할 수 있습니다.
구현¶
여기서는 간단한 재귀 방식의 FFT 구현과 역 FFT 구현을 제시합니다. 순방향 FFT와 역 FFT의 차이가 매우 작으므로 둘을 하나의 함수에 담았습니다. 복소수를 저장하기 위해 C++ STL의 complex 타입을 사용합니다.
using cd = complex<double>;
const double PI = acos(-1);
void fft(vector<cd> & a, bool invert) {
int n = a.size();
if (n == 1)
return;
vector<cd> a0(n / 2), a1(n / 2);
for (int i = 0; 2 * i < n; i++) {
a0[i] = a[2*i];
a1[i] = a[2*i+1];
}
fft(a0, invert);
fft(a1, invert);
double ang = 2 * PI / n * (invert ? -1 : 1);
cd w(1), wn(cos(ang), sin(ang));
for (int i = 0; 2 * i < n; i++) {
a[i] = a0[i] + w * a1[i];
a[i + n/2] = a0[i] - w * a1[i];
if (invert) {
a[i] /= 2;
a[i + n/2] /= 2;
}
w *= wn;
}
}
이 함수는 계수 벡터를 전달받아 DFT 또는 역 DFT를 계산하고, 그 결과를 다시 이 벡터에 저장합니다. 인자 $\text{invert}$는 순방향 DFT와 역 DFT 중 어느 것을 계산할지를 나타냅니다. 함수 내부에서는 먼저 벡터의 길이가 1인지 확인하며, 그렇다면 아무것도 할 필요가 없습니다. 그렇지 않으면 벡터 $a$를 두 벡터 $a0$과 $a1$로 나누고, 두 벡터 각각에 대해 재귀적으로 DFT를 계산합니다. 그다음 값 $wn$과, $wn$의 현재 거듭제곱을 담을 변수 $w$를 초기화합니다. 그런 다음 위의 공식들을 사용하여 결과 DFT의 값들을 계산합니다.
플래그 $\text{invert}$가 설정되어 있으면 $wn$을 $wn^{-1}$로 바꾸고, 결과의 각 값을 $2$로 나눕니다(이 작업은 재귀의 각 단계마다 수행되므로, 결국 최종 값들을 $n$으로 나누게 됩니다).
이 함수를 사용하여 두 다항식을 곱하는 함수를 만들 수 있습니다:
vector<int> multiply(vector<int> const& a, vector<int> const& b) {
vector<cd> fa(a.begin(), a.end()), fb(b.begin(), b.end());
int n = 1;
while (n < a.size() + b.size())
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<int> result(n);
for (int i = 0; i < n; i++)
result[i] = round(fa[i].real());
return result;
}
이 함수는 정수 계수를 갖는 다항식에 대해 동작하지만, 다른 타입에 대해서도 동작하도록 수정할 수 있습니다. 복소수로 연산할 때에는 약간의 오차가 생기므로, 마지막에 결과 계수들을 반올림해야 합니다.
마지막으로 두 긴 수를 곱하는 함수는 다항식을 곱하는 함수와 사실상 다르지 않습니다. 그 후에 해야 할 유일한 일은 수를 정규화하는 것입니다:
int carry = 0;
for (int i = 0; i < n; i++)
result[i] += carry;
carry = result[i] / 10;
result[i] %= 10;
}
두 수의 곱의 길이는 두 수의 길이의 합을 절대 넘지 않으므로, 벡터의 크기는 모든 자리올림 연산을 수행하기에 충분합니다.
개선된 구현: 제자리(in-place) 계산¶
효율을 높이기 위해 재귀 구현에서 반복 구현으로 전환하겠습니다. 위의 재귀 구현에서는 벡터 $a$를 명시적으로 두 벡터로 분리했습니다. 짝수 위치의 원소들은 한 임시 벡터에, 홀수 위치의 원소들은 다른 임시 벡터에 할당했습니다. 그러나 원소들을 특정한 방식으로 재배열하면 이러한 임시 벡터를 만들 필요가 없습니다(즉, 모든 계산을 벡터 $A$ 자체에서 바로 "제자리(in-place)"로 수행할 수 있습니다).
첫 번째 재귀 단계에서는 위치의 최하위 비트가 0인 원소들이 벡터 $a_0$에 할당되고, 위치의 최하위 비트가 1인 원소들이 $a_1$에 할당되었다는 점에 주목하십시오. 두 번째 재귀 단계에서도 같은 일이 일어나지만, 이번에는 두 번째 하위 비트를 기준으로 하며, 이후 단계도 마찬가지입니다. 따라서 각 계수의 위치의 비트를 뒤집고 이 뒤집힌 값을 기준으로 정렬하면 원하는 순서를 얻습니다(이를 비트 반전 순열(bit-reversal permutation)이라고 합니다).
예를 들어 $n = 8$일 때 원하는 순서는 다음과 같은 형태입니다:
실제로 첫 번째 재귀 단계(중괄호로 둘러싸인 부분)에서 벡터는 $[a_0, a_2, a_4, a_6]$과 $[a_1, a_3, a_5, a_7]$의 두 부분으로 나뉩니다. 보다시피 비트 반전 순열에서 이는 단순히 벡터를 두 절반, 즉 처음 $\frac{n}{2}$개의 원소와 마지막 $\frac{n}{2}$개의 원소로 나누는 것에 해당합니다. 그다음 각 절반에 대해 재귀 호출이 이루어집니다. 각 절반에 대한 결과 DFT가 원소들 자리에 그대로 반환된다고 합시다(즉, 각각 벡터 $a$의 앞쪽 절반과 뒤쪽 절반에 반환됩니다.
이제 두 DFT를 합쳐 전체 벡터에 대한 하나의 DFT를 만들고자 합니다. 원소들의 순서가 이상적이므로 두 DFT를 합치는 과정도 이 벡터 안에서 직접 수행할 수 있습니다. 원소 $y_0^0$과 $y_0^1$을 가져와 나비 변환을 수행할 수 있습니다. 결과로 나온 두 값의 위치는 처음 두 값의 위치와 같으므로 다음을 얻습니다:
마찬가지로 $y_1^0$과 $y_1^1$의 나비 변환을 계산하여 결과를 제자리에 넣을 수 있으며, 나머지도 같은 방식으로 처리합니다. 그 결과 다음을 얻습니다:
이렇게 하여 벡터 $a$로부터 필요한 DFT를 계산했습니다.
여기서는 첫 번째 재귀 단계에서만 DFT를 계산하는 과정을 설명했지만, 나머지 모든 단계에서도 똑같이 적용된다는 것은 분명합니다. 따라서 비트 반전 순열을 적용한 뒤에는 추가 메모리 없이 제자리(in-place)에서 DFT를 계산할 수 있습니다.
또한 이를 통해 재귀를 없앨 수도 있습니다. 가장 낮은 단계에서 시작하기만 하면 됩니다. 즉, 벡터를 쌍으로 나누고 각 쌍에 나비 변환을 적용합니다. 그러면 마지막 단계의 작업이 적용된 벡터 $a$를 얻습니다. 다음 단계에서는 벡터를 크기 $4$의 벡터들로 나누고 다시 나비 변환을 적용하며, 이를 통해 크기 $4$인 각 블록의 DFT를 얻습니다. 이런 식으로 계속됩니다. 마지막 단계에서는 $a$의 두 절반에 대한 DFT 결과를 얻은 상태이며, 나비 변환을 적용하여 전체 벡터 $a$의 DFT를 얻습니다.
using cd = complex<double>;
const double PI = acos(-1);
int reverse(int num, int lg_n) {
int res = 0;
for (int i = 0; i < lg_n; i++) {
if (num & (1 << i))
res |= 1 << (lg_n - 1 - i);
}
return res;
}
void fft(vector<cd> & a, bool invert) {
int n = a.size();
int lg_n = 0;
while ((1 << lg_n) < n)
lg_n++;
for (int i = 0; i < n; i++) {
if (i < reverse(i, lg_n))
swap(a[i], a[reverse(i, lg_n)]);
}
for (int len = 2; len <= n; len <<= 1) {
double ang = 2 * PI / len * (invert ? -1 : 1);
cd wlen(cos(ang), sin(ang));
for (int i = 0; i < n; i += len) {
cd w(1);
for (int j = 0; j < len / 2; j++) {
cd 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 (cd & x : a)
x /= n;
}
}
먼저 각 원소를 반전된 위치의 원소와 교환하여 비트 반전 순열을 적용합니다. 그런 다음 알고리즘의 $\log n - 1$개 단계에서 해당 크기 $\text{len}$인 각 블록의 DFT를 계산합니다. 이 블록들은 모두 같은 단위근 $\text{wlen}$을 가집니다. 모든 블록을 순회하며 각 블록에 나비 변환을 수행합니다.
비트 반전을 더 최적화할 수 있습니다. 이전 구현에서는 인덱스의 모든 비트를 순회하여 비트 단위로 반전된 인덱스를 만들었습니다. 그러나 비트를 다른 방식으로 반전할 수도 있습니다.
$j$에 이미 $i$를 반전한 값이 들어 있다고 가정합시다. 그러면 $i + 1$로 넘어가려면 $i$를 1 증가시켜야 하고, $j$도 1 증가시켜야 하는데, 이때 $j$는 "반전된" 수 체계에서 증가시켜야 합니다. 일반적인 이진법에서 1을 더하는 것은 끝에 연속된 모든 1을 0으로 뒤집고 그 바로 앞의 0을 1로 뒤집는 것과 같습니다. 마찬가지로 "반전된" 수 체계에서는 앞쪽에 연속된 모든 1을 뒤집고, 그다음의 0도 뒤집습니다.
따라서 다음과 같은 구현을 얻습니다:
using cd = complex<double>;
const double PI = acos(-1);
void fft(vector<cd> & a, bool invert) {
int n = a.size();
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 * PI / len * (invert ? -1 : 1);
cd wlen(cos(ang), sin(ang));
for (int i = 0; i < n; i += len) {
cd w(1);
for (int j = 0; j < len / 2; j++) {
cd 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 (cd & x : a)
x /= n;
}
}
또한 비트 반전 순열을 미리 계산해 둘 수 있습니다. 이는 모든 호출에서 크기 $n$이 같을 때 특히 유용합니다. 하지만 호출이 세 번뿐인 경우(두 다항식을 곱하는 데 필요한 횟수)에도 효과가 눈에 띕니다. 또한 모든 단위근과 그 거듭제곱을 미리 계산해 둘 수 있습니다.
수론적 변환¶
이제 목표를 조금 바꿔 보겠습니다. 여전히 두 다항식을 $O(n \log n)$ 시간에 곱하려 하지만, 이번에는 계수를 어떤 소수 $p$로 나눈 나머지로 계산하려 합니다. 물론 이 작업에는 일반적인 DFT를 사용한 뒤 결과에 모듈로 연산자를 적용할 수 있습니다. 그러나 그렇게 하면 특히 큰 수를 다룰 때 반올림 오차가 발생할 수 있습니다. 수론적 변환(number theoretic transform, NTT)은 정수만으로 동작하므로 결과가 정확함이 보장된다는 장점이 있습니다.
이산 푸리에 변환은 복소수와 1의 $n$제곱근에 기반합니다. 이를 효율적으로 계산하기 위해 근의 성질(예를 들어 거듭제곱으로 다른 모든 근을 생성하는 근이 하나 존재한다는 성질)을 폭넓게 활용합니다.
그런데 모듈러 산술에서의 1의 $n$제곱근에 대해서도 같은 성질이 성립합니다. 원시체에서의 $n$차 단위근은 다음을 만족하는 수 $w_n$입니다:
나머지 $n-1$개의 근은 근 $w_n$의 거듭제곱으로 얻을 수 있습니다.
이를 고속 푸리에 변환 알고리즘에 적용하려면 $2$의 거듭제곱인 어떤 $n$에 대해 근이 존재해야 하고, 그보다 작은 모든 거듭제곱에 대해서도 근이 존재해야 합니다. 다음과 같은 흥미로운 성질에 주목할 수 있습니다:
따라서 $w_n$이 $n$차 단위근이면, $w_n^2$은 $\frac{n}{2}$차 단위근입니다. 그러므로 2의 더 작은 모든 거듭제곱에 대해서도 필요한 차수의 근이 존재하며, 이 근들은 $w_n$을 사용해 계산할 수 있습니다.
역 DFT를 계산하려면 $w_n$의 역원 $w_n^{-1}$이 필요합니다. 그런데 법이 소수이면 역원은 항상 존재합니다.
따라서 $n$차 단위근이 존재하는 충분히 큰 법 $p$가 있다면, 복소수 근에서 필요한 모든 성질을 모듈러 연산에서도 똑같이 사용할 수 있습니다.
예를 들어 다음 값을 사용할 수 있습니다: 법 $p = 7340033$, $w_{2^{20}} = 5$. 이 법으로 충분하지 않다면 다른 쌍을 찾아야 합니다. $p = c 2^k + 1$ 꼴의 법($p$는 소수)에 대해서는 $2^k$차 단위근이 항상 존재한다는 사실을 이용할 수 있습니다. $g$가 $p$의 원시근일 때, $g^c$가 그러한 $2^k$차 단위근임을 보일 수 있습니다.
const int mod = 7340033;
const int root = 5;
const int root_1 = 4404020;
const int root_pw = 1 << 20;
void fft(vector<int> & a, bool invert) {
int n = a.size();
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) {
int wlen = invert ? root_1 : root;
for (int i = len; i < root_pw; i <<= 1)
wlen = (int)(1LL * wlen * wlen % mod);
for (int i = 0; i < n; i += len) {
int w = 1;
for (int j = 0; j < len / 2; j++) {
int u = a[i+j], v = (int)(1LL * a[i+j+len/2] * w % mod);
a[i+j] = u + v < mod ? u + v : u + v - mod;
a[i+j+len/2] = u - v >= 0 ? u - v : u - v + mod;
w = (int)(1LL * w * wlen % mod);
}
}
}
if (invert) {
int n_1 = inverse(n, mod);
for (int & x : a)
x = (int)(1LL * x * n_1 % mod);
}
}
여기서 함수 inverse는 모듈러 역원을 계산합니다(모듈러 곱셈 역원 참고). 상수 mod, root, root_pw는 법과 근을 결정하며, root_1은 mod를 법으로 하는 root의 역원입니다.
실제로 이 구현은 (매우 많은 나머지 연산 때문에) 복소수를 사용하는 구현보다 느리지만, 메모리 사용량이 적고 반올림 오차가 없다는 등의 장점이 있습니다.
임의의 법에 대한 곱셈¶
여기서는 이전 절과 같은 목표를 달성하려 합니다. 두 다항식 $A(x)$와 $B(x)$를 곱하고, 계수를 어떤 수 $M$으로 나눈 나머지로 계산하는 것입니다. 수론적 변환은 특정한 소수에 대해서만 동작합니다. 법이 원하는 형태가 아닌 경우에는 어떻게 해야 합니까?
한 가지 방법은 $c 2^k + 1$ 꼴의 서로 다른 소수들로 수론적 변환을 여러 번 수행한 다음, 중국인의 나머지 정리를 적용해 최종 계수를 계산하는 것입니다.
또 다른 방법은 다항식 $A(x)$와 $B(x)$를 각각 더 작은 두 다항식으로 나누는 것입니다
여기서 $C \approx \sqrt{M}$입니다.
그러면 $A(x)$와 $B(x)$의 곱은 다음과 같이 나타낼 수 있습니다.
다항식 $A_1(x)$, $A_2(x)$, $B_1(x)$, $B_2(x)$는 $\sqrt{M}$보다 작은 계수만을 가지므로, 등장하는 모든 곱의 계수는 $M \cdot n$보다 작으며, 이는 보통 일반적인 부동소수점 자료형으로 다루기에 충분히 작습니다.
따라서 이 방법은 (일반적인 FFT와 역 FFT를 사용하여) 더 작은 계수를 가진 다항식들의 곱을 계산해야 하며, 그 후 모듈러 덧셈과 곱셈을 사용하여 $O(n)$ 시간에 원래의 곱을 복원할 수 있습니다.
응용¶
DFT는 언뜻 보기에는 다항식 곱셈과 전혀 관계없어 보이는 매우 다양한 다른 문제에도 사용할 수 있습니다.
가능한 모든 합¶
두 배열 $a[]$와 $b[]$가 주어집니다. 가능한 모든 합 $a[i] + b[j]$를 구하고, 각 합이 몇 번 나타나는지 세어야 합니다.
예를 들어 $a = [1,~ 2,~ 3]$이고 $b = [2,~ 4]$이면 다음을 얻습니다. 합 $3$은 $1$가지 방법으로, 합 $4$도 $1$가지 방법으로, $5$는 $2$가지, $6$은 $1$가지, $7$은 $1$가지 방법으로 얻을 수 있습니다.
배열 $a$와 $b$에 대해 두 다항식 $A$와 $B$를 구성합니다. 배열의 수는 다항식의 지수 역할을 하고($a[i] \Rightarrow x^{a[i]}$), 이 항의 계수는 그 수가 배열에 나타나는 횟수가 됩니다.
그런 다음 이 두 다항식을 $O(n \log n)$ 시간에 곱하면 다항식 $C$를 얻는데, 여기서 지수는 어떤 합을 얻을 수 있는지를, 계수는 그 합을 몇 가지 방법으로 얻을 수 있는지를 알려 줍니다. 이를 예제에 적용하여 보이면 다음과 같습니다.
가능한 모든 스칼라곱¶
길이가 $n$인 두 배열 $a[]$와 $b[]$가 주어집니다. $a$와 $b$의 각 순환 이동 사이의 곱을 모두 계산해야 합니다.
크기가 $2n$인 새로운 배열 두 개를 만듭니다. $a$를 뒤집고 그 뒤에 $n$개의 0을 덧붙입니다. 그리고 $b$는 단순히 자기 자신 뒤에 덧붙입니다. 이 두 배열을 다항식으로 보고 곱한 뒤, 곱 $c$의 계수 $c[n-1],~ c[n],~ \dots,~ c[2n-2]$를 살펴보면 다음을 얻습니다.
그리고 $i \ge n$인 모든 원소에 대해 $a[i] = 0$이므로 다음이 성립합니다.
이 합이 벡터 $a$와 $b$를 왼쪽으로 $(k - (n - 1))$번 순환 이동한 것의 스칼라곱에 불과하다는 것을 쉽게 알 수 있습니다. 따라서 이 계수들이 문제의 답이며, 여전히 $O(n \log n)$ 시간에 이를 구할 수 있었습니다. 여기서 $c[2n-1]$도 $n$번째 순환 이동을 제공하지만, 이는 $0$번째 순환 이동과 같으므로 답에서 따로 고려할 필요가 없다는 점에 유의하십시오.
두 개의 띠¶
두 불리언 띠($0$과 $1$ 값으로 이루어진 순환 배열) $a$와 $b$가 주어집니다. 어떤 위치에서도 첫 번째 띠의 $1$이 두 번째 띠의 $1$과 맞닿지 않도록 첫 번째 띠를 두 번째 띠에 붙이는 모든 방법을 찾고자 합니다.
이 문제는 사실 앞의 문제와 크게 다르지 않습니다. 두 띠를 붙인다는 것은 단지 두 번째 배열에 순환 이동을 수행한다는 뜻이며, 두 배열의 스칼라곱이 $0$이면 두 띠를 붙일 수 있습니다.
문자열 매칭¶
소문자로 이루어진 두 문자열, 텍스트 $T$와 패턴 $P$가 주어집니다. 텍스트에서 패턴이 나타나는 모든 위치를 계산해야 합니다.
각 문자열에 대해 다항식을 만듭니다($T[i]$와 $P[I]$는 알파벳 $26$개 글자에 대응하는 $0$에서 $25$ 사이의 수입니다).
여기서
그리고
여기서
식 $P[m-i-1]$을 사용하면 패턴이 명시적으로 뒤집힌다는 점에 주목하십시오.
두 다항식의 곱 $C(x) = A(x) \cdot B(x)$의 $(m-1+i)$번째 계수는 패턴이 텍스트의 위치 $i$에 나타나는지를 알려 줍니다.
여기서 $\alpha_{i+j} = \frac{2 \pi T[i+j]}{26}$이고 $\beta_j = \frac{2 \pi P[j]}{26}$입니다
일치하는 경우 $T[i+j] = P[j]$이므로 $\alpha_{i+j} = \beta_j$입니다. (피타고라스 삼각 항등식을 사용하면) 다음을 얻습니다.
일치하지 않는 경우 적어도 한 문자가 다르므로 곱 $a_{i+1} \cdot b_{m-1-j}$ 중 하나는 $1$과 같지 않으며, 따라서 계수 $c_{m-1+i} \ne m$이 됩니다.
와일드카드를 포함한 문자열 매칭¶
이는 앞의 문제를 확장한 것입니다. 이번에는 패턴에 모든 문자와 일치할 수 있는 와일드카드 문자 $\*$가 포함되는 것을 허용합니다. 예를 들어 패턴 $a*c$는 텍스트 $abccaacc$에서 정확히 세 위치, 즉 인덱스 $0$, 인덱스 $4$, 인덱스 $5$에 나타납니다.
$P[m-i-1] = *$이면 $b_i = 0$으로 설정한다는 점을 제외하면 완전히 같은 다항식을 만듭니다. $x$가 $P$에 있는 와일드카드의 개수라면, $c_{m-1+i} = m - x$일 때 $T$의 인덱스 $i$에서 $P$가 일치합니다.
연습 문제¶
- SPOJ - POLYMUL
- SPOJ - MAXMATCH
- SPOJ - ADAMATCH
- Codeforces - Yet Another String Matching Problem
- Codeforces - Lightsabers (hard)
- Codeforces - Running Competition
- Kattis - A+B Problem
- Kattis - K-Inversions
- Codeforces - Dasha and cyclic table
- CodeChef - Expected Number of Customers
- CodeChef - Power Sum
- Codeforces - Centroid Probabilities