Biến đổi Fourier rời rạc

Phần này mình dịch từ [11].

Giới thiệu

Cho số \(t \in \mathbb{N}\) và \(n = 2^t\).

Gọi \(R\) là vành với đơn vị \(1\). Vành \(R\) chứa phần tử \(2^{-1}\) là nghịch đảo của phần tử \(2\). Hơn nữa vành \(R\) cũng chứa phần tử \(\zeta_{2n}\) là một nghiệm cố định nào đó của phương trình \(x^n + 1 = 0\). Đặt \(\zeta_n = \zeta_{2n}^2\).

Ví dụ 13

Trên \(\mathbb{C}\) thì \(\zeta_{2n} = e^{\pi i/n}\) là một nghiệm của phương trình \(x^n + 1 = 0\). Khi đó \(\zeta_n = \zeta_{2n}^2 = e^{2\pi i /n}\).

Bổ đề 1

Phần tử \(\zeta_{2n}\) sinh ra một nhóm vòng (cyclic group) với order \(2n\) trong nhóm nhân của vành \(R\).

Chứng minh

Do \(\zeta_{2n}\) là nghiệm của phương trình \(x^n + 1 = 0\) nên \(\zeta_{2n}^n = -1\), hay \(\zeta_{2n}^{2^t} = -1\). Từ đó suy ra \(\left(\zeta_{2n}^{2^t}\right)^2 = 1\) nên order của phần tử \(\zeta_{2n}\) là \(2^{t+1} = 2n\).

Giả sử \((f_0, \ldots, f_{n-1}) \in R^n\) là vector bất kì. Khi đó:

Biến đổi Fourier loại 1 của vector trên là vector \(n\) chiều \((\hat f_0, \ldots, \hat f_{n-1})\) thuộc \(R^n\) xác định bởi:

\[\hat f_i = \sum_{j=0}^{n-1} \zeta_n^{ij} f_j, \quad i = 0, 1, \ldots, n-1.\]

Biến đổi Fourier loại 2 của vector trên là vector \(n\) chiều \((\check f_1, \check f_3, \ldots, \check f_{2n-1})\) thuộc \(R^n\) xác định bởi:

\[\check f_i = \sum_{j=0}^{n-1} \zeta_{2n}^{ij} f_j, \quad i = 1, 3, \ldots, 2n-1.\]

Ví dụ 14

Mình lấy ví dụ với vector \((f_0, f_1, f_2) = (1, 2, 3)\) thuộc \(\mathbb{C}^3\).

Chọn nghiệm \(\zeta_{2n} = e^{\pi i / 3}\) của phương trình \(x^3 + 1 = 0\). Khi đó \(\zeta_n = \zeta_{2n}^2 = e^{2\pi i / 3}\).

Biến đổi Fourier loại 1 lúc này là:

\[\begin{split}\hat{f}_0 & = \left(e^{2\pi i/3}\right)^{0 \cdot 0} \cdot 1 + \left(e^{2\pi i/3}\right)^{0 \cdot 1} \cdot 2 + \left(e^{2\pi i/3}\right)^{0 \cdot 2} \cdot 3 \\ & = 1 + 2 + 3 = 6. \\ \hat{f}_1 & = \left(e^{2\pi i/3}\right)^{1 \cdot 0} \cdot 1 + \left(e^{2\pi i/3}\right)^{1 \cdot 1} \cdot 2 + \left(e^{2\pi i/3}\right)^{1 \cdot 2} \cdot 3 \\ & = 1 + 2 \cdot e^{2\pi i/3} + 3 \cdot 3 \cdot e^{4\pi i/3} \\ \hat{f}_2 & = \left(e^{2\pi i/3}\right)^{2 \cdot 0} \cdot 1 + \left(e^{2\pi i/3}\right)^{2 \cdot 1} \cdot 2 + \left(e^{2\pi i/3}\right)^{2 \cdot 2} \cdot 3 \\ & = 1 + 2 \cdot e^{4\pi i/3} + 3 \cdot e^{2\pi i/3}.\end{split}\]

Chú ý rằng \(e^{2\pi i} = 1\) nên biểu thức cuối của \(\hat{f}_2\) rút gọn còn \(e^{2\pi i / 3}\) (\(8 \equiv 2 \bmod 3\)).

Tiếp theo ta biểu diễn \(e^{4\pi i / 3}\) theo \(e^{2\pi i / 3}\) như sau. Ta xét

\[\begin{split}e^{2zi} & = \cos 2z + i\sin 2z \\ & = 2\cos^2 z - 1 + 2i \sin z \cos z \\ & = 2\cos z (\cos z + i\sin z) - 1 \\ & = 2\cos z \cdot e^{zi} - 1.\end{split}\]

Thay \(z = 2\pi / 3\) vào kết quả trên ta được

\[e^{4\pi i / 3} = 2\cos \frac{2\pi}{3} \cdot e^{2\pi i/3} - 1 = -e^{2\pi i / 3} - 1.\]

Như vậy kết quả của biến đổi Fourier rời rạc bên trên là

\[(\hat{f}_0, \hat{f}_1, \hat{f}_2) = (6, -e^{2\pi i / 3} - 2, e^{2\pi i / 3} - 1).\]

Ở ví dụ trên chúng ta tính biến đổi Fourier cho vector có độ dài không phải lũy thừa của \(2\). Sau đây chúng ta sẽ xem xét một ví dụ với vector có độ dài là lũy thừa của \(2\) và ở phần sau sẽ tối ưu tính toán với biến đổi Fourier nhanh.

Ví dụ 15

Mình lấy ví dụ với vector \((f_0, f_1, f_2, f_3) = (1, 2, 3, 4)\) thuộc \(\mathbb{C}^4\).

Chọn nghiệm \(\zeta_{2n} = e^{i\pi / 4}\) của phương trình \(x^4 + 1 = 0\). Khi đó \(\zeta_n = \zeta_{2n}^2 = e^{i \pi / 2} = i\).

Biến đổi Fourier loại 1 lúc này là:

