¶ Biến đổi Fourier nhanh - Fast Fourier transform
Người viết:
- Nguyễn Hoàng Vũ - Trường Đại học Công nghệ, ĐHQGHN
Reviewer:
- Trần Xuân Bách - Đại học Chicago
- Nguyễn Minh Nhật - Trường THPT chuyên Khoa học Tự nhiên, ĐHQGHN
Các bài toán tổ hợp ngày càng xuất hiện nhiều trong các cuộc thi lập trình thi đấu và thường xuyên nắm giữ các vị trí khó nhất. Bài viết này sẽ giới thiệu về một công cụ quan trọng để giải các bài toán tổ hợp, đó là Biến đổi Fourier nhanh - Fast Fourier transform, hay còn được viết tắt là FFT.
¶ 1. Các kiến thức cần biết
¶ 1.1. Số phức (Complex number)
¶ Định nghĩa
Số phức là các số có dạng trong đó và . Số được gọi là đơn vị ảo. Tập hợp tất cả các số phức được kí hiệu là .
Ta có thể biểu diễn số phức trên mặt phẳng toạ độ bằng vectơ .
- Phần thực (Real) và phần ảo (Imaginary): .
- Số phức liên hợp (Conjugate): .
- Môđun (Modulus): .
- Acgumen (Argument): Góc có hướng từ trục thực đến vectơ , tức góc trong hình dưới đây.
- Dạng lượng giác (Polar form): .
Về biểu diễn , các bạn có thể tìm hiểu ở Công thức Euler (Wikipedia).
¶ Các phép toán
Xét hai số phức và :
- Phép cộng: .
- Phép trừ: .
- Phép nhân: .
- Phép chia: .
Phép cộng (trừ) hai số phức tương đương với phép cộng (trừ) hai vectơ biểu diễn chúng. Khi nhân hai số phức, ta nhân môđun của chúng và cộng acgumen của chúng.
¶ Căn đơn vị (Root of unity)
Xét một số nguyên dương :
- Căn đơn vị cấp là các số phức thoả mãn .
- Công thức tổng quát cho các căn đơn vị cấp : với từ đến . Kí hiệu thì ta còn có thể viết dưới dạng .
- Các căn đơn vị nằm trên đường tròn đơn vị tạo thành một đa giác đều đỉnh.
¶ Một số tính chất của căn đơn vị
- .
- Với chẵn, , từ đây suy ra .
- Với chẵn, .
¶ 1.2. Ma trận
Các bạn có thể tham khảo bài viết Nhân ma trận (VNOI).
¶ 2. Thuật toán FFT
¶ 2.1. Bài toán
Cho hai dãy và . Tính dãy được cho bởi công thức: với từ đến .
Dãy được định nghĩa như trên được gọi là tích chập (convolution) của hai dãy và .
¶ 2.2. Ý tưởng
Biểu thức của dãy phía trên gợi lại cho chúng ta về phép nhân đa thức. Cụ thể, định nghĩa:
thì chính là dãy hệ số của đa thức .
Có thể thấy rằng việc tính đa thức theo định nghĩa sẽ có độ phức tạp thời gian , không đủ nhanh khi độ dài khá lớn, ta cần một hướng tiếp cận khác. Giả sử ta tính được giá trị của và tại điểm khác nhau:
thì giá trị của tại các điểm tương ứng là: .
Ta có định lý quan trọng sau đây:
Định lý nội suy đa thức. Cho cặp số phức thoả mãn với mọi , tồn tại duy nhất một đa thức có bậc không quá thoả mãn với mọi .
Ví dụ:
- Có duy nhất một đường thẳng (đa thức bậc ) đi qua hai điểm bất kì trên mặt phẳng.
- Có duy nhất một parabol (đa thức bậc ) đi qua ba điểm bất kì trên mặt phẳng.
Ý tưởng của thuật toán FFT là chọn ra một tập điểm sao cho ta có thể tính nhanh giá trị của đa thức và trên đó, đồng thời có thể khôi phục được đa thức dựa trên .
¶ 2.3. Thuật toán
¶ Biến đổi xuôi
Ta có một đa thức . Không mất tính tổng quát, giả sử là một luỹ thừa của hay với . Nếu không phải là một lũy thừa của , ta thêm các số hạng bị thiếu và cho các hệ số bằng .
Tập số mà FFT chọn để tính là tập các căn đơn vị cấp , tức (nhắc lại ).
Định nghĩa. Cho một dãy , Biến đổi Fourier nhanh - Fast Fourier Transform - FFT là bất kì thuật toán nào tính dãy trong thời gian . Phép biến đổi này bản thân nó được gọi là Biến đổi Fourier rời rạc - Discrete Fourier Transform - DFT.
Bài viết này sẽ giới thiệu về một thuật toán FFT thông dụng là thuật toán Cooley-Tukey.
Ta sẽ biểu diễn bài toán dưới dạng ma trận:
Trước khi đi vào trường hợp tổng quát, ta sẽ xem xét một ví dụ với :
Các số mũ trên bảng đã được lấy dư cho . Ta sẽ đưa các cột màu đỏ (chỉ số chẵn) về bên trái, các cột màu xanh (chỉ số lẻ) về bên phải và chia ma trận thành bốn ma trận con như dưới đây. Chú ý rằng việc này cũng thay đổi thứ tự của vectơ hệ số.
Nhận thấy rằng các hệ số tương ứng ở hai ma trận con bên trái bằng nhau và các hệ số tương ứng ở hai ma trận con bên phải trái dấu, do tính chất .
Vậy biểu thức cần tính có dạng:
Để ý rằng:
Vậy ta chỉ cần tính và là đủ để tính kết quả. Mặt khác, có thể thấy tính và tương đương với tính và .
Xét trường hợp tổng quát, với , ta có (chú ý ở đây là chỉ số, không phải đơn vị ảo):
Nếu coi và thì hai biểu thức trên có thể viết lại thành:
Vậy ta cần tính và là hai bài toán với kích thước giảm đi một nửa.
Độ phức tạp. Thuật toán FFT là một thuật toán chia để trị nên ta dễ thấy nó có độ phức tạp .
¶ Biến đổi ngược
Bổ đề. .
Chứng minh:
Đặt và .
Ta có:
Xét , ta có:
bởi là một căn đơn vị.
Từ đây ta suy ra với .
Ngoài ra, nếu ta sử dụng thay cho ở lần thứ hai, kết quả ta nhận được sẽ là (bạn đọc tự chứng minh).
Ta sẽ cài đặt chung cả biến đổi xuôi và ngược:
using cd = complex<long double>;
void fft(vector<cd> &a, bool invert) {
/// invert = true tương ứng với biến đổi ngược
int n = a.size();
if (n == 1) return;
vector<cd> a0, a1;
for (int i = 0; i < n / 2; i++) {
a0.push_back(a[2 * i]);
a1.push_back(a[2 * i + 1]);
}
fft(a0, invert); fft(a1, invert);
cd w = 1, wn = polar(1.0L, acos(-1.0L) / n * (invert ? -2 : 2));
/// polar(r, t) = r * exp(it) và acos(-1.0L) = pi
/// thay wn = wn^-1 ở biến đổi ngược
for (int i = 0; i < n / 2; i++) {
a[i] = a0[i] + w * a1[i];
a[i + n / 2] = a0[i] - w * a1[i];
/// Ta sẽ chia 2 ở mỗi tầng đệ quy thay cho việc chia n ở cuối
if (invert) {
a[i] /= 2; a[i + n / 2] /= 2;
}
w *= wn;
}
}
Dưới đây là cài đặt để tính tích chập của hai dãy số:
vector<int> conv(const vector<int> &a, const vector<int> &b) {
if (a.empty() || b.empty()) return {};
vector<cd> fa(a.begin(), a.end());
vector<cd> fb(b.begin(), b.end());
int n = 1;
while (n < int(a.size() + b.size()) - 1) 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> res(n);
for (int i = 0; i < n; i++)
res[i] = int(real(fa[i]) + 0.5);
return res;
}
¶ Cài đặt khử đệ quy
Ta xét các tầng đệ quy với :
- Tầng :
- Tầng :
- Tầng :
Vậy nếu ta sắp xếp dãy lại thành thì mỗi lần gọi đệ quy ở tầng thứ sẽ thực hiện trên một đoạn dài .
Viết lại dãy chỉ số dưới dạng nhị phân:
Có thể thấy rằng nếu ta đảo ngược thứ tự bit thì, ví dụ thì ta nhận được dãy tăng dần từ đến .
Gọi là số nhận được sau khi đảo ngược thứ tự bit của , ta có với thể hiện phép toán bitwise OR (bạn đọc tự chứng minh).
Ta có cài đặt sau:
void fft(vector<cd> &a, bool invert) {
int n = a.size(), L = __builtin_ctz(n);
vector<int> rev(n);
for (int i = 0; i < n; i++) {
rev[i] = (rev[i >> 1] | (i & 1) << L) >> 1;
if (i < rev[i]) swap(a[i], a[rev[i]]);
}
for (int len = 2; len <= n; len <<= 1) {
cd wlen = polar(1.0L, acos(-1.0L) / len * (invert ? -2 : 2));
for (int i = 0; i < n; i += len) {
cd w = 1;
for (int j = 0; j < len / 2; j++) {
cd u = a[i + j];
cd v = a[i + j + len / 2] * w;
a[i + j] = u + v;
a[i + j + len / 2] = u - v;
w *= wlen;
}
}
}
if (invert) {
for (auto &x : a) x /= n;
}
}
¶ Vấn đề về độ chính xác
Ở cài đặt phía trên ta có viết w *= wlen để tính luỹ thừa của căn đơn vị. Việc nhân nhiều lần sẽ ảnh hưởng rất lớn đến độ chính xác của thuật toán vì ta đang thực hiện tính toán trên số thực.
Nhận xét rằng với mỗi ta chỉ cần tính .
Định nghĩa một mảng sao cho với mỗi và với mỗi ta có .
Cài đặt để tính mảng :
vector<cd> root(n);
root[1] = 1;
for (int k = 2; k < n; k *= 2) {
cd z = polar(1.0l, acos(-1.0l) / k * (invert ? 1 : -1));
for (int j = k / 2; j < k; j++) {
root[2 * j] = root[j];
root[2 * j + 1] = root[j] * z;
}
}
Sau khi có mảng ta có thể sửa cài đặt FFT như sau:
for (int k = 1; k < n; k *= 2)
for (int i = 0; i < n; i += 2 * k)
for (int j = 0; j < k; j++) {
cd z = root[j + k] * a[i + j + k];
a[i + j + k] = a[i + j] - z;
a[i + j] += z;
}
Thử nghiệm với , sai số chỉ rơi vào khoảng .
¶ FFT hai dãy cùng một lúc
Ta có thể tính của hai dãy và cùng một lúc bằng cách tính của dãy với .
Gọi lần lượt là các biến đổi tương ứng. Ta có:
kết hợp với suy ra:
Sử dụng công thức này ta có thể tính tích chập chỉ dùng hai lần gọi hàm fft:
vector<int> conv(const vector<int> &a, const vector<int> &b) {
if (a.empty() || b.empty()) return {};
int n = 1;
while (n < int(a.size() + b.size()) - 1) n <<= 1;
vector<cd> in(n), out(n);
for (int i = 0; i < int(a.size()); i++)
in[i].real(a[i]);
for (int i = 0; i < int(b.size()); i++)
in[i].image(b[i]);
fft(in, false);
for (int i = 0; i < n; i++)
in[i] *= in[i];
for (int i = 0; i < n; i++) {
/// (n - i) mod n
int j = -i & (n - 1);
out[i] = in[i] - conj(in[j]);
}
fft(out, true)
vector<int> res(n);
/// ở trên ta không chia cho 4i nên kết quả sẽ nằm trong phần ảo
for (int i = 0; i < n; i++)
res[i] = int(imag(out[i]) / 4 + 0.5);
return res;
}
¶ 3. Thuật toán NTT
Xét bài toán nhân đa thức nhưng lần này ta muốn các hệ số của đa thức chia lấy dư cho một số nguyên tố . Nếu ta sử dụng thuật toán FFT thông thường có thể gây ra sai số lớn vì hệ số của đa thức kết quả có thể rất lớn. Thuật toán NTT cho phép ta tính toán chỉ dùng số nguyên, từ đó kết quả luôn đảm bảo chính xác.
Thuật toán FFT dựa trên các tính chất của căn đơn vị. Các tính chất này cũng xuất hiện trên căn đơn vị trong số học modulo. Cụ thể ta gọi căn đơn vị cấp modulo là một số nguyên thoả mãn:
Điều kiện thứ hai có thể được viết lại thành: với mọi .
Các căn đơn vị cấp khác được biểu diễn bằng một luỹ thừa của .
Để áp dụng thuật toán FFT, ta cần các căn đơn vị cho các luỹ thừa nhỏ hơn của . Ta có thể chứng minh tính chất sau:
Từ đây ta có nếu là căn đơn vị cấp , thì là căn đơn vị cấp và do đó ta tính được căn đơn vị cho các luỹ thừa của nhỏ hơn.
Để tính biến đổi ngược ta cần nghịch đảo modulo tồn tại, điều này là hiển nhiên vì là số nguyên tố.
Ta chứng minh được rằng với một số nguyên tố có dạng , tồn tại với là một căn nguyên thuỷ modulo . Lấy ví dụ với có căn nguyên thuỷ và .
Cài đặt với (các hằng số root, root_pw thể hiện căn đơn vị và số mũ tương ứng, root_1 là nghịch đảo của căn đơn vị, hàm inverse tính nghịch đảo modulo của một số):
const int mod = 998244353;
const int root = 15311432;
const int root_1 = 469870224;
const int root_pw = 1 << 23;
void fft(vector<int> & a, bool invert) {
int n = a.size(), L = __builtin_ctz(n);
vector<int> rev(n);
for (int i = 0; i < n; i++) {
rev[i] = (rev[i >> 1] | (i & 1) << L) >> 1;
if (i < rev[i]) swap(a[i], a[rev[i]]);
}
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];
int v = 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 = 1ll * w * wlen % mod;
}
}
}
if (invert) {
int n_1 = inverse(n, mod);
for (int & x : a)
x = 1ll * x * n_1 % mod;
}
}
¶ Nhân đa thức với modulo bất kì
Có thể thấy rằng thuật toán NTT chỉ hoạt động với nếu căn đơn vị tồn tại. Với các modulo không thoả mãn ta có hai cách sau:
¶ Sử dụng định lý thặng dư Trung Hoa
Nếu kết quả của phép nhân đa thức có hệ số nhỏ hơn , với là hai số nguyên tố có dạng , ta có thể thực hiện NTT trên hai modulo này và dùng định lý thặng dư Trung Hoa để khôi phục kết quả.
¶ Chia nhỏ đa thức
Ta cần tính với hệ số modulo , thực hiện tách hai đa thức như sau:
với .
Khi đó tích của hai đa thức được biểu diễn thành:
Sử dụng kĩ thuật tương tự với FFT hai dãy cùng một lúc, ta có thể tính biểu thức trên với chỉ lần gọi FFT.
Cài đặt:
using ll = long long;
vector<int> convMod(const vector<int> &a, const vector<int> &b, int M) {
if (a.empty() || b.empty()) return {};
vector<int> res(a.size() + b.size() - 1);
int B = 32 - __builtin_clz(res.size());
int n = 1 << B, cut = sqrt(M);
vector<cd> L(n), R(n), outs(n), outl(n);
for (int i = 0; i < int(a.size()); i++)
L[i] = cd(a[i] / cut, a[i] % cut);
for (int i = 0; i < int(b.size()); i++)
R[i] = cd(b[i] / cut, b[i] % cut);
fft(L, false); fft(R, false);
for (int i = 0; i < n; i++) {
/// j = (n - i) % n
int j = -i & (n - 1);
outl[i] = (L[i] + conj(L[j])) * R[i] / 2.0l;
outs[i] = (L[i] - conj(L[j])) * R[i] / 2il;
}
fft(outl, true), fft(outs, true);
for (int i = 0; i < int(res.size()); i++) {
ll av = ll(real(outl[i]) + 0.5);
ll cv = ll(imag(outs[i]) + 0.5);
ll bv = ll(imag(outl[i]) + 0.5) + ll(real(outs[i]) + 0.5);
res[i] = ((av % M * cut + bv) % M * cut + cv) % M;
}
return res;
}