Fast Fourier Transform - Nhân đa thức nhanh

Bài toán nhân đa thức

Nhân hai đa thức là một thao tác cơ bản trong lập trình, nhất là trong mật mã. Giả sử chúng ta có hai đa thức:

• $A(x) = a_0 + a_1x + a_2x^2 + \dots + a_{n-1}x^{n-1}$

• $B(x) = b_0 + b_1x + b_2x^2 + \dots + b_{n-1}x^{n-1}$

Cách tiếp cận tự nhiên nhất mà chúng ta được học ở trường phổ thông là phép nhân phân phối (Elementary Multiplication). Để tìm hệ số của đa thức tích $C(x) = A(x) \cdot B(x)$, ta lấy từng hạng tử của $A$ nhân với mọi hạng tử của $B$. Trong khoa học máy tính, đây là độ phức tạp $O(N^2)$với $N,M$ lần lượt là bậc của đa thức $A,B$. Với $N,M \ge 1.000.000$ (trong các bài toán xử lý tín hiệu hoặc mã hóa), số phép toán lên tới 1.000 tỷ, khiến các hệ thống hiện đại cũng phải "đầu hàng" về mặt thời gian thực.

Sự xuất hiện của FFT (Fast Fourier Transform), cụ thể là thuật toán Cooley-Tukey, đã thay đổi hoàn toàn cuộc chơi này bằng cách đưa độ phức tạp về mức $O(N \log N)$ (để cho đơn giản, ta cứ giả định 2 đa thức $A,B$ có cùng bậc $N$, với $N$ là bậc lớn nhất trong $N$ và $M$). Ví dụ đưa 2 đa thức về cùng bậc 3 như sau:

Hệ số (Coefficient) vs. Điểm (Point-Value)

Thông thường, chúng ta định nghĩa đa thức qua các hệ số của nó. Nhưng thực ra, chúng ta vẫn có thể định nghĩa một đa thức thông qua các điểm phân biệt mà đồ thị của nó đi qua. Theo định lý cơ bản về đa thức: Một đa thức bậc $n-1$ được xác định duy nhất bởi tập hợp $n$ cặp điểm $(x_i, y_i)$ phân biệt.

Ví dụ:

• Qua 2 điểm, ta xác định được duy nhất 1 đường thẳng (bậc 1).

• Qua 3 điểm, ta xác định được duy nhất 1 đường parabol (bậc 2).

Hãy tưởng tượng bạn có hai đa thức $A$ và $B$ đã được chuyển sang dạng điểm tại cùng một tập hợp các giá trị $\{x_0, x_1, \dots, x_{2n-1}\}$.

Để tìm đa thức tích $C = A \cdot B$ ở dạng điểm, chúng ta chỉ cần tính:

Phép toán này chỉ tốn đúng $2n$ phép nhân, tức là độ phức tạp $O(N)$, yeah!!!

Vậy tại sao nhanh như vậy mà ta không dùng nó ngay từ đầu? Đó là vì chi phí chuyển đổi:

1. Chuyển từ Hệ số $\rightarrow$ Điểm (Evaluation): Nếu tính toán thông thường (thế $x$ vào đa thức), ta mất $O(N^2)$.

2. Chuyển từ Điểm $\rightarrow$ Hệ số (Interpolation): Dùng các phương pháp như Newton hay Lagrange cũng mất $O(N^2)$.

FFT cho phép chúng ta thực hiện cả hai quá trình chuyển đổi trên chỉ với độ phức tạp $O(N \log N)$.

Thuật toán Cooley-Tukey

Để chuyển một đa thức từ dạng Hệ số sang dạng Điểm trong $O(N \log N)$, thuật toán Cooley-Tukey sử dụng chiến lược Chia để trị (Divide and Conquer) dựa trên một lựa chọn toán học cực kỳ thông minh: Căn đơn vị phức (Complex Roots of Unity).

Căn đơn vị phức

Nếu ta chọn các điểm $x$ ngẫu nhiên, ta vẫn mất $O(N^2)$ để tính giá trị đa thức. Tuy nhiên, nếu ta chọn $x$ là các nghiệm của phương trình $x^N = 1$, chúng ta có những tính chất đối xứng tuyệt vời. Các điểm này nằm trên vòng tròn đơn vị trong mặt phẳng phức, được xác định bởi công thức:

Với $k = 0, 1, \dots, N-1$.

Tính chất: $(\omega_N^k)^2 = \omega_{N/2}^k$. Điều này cho phép ta thu nhỏ bài toán kích thước $N$ thành bài toán kích thước $N/2$ một cách đồng nhất.

Chiến lược Chia để trị: Tách Chẵn và Lẻ

Giả sử ta có đa thức $A(x)$ bậc $N-1$(với $N$ là lũy thừa của 2).

$A(x) = a_0 + a_1x + a_2x^2 + a_3x^3 + a_4x^4 + a_5x^5 + \dots$

Chúng ta nhóm các số hạng có số mũ chẵn vào một nhóm, và các số hạng có số mũ lẻ vào một nhóm:

• Nhóm chẵn: $a_0 + a_2x^2 + a_4x^4 + \dots$