\[\begin{split}\hat{f}_0 & = 1 \cdot 1 + 2 \cdot 1 + 3 \cdot 1 + 4 \cdot 1 = 10, \\ \hat{f}_1 & = 1 \cdot 1 + 2 \cdot i + 3 \cdot (-1) + 4 \cdot (-i) = -2 - 2i, \\ \hat{f}_2 & = 1 \cdot 1 + 2 \cdot (-1) + 3 \cdot 1 + 4 \cdot (-1) = -2, \\ \hat{f}_3 & = 1 \cdot 1 + 2 \cdot (-i) + 3 \cdot (-1) + 4 \cdot i = -2 + 2i.\end{split}\]

Đối với vector độ dài \(n = 2^2\) thì \(t = 2\), chúng ta cần \(n - 1\) phép cộng và \(n\) phép nhân với lũy thừa của \(\zeta_{n}\) để tính mỗi \(\hat{f}_i\), với \(i = 0, 1, 2, 3\). Tổng quát, ta cần tất cả \((n + n - 1) \cdot n\) phép tính cộng và nhân với lũy thừa của \(\zeta_n\). Độ phức tạp để tính biến đổi Fourier rời rạc là \(O(n^2)\).

Nhận xét 8

Các phần tử \(\hat f_i\) là giá trị của đa thức \(F(x) = \sum\limits_{j=0}^{n-1} f_j x^j \in R[x]\) tại điểm \(x = \zeta_n^i\), các phần tử \(\check f_i\) là giá trị của đa thức \(F(x)\) tại \(x = \zeta_{2n}^i\) khi \(i\) lẻ.

Bổ đề 2

Ta có biểu diễn ngược như sau:

(1)\[f_i = n^{-1} \sum_{j=0}^{n-1} \hat f_j \zeta_n^{-ij},\]

(2)\[\begin{split}f_i = n^{-1} \sum_{\substack{1 \leqslant j \leqslant 2n-1 \\ j \ \text{lẻ}}} \check f_j \zeta_{2n}^{-ij},\end{split}\]

với \(n^{-1} = (2^{-1})^t \in R\).

Chứng minh

Đặt \(k \in \mathbb{Z}\). Khi đó \(\zeta_{2n}^{-k} = \zeta_{2n}^{2nr - k}\) với \(r \in \mathbb{Z}\), \(2nr - k \geqslant 0\).

Ngoài ra, vì order của phần tử \(\zeta_n\) trong phép nhân là \(n\) nên ta có:

(3)\[\sum_{j = 0}^{n - 1} \zeta_n^{lj} = n \quad \text{khi} \quad l \equiv 0 \pmod n,\]

(4)\[\sum_{j = 0}^{n - 1} \zeta_n^{lj} = 0 \quad \text{khi} \quad l \not\equiv 0 \pmod n.\]

Phương trình (1) thỏa mãn vì:

\[\sum_{j=0}^{n-1} \hat f_j \zeta_n^{-ij} = \sum_{j=0}^{n-1} \sum_{k=0}^{n-1} f_k \zeta_n^{jk} \zeta_n^{-ij} = \sum_{k=0}^{n-1} f_k \sum_{j=0}^{n-1} \zeta_n^{(k-i)j} = n f_i\]

theo (3) và (4), cụ thể là khi \(k - i \equiv 0 \bmod n\).

Đẳng thức (2) có được từ

\[\begin{split}\sum_{\substack{1 \leqslant j \leqslant 2n - 1 \\ j \ \text{lẻ}}} \check f_j \zeta_{2n}^{-ij} & = \sum_{k=0}^{n-1} \check f_{2k+1} \zeta_{2n}^{-i(2k+1)} \\ & = \sum_{k=0}^{n-1} \sum_{l=0}^{n-1} f_l \zeta_{2n}^{l(2k+1)} \zeta_{2n}^{-i(2k+1)} \\ & = \sum_{l=0}^{n-1} f_l \zeta_{2n}^{l-i} \sum_{k=0}^{n-1} \zeta_n^{(l-i) k}.\end{split}\]

Từ (3) và (4) suy ra được đẳng thức (2), cũng là khi \(l - i \equiv 0 \bmod n\).

Việc tính toán biến đổi Fourier rời rạc bằng (1) và (2) cần \(O(n^2)\) phép tính cộng và nhân trong vành \(R\). Phần sau sẽ trình bày phương pháp tính biến đổi Fourier rời rạc với độ phức tạp \(O(n \log n)\). Phương pháp này được gọi là biến đổi Fourier nhanh (hay Fast Fourier Transform, FFT, быстрое преобразование Фурье).

Biến đổi Fourier nhanh

Định lí 6

Biến đổi Fourier rời rạc \((\hat f_0, \ldots, \hat f_{n-1})\) có thể được tính bằng:

• \(nt\) phép cộng trong vành \(R\);

• \(nt\) phép nhân với lũy thừa của \(\zeta_n\) trong vành \(R\).

Biến đổi Fourier rời rạc \((\check f_1, \ldots, \check f_{2n-1})\) có thể được tính bằng:

• \(nt\) phép cộng trong vành \(R\);

• \(nt\) phép nhân với lũy thừa của \(\zeta_{2n}\) trong vành \(R\).

Nếu có \((\hat f_i)\) và \((\check f_i)\) thì có thể tính được vector \((f_0, \ldots, f_{n-1})\) bằng:

• \(nt\) phép cộng trong vành \(R\);

• \(nt\) phép nhân với lũy thừa của \(\zeta_{2n}\) trong vành \(R\);

• \(n\) phép nhân với \(n^{-1} \in R\).

Chứng minh

Đặt

\[\begin{split}F(x) = \sum_{j=0}^{n-1} f_i x^i & = \sum_{\substack{0 \leqslant j \leqslant n-1 \\ j \ \text{chẵn}}} f_j x^j + \sum_{\substack{0 \leqslant j \leqslant n-1 \\ j \ \text{lẻ}}} f_j x^j \\ & = F_0(x^2) + x F_1(x^2),\end{split}\]

với \(\deg F_0(x), \deg F_1(x) < \dfrac{n}{2} = 2^{t-1}\).

Khi đó

(5)\[F(\zeta_n^i) = F_0(\zeta_n^{2i}) + \zeta_n^i F_1(\zeta_n^{2i}),\]

với \(i = 0, \ldots, n-1\).

Đặt \(\zeta_{n/2} = \zeta_n^2\). Khi đó

\[\left\{ \zeta_n^{2i} : 0 \leqslant i \leqslant n - 1 \right\} = \left\{ \zeta_{n/2}^i : 0 \leqslant i \leqslant \dfrac{n}{2} - 1\right\}.\]

Dựa trên quy nạp ta sẽ chứng minh biến đổi Fourier rời rạc loại 1, tức là tính \((F(\zeta_n^0), \ldots, F(\zeta_n^{n-1}))\), theo công thức (5) bằng \(n\) phép cộng trong vành \(R\) và \(n\) phép tính nhân với lũy thừa của \(\zeta_n\) trong vành \(R\) nếu ta đã biết các giá trị \(F_0(\zeta_{n/2}^i)\) và \(F_1(\zeta_{n/2}^i)\) với \(i = 0, \ldots, \dfrac{n}{2} - 1\).

Khi \(t = 1\) và \(n = 2 = 2^t\) (bước cơ sở quy nạp) thì để tính \((\hat f_0, \hat f_1)\) ta cần tìm \(F_0 + \zeta_2^i F_1\) với \(F_0\) và \(F_1\) là các phần tử thuộc vành \(R\), \(i = 0, 1\). Như vậy, với \(n = 2\) thì cần:

• \(2 = nt\) phép nhân với lũy thừa của \(\zeta_n = \zeta_2\) trong vành \(R\), ứng với \(\zeta_2^i F_1\) khi \(i = 0, 1\);

• \(2 = nt\) phép cộng trong vành \(R\), ứng với \(F_0 {\color{red}+} \zeta_2^i F_1\).

Giả sử với mọi \(j < t\), để tính biến đổi Fourier rời rạc loại 1 trên vector \(2^j\) chiều ta cần:

• \(2^j \cdot j\) phép nhân với lũy thừa của \(\zeta_{2^j} = \left(\zeta_n\right)^{2^{t-j}}\) trong vành \(R\);

• \(2^j \cdot j\) phép cộng trong vành \(R\).

Khi đó, nếu \(j = t\) thì việc tính \((\hat f_0, \ldots, \hat f_{n-1})\) theo công thức (5) bao gồm tính \(F_0(\zeta_{n/2}^i)\) và \(F_1(\zeta_{n/2}^i)\) với \(i = 0, \ldots, n-1\). Nói cách khác là tính biến đổi Fourier cho vector độ dài \(n / 2\) gồm các hệ số \(F_0(x)\) và \(F_1(x)\) cộng thêm:

• \(n\) phép cộng trong vành \(R\);

• \(n\) phép nhân với lũy thừa của \(\zeta_n\) trong vành \(R\).

Khi đó, theo giả thiết quy nạp, việc tính vector \((\hat f_0, \ldots, \hat f_{n-1})\) cần không nhiều hơn:

• \(n\) phép cộng trong vành \(R\);

• cộng với \(n\) phép nhân trong vành \(R\);

• cộng với \(2 \cdot 2^{t-1} \cdot (t-1)\) phép nhân với lũy thừa của \(\zeta_n\) trong vành \(R\).

Như vậy cần tổng cộng \(n + n(t - 1) = nt\) phép cộng trong vành \(R\) và \(n + n(t-1) = nt\) phép nhân với lũy thừa của \(\zeta_n\) trong vành \(R\). Như vậy mệnh đề đầu tiên của định lý được chứng minh.

Mệnh đề thứ hai chứng minh tương tự, ta thay \(\zeta_n\) thành \(\zeta_{2n}\).

Mệnh đề thứ ba suy ra từ hai mệnh đề trên và Bổ đề 2.

Ví dụ 16

Xét vector \((f_0, f_1, f_2, f_3) = (1, 2, 3, 4)\) thuộc \(\mathbb{C}^4\) như bên trên.

Ở đây, \(n = 4 = 2^2\) và \(t = 2\).

Chọn nghiệm \(\zeta_{8} = e^{i\pi / 4}\) của phương trình \(x^4 + 1 = 0\). Khi đó \(\zeta_4 = \zeta_{8}^2 = e^{i \pi / 2} = i\).

Ta chia vector thành hai phần:

• vector con gồm các phần tử ở vị trí chẵn là \((f_0, f_2) = (1, 3)\);

• vector con gồm các phần tử ở vị trí lẻ là \((f_1, f_3) = (2, 4)\).

Đặt \(\zeta_2 = \zeta_{4}^2 = -1\).

Khi đó, với vector con \((f_0, f_2) = (1, 3)\) ta tính:

\[\begin{split}F_0 & = f_0 + \zeta_{2}^0 \cdot f_2 = 1 + 1 \cdot 3 = 4, \\ F_1 & = f_0 + \zeta_{2}^1 \cdot f_2 = 1 + (-1) \cdot 3 = -2.\end{split}\]

Tương tự, với vector con \((f_1, f_3) = (2, 4)\) ta tính

\[\begin{split}F_2 & = f_1 + \zeta_{2}^0 \cdot f_3 = 2 + 1 \cdot 4 = 6, \\ F_3 & = f_1 + \zeta_{2}^1 \cdot f_3 = 2 + (-1) \cdot 4 = -2.\end{split}\]

Kết hợp hai vector \((F_0, F_2)\) và \((F_3, F_4)\) ta tính được biến đổi Fourier nhanh.

\[\begin{split}\hat{f}_0 & = F_0 + \zeta_4^0 \cdot F_2 = 4 + 1 \cdot 6 = 10, \\ \hat{f}_2 & = F_0 + \zeta_4^2 \cdot F_2 = 4 + (-1) \cdot 6 = -2, \\ \hat{f}_1 & = F_1 + \zeta_4^1 \cdot F_3 = -2 + i \cdot (-2) = -2 - 2 i, \\ \hat{f}_3 & = F_1 + \zeta_4^3 \cdot F_3 = -2 + (-i) \cdot (-2) = -2 + 2i.\end{split}\]