• Nhóm lẻ: $a_1x + a_3x^3 + a_5x^5 + \dots = x(a_1 + a_3x^2 + a_5x^4 + \dots)$

Bây giờ, thay vì giữ nguyên biến $x$, chúng ta đặt $y = x^2$. Khi đó:

• Đa thức chứa các hệ số chẵn trở thành: $A_{even}(y) = a_0 + a_2y + a_4y^2 + \dots$

• Đa thức chứa các hệ số lẻ trở thành: $A_{odd}(y) = a_1 + a_3y + a_5y^2 + \dots$

Khi thay $y = x^2$ ngược lại vào hai đa thức trên, ta thấy:

• $A_{even}(x^2) = a_0 + a_2x^2 + a_4x^4 + \dots$ (Chính là phần chẵn của $A(x)$)

• $x \cdot A_{odd}(x^2) = x(a_1 + a_3x^2 + a_5x^4 + \dots) = a_1x + a_3x^3 + a_5x^5 + \dots$ (Chính là phần lẻ của $A(x)$)

Vậy, ta có biểu diễn cuối cùng:

Cấu trúc cánh bướm

Trong FFT, chúng ta không thay $x$ bằng các số ngẫu nhiên, mà thay $x$ bằng các căn đơn vị phức $\omega_N^k$. Chìa khóa nằm ở tính chất: khi bình phương một căn đơn vị, bậc của nó giảm đi một nửa.

Chúng ta có công thức như đã đề cập ở trên: $(\omega_N^k)^2 = \omega_{N/2}^k$.

Để tính $A(\omega_N^k)$ cho mọi $k \in [0, N-1]$, ta sử dụng công thức tách chẵn-lẻ:

Khi tính toán cho nửa sau của mảng (tức là khi chỉ số $k$ nằm trong khoảng $N/2 \leq k < N$), ta đặt $k = m + N/2$ (với $0 \leq m < N/2$). Lúc này, công thức trở thành:

• Với nửa đầu ($0 \leq k < N/2$):

• Với nửa sau ($k = m + N/2$, $0 \leq m < N/2$):

Thay $k = m + N/2$ vào công thức tổng quát, ta có:

Tại đây, có 2 biến đổi quan trọng:

• Bình phương: $(\omega_N^{m + N/2})^2 = \omega_N^{2m + N} = \omega_N^{2m} \cdot \omega_N^N$. Vì $\omega_N^N = 1$, nên nó bằng $\omega_N^{2m} = \omega_{N/2}^m$.

• Số hạng xoay (Twiddle factor): $\omega_N^{m + N/2} = \omega_N^m \cdot \omega_N^{N/2}$. Mà $\omega_N^{N/2} = -1$ (đây là tính chất chu kỳ của số phức trên vòng tròn đơn vị), nên nó bằng $-\omega_N^m$.

Thay vào biểu thức trên, ta được kết quả cuối cùng:

Cả hai nửa đầu và nửa sau đều dùng chung đúng hai giá trị là $A_{even}(m)$ và $A_{odd}(m)$.

• Để tính nửa đầu, lấy Tổng.

• Để tính nửa sau, lấy Hiệu.

Thay vì phải chạy vòng lặp để tính toàn bộ $A(x)$ mất $O(N^2)$, giờ đây với mỗi cặp $(k, m+N/2)$, ta chỉ cần tính $A_{even}$ và $A_{odd}$ một lần, sau đó cộng/trừ là ra được ngay 2 kết quả cho tầng trên. Đây chính là lý do giúp thuật toán giảm độ phức tạp xuống $O(N \log N)$.

Kết quả: Mỗi bước tính toán ở tầng dưới cung cấp dữ liệu cho hai giá trị ở tầng trên chỉ bằng một phép nhân và hai phép cộng/trừ. Đây chính là cấu trúc Butterfly nổi tiếng giúp giảm số lượng phép tính từ $N^2$ xuống còn $N \log_2 N$.

Bức ảnh này giải thích chính xác cơ chế "phép màu" đã được đề cập:

• Đầu vào (bên trái): Hai giá trị $A[k]$ và $A[k+N/2]$ từ tầng tính toán trước đó.

• Hệ số xoay (Twiddle Factor - $W_N^k$): Giá trị phức được nhân với đầu vào thứ hai trước khi kết hợp.

• Điểm kết hợp (Cánh bướm):

  • Phép cộng (trên): Tạo ra đầu ra thứ nhất: $A[k]_{next} = A[k] + W_N^k A[k+N/2]$.

  • Phép trừ (dưới): Tạo ra đầu ra thứ hai: $A[k+N/2]_{next} = A[k] - W_N^k A[k+N/2]$.

Cấu trúc đối xứng này cho phép máy tính tái sử dụng các kết quả tính toán trung gian, biến đổi độ phức tạp từ $O(N^2)$ thành $O(N \log N)$.

Quy trình nhân đa thức

Để nhân hai đa thức $A$ và $B$, quy trình gồm 3 bước:

Bước 1: Mở rộng và Biến đổi (Padding & Forward FFT)