Ở đây, hệ số \(\hat{f}_i\) ở vị trí lẻ sẽ được tính nhờ các \(F_i\) ở vị trí lẻ cùng với lũy thừa lẻ của \(\zeta_n\), và ngược lại, hệ số \(\hat{f}_i\) ở vị trí chẵn sẽ được tính nhờ các \(F_i\) ở vị trí chẵn cùng với lũy thừa chẵn của \(\zeta_n\).

Như vậy kết quả là

\[(\hat{f}_0, \hat{f}_1, \hat{f}_2, \hat{f}_3) = (10, -2 - 2i, -2, -2 + 2i),\]

giống với ví dụ ở trên. Tuy nhiên chúng ta chỉ dùng \(4\) phép cộng và \(4\) phép nhân với \(\zeta_{2}\) để tính các \(F_i\), cộng thêm \(4\) phép cộng và \(4\) phép nhân với \(\zeta_4\) để tính các \(\hat{f}_i\). Như vậy tổng cộng ta dùng \(8 = n t\) phép cộng, và \(8 = n t\) phép nhân cho các lũy thừa của \(\zeta_4\) (vì \(\zeta_{2}\) cũng là lũy thừa của \(\zeta_4\)).

Biến đổi Fourier và phép nhân đa thức

Bổ đề 3

Một phép nhân trong vành \(R[x]/(x^{2^t} + 1)\) có thể được tính bởi:

• \(n = 2^t\) phép nhân trong vành \(R\);

• \(3nt\) phép cộng trong vành \(R\);

• \(3nt\) phép nhân với lũy thừa của \(\zeta_{2n}\) trong vành \(R\);

• \(n\) phép nhân với \(n^{-1}\) trong vành \(R\).

Nhận xét 9

Bổ đề không áp dụng khi \(R = \mathbb{Z}\) và \(R = \mathbb{Q}\) vì hai vành này không có phần tử \(\zeta_{2n}\) (căn bậc \(n\) của đơn vị).

Sau đây ta chứng minh Bổ đề 3.

Chứng minh

Đặt \(\sum\limits_{i=0}^{n-1} f_i x^i\) và \(\sum\limits_{i=0}^{n-1} g_i x^i\) với \(F, G \in R[x]\) và là biểu diễn của phần tử trong lớp nhân tử \(R[x]/(x^{2^t} + 1)\).

Đặt \(H = \sum\limits_{i=0}^{n-1} h_i x^i \in R[x]\) sao cho \(F \cdot G \equiv H \pmod{x^n+1}\). Nói cách khác \(H\) là tích \(F \cdot G\) trong vành và ta sẽ tính \(H\).

Biến đổi Fourier loại 2 cho vector \((f_0, \ldots, f_{n-1})\) và \((g_0, \ldots, g_{n-1})\) cho ta đẳng thức:

\[\check f_i \cdot \check g_i = F(\zeta_{2n}^i) \cdot G(\zeta_{2n}^i) = H(\zeta_{2n}^i) = \check h_i\]

với \(i\) lẻ, \(1 \leqslant i \leqslant 2n-1\), vì \(\zeta_{2n}^n + 1 = 0\).

Như vậy nếu ta biết tất cả \(\check f_i\) và \(\check g_i\) thì có thể tính mọi \(\check h_i\) với \(n\) phép nhân trong vành \(R\).

Theo Định lí 6, các phần tử \(\check f_i\) và \(\check g_i\) có thể tính với \(2tn\) phép cộng trong vành \(R\) và \(2tn\) phép nhân với lũy thừa của \(\zeta_{2n}\) trong vành \(R\).

Cũng theo Định lí 6, các phần tử \(h_i\) có thể tính, với điều kiện đã biết \(\check h_i\), với

• \(tn\) phép cộng trong vành \(R\);

• \(tn\) phép nhân với lũy thừa của \(\zeta_{2n}\) trong vành \(R\)

• \(n\) phép nhân với phần tử \(n^{-1}\) trong vành \(R\).

Ta có điều phải chứng minh.

Hệ quả 1

Đặt \(T\) là vành giao hoán với đơn vị, \(2^{-1} \in T\), \(\zeta_{4n} = \zeta_{2^{t + 2}} \in T\) là nghiệm phương trình \(x^{2n} + 1 = 0\). Nếu \(F(x), G(x) \in T[x]\), \(\deg F(x) < n\), \(\deg G(x) < n\) thì tích \(F(x) \cdot G(x)\) có thể tính với:

• \(2n\) phép nhân trong \(T\);

• \(6n(t+1)\) phép cộng trong \(T\);

• \(6n(t+1)\) phép nhân với lũy thừa của \(\zeta_{4n}\) trong \(T\);

• \(2n\) phép nhân với \((2^{-1})^{t+1}\) trong \(T\).

Chứng minh

Vì chặn trên của bậc đa thức \(F(x)\) và \(G(x)\) nên bậc của \(F(x) \cdot G(x)\) sẽ nhỏ hơn \(2n\). Khi chia \(F(x) \cdot G(x)\) cho \(2^{2n} + 1\) ta thu được chính đa thức \(F(x) \cdot G(x)\). Khi đó tích \(F(x) \cdot G(x)\) trong \(R[x]\) có thể tính với một phép nhân trong \(R[x] / (x^{2n} + 1)\). Từ Bổ đề 1, ta thay \(n\) thành \(2n\) và \(t\) thay thành \(t+1\), kết hợp Bổ đề 3 ta có điều phải chứng minh.

Đặt \(T\) là vành giao hoán với đơn vị \(1\) và có phần tử \(2^{-1}\).

Ta hiểu phép cộng trong \(T\) là cả phép cộng lẫn phép trừ (cộng cho số đối). Ta cần bổ đề cơ sở sau.

Bổ đề 4

Nếu \(t \geqslant 2\) thì một phép nhân trong vành \(T[x] / (x^{2^t} + 1)\) có thể thực hiện với:

• \(O(2^t \cdot t)\) phép nhân trong \(T\);

• \(O(2^t \cdot t \cdot \log t)\) phép cộng trong \(T\).

Nhận xét 10

Kết quả trên có thể áp dụng với vành số nguyên \(\mathbb{Z}\) nếu ta xem \(T\) là vành như sau:

\[T = \left\{\frac{m}{2^k} : m \in \mathbb{Z}, k \in \mathbb{Z}_{\geqslant 0}\right\}, \quad T \supseteq \mathbb{Z}.\]

Trên máy tính, mỗi phần tử của \(T\) được lưu chính xác, ví dụ dưới dạng cặp \((m, k)\).

Ta sẽ xem xét định lí quan trọng tiếp theo rồi quay lại chứng minh Bổ đề 4.

Định lí 7

Phép nhân hai đa thức có bậc không quá \(n\) trong vành \(T[x]\), với \(n \geqslant 3\), được tính với \(M(n) = O(n \log n)\) phép nhân trong \(T\) và \(A(n) = O(n \log n \log \log n)\) phép cộng trong \(T\).

Chứng minh

Gọi \(t \in \mathbb{N}\) là số sao cho \(2^{t-1} \leqslant 2n < 2^t\). Dễ thấy \(t \geqslant 2\) và \(2n < 2^t < 4n\).

Do đó phép nhân hai đa thức có bậc không quá \(n\) theo modulo \(x^{2^t} + 1\) là phép nhân thông thường, kết quả không thay đổi sau phép modulo.

Theo Bổ đề 4 thì ta có thể tính tích với \(M(n) = O(2^t \cdot t) = O(n \log n)\) phép nhân trong \(T\) và \(A(n) = O(2^t \cdot t \log t) = O(n \cdot \log n \cdot \log \log n)\) phép cộng trong \(T\).

Tiếp theo chúng ta quay lại chứng minh Bổ đề 4.

Chứng minh

Đặt \(F = F(x)\) và \(G = G(x)\) là các đa thức bậc không quá \(2^t - 1\). Đặt

\[H = H(x) = \sum_{i=0}^{2^t - 1} H_i x^i\]

với \(H = FG \pmod{x^{2^t} + 1}\). Ta cần tính các hệ số \(H_0\), \(H_1\), ..., \(H_{2^t - 1}\) theo các hệ số của \(F\) và \(G\).

Đặt \(k\) là tham số tự nhiên với \(1 \leqslant k < t\) mà ta sẽ chọn ở dưới.

Ta biểu diễn \(F\) và \(G\) dưới dạng

\[F = \sum_{i = 0}^{2^{t - k} - 1} f_i(x) x^{i \cdot 2^k}, \quad G = \sum_{i = 0}^{2^{t - k} - 1} g_i(x) x^{i \cdot 2^k}\]

với

• \(f_i(x), g_i(x) \in T[x]\);

• \(\deg f_i(x) \leqslant 2^k - 1\);

• \(\deg g_i(x) \leqslant 2^k - 1\).

Thuật toán tính \(F \cdot G\) được thực hiện theo các bước sau.

Bước 1. Nhân \(\sum\limits_{i=0}^{2^{t-k} - 1} f_i(x) Y_i\) với \(\sum\limits_{i=0}^{2^{t-k} - 1} g_i(x) Y_i\) trên vành \(T[x, Y] / (Y^{2^{t - k} - 1})\). Kết quả được biểu diễn bởi \(\widetilde{H} = \sum\limits_{i=0}^{2^{t-k} - 1} h_i(x) Y^i\).

Bước 2. Thay \(Y = x^{2^k}\) vào \(\widetilde{H}\) vào tính modulo \(x^{2^t} + 1 = Y^{2^{t - k}} + 1\). Khi đó

\[F(x) \cdot G(x) = \sum_{l=0}^{2^{t-k} - 1} f_l(x) x^{2^k l} \cdot \sum_{j=0}^{2^{t-k} - 1} g_j(x) x^{2^k j} \equiv \sum_{i=0}^{2^{t-k} - 1} h_i(x) x^{2^k i} \pmod{(x^{2^k})^{2^{t-k}} + 1}.\]

Từ bước 2 ta sẽ tính được \(H(x) = \sum\limits_{i=0}^{2^t - 1} H_i x^i\). Ta cần hiểu tại sao ở bước 2 dãy \(\{ h_i(x) \}\) theo modulo \(x^{2^t} + 1\) tìm được hệ số \(H_i\) của \(H(x)\).

Ở bước 1 ta nhân

\[\sum_{l=0}^{2^{t-k} - 1} f_l(x) Y^l \cdot \sum_{j=0}^{2^{t-k}-1} g_j(x) Y^j \equiv \sum_{i=0}^{2^{t-k} - 1} h_i(x) Y^i \pmod{Y^{2^{t-k}} + 1}.\]

Ở đây \(l + j \leqslant 2^{t - k + 1} - 2 < 2^{t-k+1} - 1\). Khi đó có hai trường hợp có thể xảy ra.

Trường hợp 1. Nếu \(0 \leqslant l + j \leqslant 2^{t-k} - 1\) thì

\[Y^l \cdot Y_j \equiv Y^{l + j} \pmod{Y^{2^{t - k}} + 1}.\]

Trường hợp 2. Nếu \(2^{t - k} \leqslant l + j \leqslant 2 \cdot 2^{t - k} - 2\) thì

\[Y^{l + j} \equiv -Y^i \pmod{Y^{2^{t-k}} + 1},\]

với \(i = l + j - 2^{t-k}\).

Khi đó

\[\begin{split}h_i(x) = \sum_{\substack{l, j \\ j + j = i}} f_l(x) g_j(x) - \sum_{\substack{l, j \\ l + j = i + 2^{t - k}}} f_l(x) g_j(x),\end{split}\]

với \(i = 0, \ldots, 2^{t - k} - 1\).

Từ đây suy ra

\[\deg h_i(x) \leqslant \deg f_l(x) = \deg g_j(x) \leqslant 2^{k+1} - 2 < 2^{k+1} - 1.\]

Vì \(k < t\) nên \(\deg h_i < 2^t - 1\).