Trước khi bắt đầu, có một chi tiết kỹ thuật cực kỳ quan trọng: Xử lý bậc của đa thức kết quả.

Nếu $A$ có bậc $n-1$ và $B$ có bậc $m-1$, thì đa thức tích $C$ sẽ có bậc $(n-1) + (m-1)$.

• Zero Padding: Chúng ta phải thêm các hệ số $0$ vào $A$ và $B$ sao cho tổng số hệ số của mỗi đa thức ít nhất bằng bậc của $C+1$, và phải là một lũy thừa của 2 (ví dụ: $2^k$).

• Thực thi: Áp dụng FFT Cooley-Tukey để chuyển $A$ và $B$ từ dạng hệ số sang dạng điểm (tại các căn đơn vị $\omega_N^k$).

  • $A \xrightarrow{FFT} \{A(\omega_N^0), A(\omega_N^1), \dots, A(\omega_N^{N-1})\}$

  • $B \xrightarrow{FFT} \{B(\omega_N^0), B(\omega_N^1), \dots, B(\omega_N^{N-1})\}$

Bước 2: Nhân trong miền giá trị (Point-wise Multiplication)

Với các giá trị điểm đã có, ta tính giá trị của đa thức tích $C$ tại từng điểm tương ứng:

Chỉ cần đúng $N$ phép nhân số phức, chúng ta đã có đầy đủ đại diện dạng điểm của đa thức tích.

Bước 3: Biến đổi ngược (Inverse FFT - IFFT)

Bây giờ, chúng ta đang có các giá trị $C(\omega_N^k)$, nhưng mục tiêu là phải xác định kết quả là các hệ số $\{c_0, c_1, \dots, c_{N-1}\}$. Lúc này, ta cần đến quá trình biến đổi ngược Inverse FFT. Thú vị thay, thuật toán cho IFFT gần như y hệt FFT. Sự khác biệt duy nhất nằm ở hai điểm:

• Thay thế hệ số xoay $\omega_N^k$ bằng nghịch đảo của nó (tương đương với việc đổi dấu phần ảo trong số phức).

• Sau khi tính xong, tất cả kết quả phải được chia cho $N$

Tổng kết độ phức tạp

Hãy nhìn lại toàn bộ hành trình:

1. FFT cho A và B: $2 \times O(N \log N)$

2. Nhân điểm: $O(N)$

3. IFFT cho C: $O(N \log N)$

Tổng cộng: $\mathbf{O(N \log N)}$.

Phân tích Pseudo Code

Để cài đặt FFT Cooley-Tukey, chúng ta sử dụng cấu trúc đệ quy chia mảng hệ số thành hai nửa: vị trí chẵn và vị trí lẻ.

1. Hàm FFT (Hệ số $\rightarrow$ Điểm)

HÀM FFT(Mảng_Hệ_Số A):
    N = Độ dài của A
    NẾU N == 1: 
        TRẢ VỀ A  // Trường hợp cơ bản: Đa thức bậc 0 là chính nó

    // BƯỚC 1: CHIA (Divide)
    A_chẵn = [A[0], A[2], ..., A[N-2]]
    A_lẻ   = [A[1], A[3], ..., A[N-1]]

    // BƯỚC 2: TRỊ (Conquer - Đệ quy)
    Y_chẵn = FFT(A_chẵn)
    Y_lẻ   = FFT(A_lẻ)

    // BƯỚC 3: KẾT HỢP (Combine - Cấu trúc Cánh bướm)
    Mảng_Kết_Quả Y = Mảng rỗng kích thước N
    Hệ_số_xoay_gốc = e^(i * 2π / N)
    w = 1

    CHO k CHẠY TỪ 0 ĐẾN N/2 - 1:
        // Tính toán dựa trên tính đối xứng và chu kỳ
        Y[k]       = Y_chẵn[k] + w * Y_lẻ[k]
        Y[k + N/2] = Y_chẵn[k] - w * Y_lẻ[k]
        
        w = w * Hệ_số_xoay_gốc // Xoay dần trên vòng tròn đơn vị

    TRẢ VỀ Y

2. Quy trình nhân hai đa thức $P_1$ và $P_2$

HÀM NHÂN_ĐA_THỨC(P1, P2):
    1. Padding: Thêm các số 0 vào P1, P2 sao cho độ dài N là lũy thừa của 2 
       (N phải đủ lớn để chứa bậc của đa thức kết quả).

    2. Chuyển đổi sang dạng điểm:
       Điểm_P1 = FFT(P1)
       Điểm_P2 = FFT(P2)

    3. Nhân trực tiếp (Point-wise):
       Điểm_Kết_Quả = Mảng chứa (Điểm_P1[i] * Điểm_P2[i]) với mọi i từ 0 đến N-1

    4. Chuyển ngược về dạng hệ số:
       Hệ_Số_Kết_Quả = IFFT(Điểm_Kết_Quả) 
       (IFFT dùng thuật toán y hệt FFT nhưng đổi dấu số phức và chia kết quả cho N)

    TRẢ VỀ Hệ_Số_Kết_Quả

Demo