Bước 2 thực hiện không nhiều hơn \(2^t\) phép cộng trên \(T\) vì nếu chúng ta đã biết các giá trị \(h_i(x) = \sum\limits_{j = 0}^{2^{k+1} - 1} h_{ij} x^j\) thì khi thay \(Y = x^{2^k}\) ta có biểu thức

\[\sum_{i = 0}^{2^{t-k} - 1} \sum_{j=0}^{2^{k+1} - 1} h_{ij} x^{j + 2^k i} \pmod{x^{2^t} + 1}.\]

Phép tính modulo \(x^{2^t} + 1\) ta sẽ được các biểu thức với \(j + 2^k i \geqslant 2^t\).

Khi đó, với \(j + 2^k i = r 2^t + l\), \(0 \leqslant l < 2^t\), nếu thay \(x^{j + 2^k i}\) thành \((-1)^r \cdot x^l\) thì đại lượng \((-1)^r \cdot h_{ij}\) sẽ được thêm vào hệ số của \(x^l\) (tức là thực hiện một phép cộng hoặc trừ). Các phép cộng như vậy sẽ không lớn hơn số lượng số hạng, tức là không lớn hơn \(2^{t - k} \cdot 2^{k + 1} = 2^{t + 1}\).

Ở đây, khi \(0 \leqslant i \leqslant 2^{t - k + 1}\) ta có bất đẳng thức

\[j + 2^k \cdot i \leqslant 2^{k + 1} - 1 + 2^{t - 1} - 2^k = 2^k + 2^{t - 1} - 1 \leqslant 2^t - 1\]

vì \(k \leqslant t - 1\). Với những số \(i\) như vậy thì phép tính modulo là không cần thiết.

Lúc này còn các giá trị \(i\) trong đoạn \(2^{t - k - 1} \leqslant i \leqslant 2^{t - k} - 1\), số lượng các giá trị này không nhiều hơn \(2^{t - k - 1}\). Khi đó, cặp \((i, j)\) ở bước 2 sẽ triệt tiêu theo modulo, và có không nhiều hơn \(2^{t - k - 1} \cdot 2^{k + 1} = 2^t\) cặp. Như vậy ta đã chứng minh được ở bước 2 thực hiện không nhiều hơn \(2^t\) phép cộng trong \(T\).

Bây giờ xét bước 1. Phép nhân ta không thực hiện trong vành \(T[x, Y] / (Y^{2^{t - k}} - 1)\) mà trong vành \(R[Y] / (Y^{2^{t - k}} - 1)\) với \(R = T[x] / (x^{2^{k+1}} + 1)\).

Trong vành \(R\) có phần tử \(\zeta_{2^{k+2}} \equiv x \pmod{x^{2^{k+1}} + 1}\) là nghiệm phương trình \(X^{2^{k+1}} + 1 = 0\).

Chọn \(k = \left[\dfrac{t}{2}\right] \geqslant 1\). Khi đó

\[k \geqslant \frac{t-1}{2}, \quad k \leqslant \frac{t}{2} < t.\]

Vì \(2^{t - k + 1} \leqslant 2^{k+2}\) nên trong \(R\) có phần tử \(\zeta_{2^{t-k+1}}\) là một lũy thừa của phần tử \(\zeta_{2^{k+2}}\).

Một phép nhân cho lũy thừa của phần tử \(\zeta_{2^{t-k+1}}\) trong vành \(R\) được thực hiện bởi \(2^{k+1}\) phép cộng trong \(T\).

Vì \(x \equiv \zeta_{2^{k+2}} \pmod{x^{2^{k+1}} + 1}\) nên

\[R = \left\{a_0 + a_1 \zeta_{2^{k+2}} + \cdots + a_{2^{k+1}} \zeta_{2^{k+2}}^{2^{k+1} - 1} : a_i \in T, i = 0, 1, \ldots, 2^{k+1} - 1\right\}.\]

Phép nhân với \(\zeta_{2^{t-k+1}}^j = \zeta_{2^{k+2}}^{j}\) với đẳng thức \(\zeta_{2^{k+2}}^{2^{k+1}} = -1\) sẽ cho kết quả là hoán vị các hệ số \(a_0\), \(a_1\), ..., \(a_{2^{k+1} - 1}\) và một số phần tử sẽ đổi dấu.

Ta dùng Bổ đề 3 để đánh giá độ phức tạp của một phép nhân trong vành \(R[Y] / (Y^{2^{t-k}} + 1)\).

Ta cần thực hiện:

• \(2^{t-k}\) phép nhân trong \(R\);

• \(3 \cdot (t - k) \cdot 2^{t-k}\) phép cộng trong \(R\) (ở đây có \(3 \cdot (t - k) \cdot 2^{t-k} \cdot 2^{k+1}\) phép cộng trong \(T\) vì các phần tử của vành \(R\) được biểu diễn bởi đa thực bậc không quá \(2^{k+1} - 1\) trong \(T[x]\));

• \(3 \cdot (t - k) \cdot 2^{t-k}\) phép nhân với lũy thừa của \(\zeta_{2^{t-k+1}}\);

• \(2^{t-k}\) phép nhân với \((2^{-1})^{t-k}\) (lần nữa, cần \(2^{t-k} \cdot 2^{k+1} = 2^{t+1}\) phép nhân trong \(T\) cho phần tử \((2^{-1})^{t-k}\)).

Tổng cộng ta cần:

• \(2^{t-k}\) phép nhân trong \(R\);

• \(6 \cdot (t - k) \cdot 2^{k+1}\) phép cộng trong \(T\);

• \(2^{t+1}\) phép nhân cho \((2^{-1})^{t-k}\) trong \(T\).

Ta tính thêm \(2^t\) phép cộng trong \(T\) ở bước 2, đặt

(6)\[k(t) = k + 1 = \left[\frac{t}{2}\right] + 1 \geqslant 2.\]

Khi đó, một phép nhân trong \(T[x] / (x^{2^{k+1}} + 1)\) được thực hiện bởi \(2^{t - k(t) + 1}\) phép nhân trong vành \(R = T[x] / (x^{2^{k(t)}} + 1)\), cộng với không nhiều hơn \(12 \cdot t \cdot 2^t\) phép cộng trong \(T\), cộng với \(2^{t+1}\) phép nhân trong \(T\).

Khi \(t \geqslant 3\) thì \(k(t) < t\). Ta chuyển phép nhân ở trên trong vành \(T[x] / (x^{2^t} + 1)\) thành phép nhân trong vành \(T[x] / (x^{2^{k(t)}} + 1)\), và đi xuống tiếp cho tới khi gặp vành \(T[x] / (x^4 + 1)\). Ở bước cuối này phép nhân bất kì đều cần \(O(1)\) phép cộng và phép nhân trong \(T\).

Ta kí hiệu \(M_1(2^j)\) và \(A_1(2^j)\) là số lượng phép nhân và cộng trong \(T\), tương ứng với số lượng cần thiết để thực hiện phép nhân ở trên trong vành \(T[x] / (x^{2^j} + 1)\).

Khi đó, với \(t \geqslant 3\) ta có bất đẳng thức

\[\begin{split}M_1(2^t) & \leqslant 2^{t - k(t) + 1} \cdot M_1(2^{k(t)}) + 2^{t+1}, \\ A_1(2^t) & \leqslant 2^{t - k(t) + 1} \cdot A_1(2^{k(t)}) + 12 \cdot t \cdot 2^t,\end{split}\]

với \(k(t)\) định nghĩa ở (6).

Đặt \(\alpha(j) = \dfrac{A_1(2^j)}{2^j}\) với \(j = 2, 3, \ldots\)

Khi đó \(\alpha(j) \leqslant 2\alpha(k(t)) + 12t\).

Giả sử ta có bất đẳng thức sau với \(2 \leqslant j \leqslant t\):

\[\alpha(j) \leqslant c \cdot j \log j\]

với hằng số tuyệt đối \(c\) nào đó.

Khi đó ta có bất đẳng thức

\[\alpha(t) \leqslant 2c \cdot k(t) \cdot \log(k(t)) + 12t \leqslant 2c \left(\frac{t}{2} + 1\right) \log\left(\frac{t}{2} + 1\right) + 12t < c \cdot t\log t\]

khi hằng số \(c\) đủ lớn.

Như vậy

\[A_1(2^t) \leqslant 2^t \cdot c \cdot t\log t = O(2^t \cdot t\log t)\]

và ta đã chứng minh được ý thứ hai của Bổ đề 4 về phép cộng.

Bây giờ, lại đặt \(\beta(j) = \dfrac{M_1(2^j)}{2^j}\) với \(j = 2, 3, \ldots\)

Ta có bất đẳng thức

\[\beta(t) \leqslant 2\beta(k(t)) + 2.\]

Từ đây suy ra

\[\beta(t) \leqslant 2(2\beta(k(k(t))) + 2) + 2 = 2^2\beta(k(k(t))) + 2^2 + 2,\]

và tương tự như vậy.

Vì \(k(t) \leqslant \dfrac{t}{2} + 1\) nên

\[k(k(t)) \leqslant \dfrac{1}{2}\left(\dfrac{t}{2} + 1\right) + 1 = \dfrac{t}{2^2} + 1 + \dfrac{1}{2}; \ldots; k(k(\ldots(k(t))\ldots)) < \dfrac{t}{2^j} + 2\]

với mọi \(j \geqslant 1\). Vì vậy khi \(j = \left[\log_2 t\right]\) ta có bất đẳng thức

\[\beta(t) \leqslant 2^j \cdot c_1 + 2 + 2^2 + \cdots + 2^j < 2^j \cdot c_1 + 2^{j+1} \leqslant c_2 t\]

với \(c_1\), \(c_2\) là các hằng số tuyệt đối.

Như vậy \(M_1(t) = O(t \cdot 2^t)\) và ta đã chứng minh xong ý đầu của Bổ đề 4.

Ví dụ nhân đa thức với Sympy

Để ví dụ, mình sẽ viết chương trình Python nhân hai đa thức bậc \(4\) trên \(\mathbb{Q}\) trong modulo \(x^4 + 1\).

Thuật toán ở các chứng minh trên là thuật toán Cooley-Tukey FFT.

Đầu tiên mình sẽ viết hàm đảo ngược bit (dùng cho radix-2 decimal-in-time).

def _reverse_bit(n: int, nbits: int) -> int:
  m = 0
  for _ in range(nbits):
      m = (m << 1) + (n & 1)
      n = n >> 1
  return m

assert _reverse_bit(3, 3) == 6
assert _reverse_bit(1, 3) == 4

FFT là thuật toán chia để trị -- để tính toán FFT cho vector độ dài \(n = 2^t\), ta sẽ tính trên hai vector độ dài \(n / 2\) và kết hợp kết quả lại.

from sympy import pi, sin, cos, I, simplify
from sympy.codegen.cfunctions import log2

def _mult_W(v: list[int]) -> list[int]:
  n = len(v)
  w = cos(2 * pi / n) + I * sin(2 * pi / n)
  u = [0] * n
  for i in range(0, n // 2):
      u[i] = v[i] + w**i * v[i + n // 2]
      u[i + n // 2] = v[i] + w**(i + n // 2) * v[i + n // 2]
  return u

def _fft(v: list[int]) -> list[int]:
    n = len(v)

    if n == 2:
        return _mult_W(v)
    else:
        u = _fft(v[:n // 2]) + _fft(v[n // 2:])
        return _mult_W(u)

def fft(v: list[int]) -> list[int]:
    n = len(v)
    b = int(log2(n))
    idx = [_reverse_bit(i, b) for i in range(n)]
    v = [v[idx[i]] for i in range(n)]

    return list(map(simplify, _fft(v)))

def _mult_iW(v: list[int]) -> list[int]:
    n = len(v)
    w = cos(-2 * pi / n) + I * sin(-2 * pi / n)
    u = [0] * n
    for i in range(0, n // 2):
        u[i] = v[i] + w**i * v[i + n // 2]
        u[i + n // 2] = v[i] + w**(i + n // 2) * v[i + n // 2]
    return u

def _ifft(v: list[int]) -> list[int]:
    n = len(v)
    if n == 2:
        return _mult_iW(v)
    else:
        u = _ifft(v[:n // 2]) + _ifft(v[n // 2:])
        return _mult_iW(u)

def ifft(v: list[int]) -> list[int]:
    n = len(v)
    b = int(log2(n))
    idx = [_reverse_bit(i, b) for i in range(n)]
    v = [v[idx[i]] for i in range(n)]

    return [simplify(i / n) for i in _ifft(v)]

Tiếp theo mình kiểm tra các hàm đã viết bên trên.

assert fft([2, 3, 0, 0]) == [5, 2 + 3*I, -1, 2 - 3*I]
assert ifft([5, 2 + 3*I, -1, 2 - 3*I]) == [2, 3, 0, 0]

import random
v = [random.randint(-4, 4) for _ in range(8)]
assert ifft(fft(v)) == v, fft(v)

Ví dụ, để nhân hai đa thức

\[\begin{split}a(x) & = 1 + 2x - x^2 + 3x^3, \\ b(x) & = -1 - 4x + 3x^2 - 2x^3,\end{split}\]

trong modulo \(x^4 + 1\), đầu tiên mình viết hệ số của mỗi đa thức thành vector theo số mũ tăng dần, như vậy

\[\bm{a} = (1, 2, -1, 3), \quad \bm{b} = (-1, -4, 3, -2).\]

Tiếp theo, mình sẽ thêm vào sau các vector \(\bm{a}\) và \(\bm{b}\) các số \(0\) để đạt độ dài là lũy thừa của \(2\) nhưng lớn hơn bậc của modulo. Ở đây modulo là \(x^4 + 1\) nên mình sẽ thêm các số \(0\) để đạt độ dài \(4 \cdot 2 = 8\).

\[\bm{a}' = (1, 2, -1, 3, 0, 0, 0, 0), \quad \bm{b}' = (-1, 4, 3, -2, 0, 0, 0, 0).\]

Tiếp theo, mình áp dụng FFT cho hai vector \(\bm{a}'\) và \(\bm{b}'\).

a = [1, 2, -1, 3, 0, 0, 0, 0]
b = [-1, -4, 3, -2, 0, 0, 0, 0]

fa = fft(a)
fb = fft(b)

Đặt

\[\begin{split}\bm{f}_a & = \mathsf{FFT}(\bm{a}') = (a_0, a_1, \ldots, a_7), \\ \bm{f}_b & = \mathsf{FFT}(\bm{b}') = (b_0, b_1, \ldots, b_7).\end{split}\]

Khi đó ta tính vector \(\bm{f}_c\) theo là tích theo cặp của \(\bm{f}_a\) và \(\bm{f}_b\), nghĩa là

\[\bm{f}_c = (c_0, c_1, \ldots, c_7), \ c_i = a_i \cdot b_i.\]

Tiếp theo, ta tính \(\bm{c}'\) là biến đổi Fourier rời rạc ngược của \(\bm{f}_c\):

\[\bm{c}' = \mathsf{FFT}^{-1}(\bm{f}_c) = (C_0, C_1, \ldots, C_7).\]

Khi đó, ta tính vector \(\bm{c}\) độ dài \(4\) theo công thức

\[\bm{c} = (C_i - C_{i+4})_{i = 0}^3.\]

fc = [i * j for i, j in zip(fa, fb)]
fc = list(map(simplify, ifft(fc)))
c = [fc[i] - fc[i + 4] for i in range(4)]
print(c)

Kết quả là vector \((18, -17, 2, 5)\).

Kiểm tra kết quả bằng tay, mình có

\[\begin{split}a(x) \cdot b(x) & = (1 + 2x - x^2 + 3x^3) \cdot (-1 - 4x + 3x^2 - 2x^3) \\ & = -1 - 4x + 3x^2 - 2x^3 \\ & - 2x - 8x^2 + 6x^3 - 4x^4 \\ & + x^2 + 4x^3 - 3x^4 + 2x^5 \\ & - 3x^3 - 12x^4 + 9x^5 - 6x^6 \\ & = -1 - 6x - 4x^2 + 5x^3 - 19x^4 + 11x^5 - 6x^6.\end{split}\]

Trong modulo \(x^4 + 1\) ta có thể thay thế \(x^4\) thành \(-1\), như vậy kết quả trên tương đương

\[\begin{split}a(x) \cdot b(x) & = -1 - 6x - 4x^2 + 5x^3 - 19x^4 + 11x^5 - 6x^6 \\ & = -1 - 6x - 4x^2 + 5x^3 - 19 \cdot (-1) + 11x \cdot (-1) - 6x^2 \cdot (-1) \\ & = 18 - 17x + 2x^2 + 5x^3 \bmod x^4 + 1.\end{split}\]

Nếu viết hệ số của kết quả dưới dạng vector theo số mũ tăng dần thì mình có \((18, -17, 2, 5)\), bằng đúng phép tính với thuật toán Cooley-Tukey.

Một lý do khiến mình viết code trong khi nhiều thư viện Python như numpy, sympy, sagemath đã hỗ trợ là vì trong [11] (tài liệu tham khảo chính của bài viết này) thì DFT thuận sử dụng \(\zeta_{2n}\) là nghiệm của \(x^n + 1\) nên nếu xét trên \(\mathbb{C}\) thì \(\zeta_{2n} = e^{i \pi / n}\). Trong khi đó các tài liệu khác như wikipedia, numpy, sympy, ... sử dụng \(e^{-i \pi / n}\) để tính DFT thuận. Ngược lại, [11] sử dụng \(e^{-i \pi /n }\) để tính DFT ngược, trong khi các tài liệu khác dùng \(e^{i \pi / n}\). Ở đây chúng ta có thể thấy dù là cách nào thì để thực hiện phép nhân nhanh đa thức modulo \(x^{2^t} + 1\) đều cần cả DFT thuận và ngược, nên kết quả không thay đổi nhờ vào tính chất của DFT.