Biến đổi Fourier Nhanh (FFT) là một thuật toán tối ưu hóa phép nhân đa thức, giúp giảm độ phức tạp từ \(O(N^2)\) xuống \(O(N \log N)\). Để hiểu rõ về FFT, chúng ta cần tìm hiểu các khái niệm cơ bản về đa thức và số phức.

Đa thức

Định nghĩa

Một đa thức là biểu thức toán học có dạng tổng của các đơn thức, mỗi đơn thức gồm một hệ số và một biến số được nâng lên lũy thừa nguyên không âm. Cụ thể, một đa thức \(f(x)\) có thể được biểu diễn như sau:

\[f(x) = \sum_{i=0}^{n} a_ix^i = a_nx^n + a_{n-1}x^{n-1} + \dots + a_1x + a_0\]

Trong đó, \(a_i\) là các hệ số, và \(x\) là biến số. Nếu \(a_n \neq 0\), thì \(n\) được gọi là bậc của đa thức, ký hiệu là \(\text{deg}(f) = n\).

Biểu diễn đa thức và các phép toán cơ bản

1. Biểu diễn bằng hệ số

Một đa thức \(f(x) = \sum_{i=0}^{n} a_ix^i\) có thể được biểu diễn đơn giản bằng một mảng (hoặc vector) các hệ số của nó: \(\mathbf{a} = [a_0, a_1, \dots, a_n]\).

Cộng đa thức

Để cộng hai đa thức \(A(x) = \sum_{i=0}^n a_ix^i\) và \(B(x) = \sum_{i=0}^n b_ix^i\) (giả sử chúng có cùng bậc \(n\) bằng cách thêm các hệ số 0 nếu cần), ta cộng các hệ số tương ứng:

\[C(x) = A(x) + B(x) = \sum_{i=0}^n (a_i + b_i)x^i\]

Độ phức tạp: \(O(n)\).

Nhân đa thức

Phép nhân hai đa thức \(A(x)\) (bậc \(n\)) và \(B(x)\) (bậc \(m\)) cho ra một đa thức \(C(x)\) có bậc \(n+m\). Các hệ số của \(C(x)\) được tính bằng tích chập của các hệ số của \(A(x)\) và \(B(x)\):

\[C(x) = A(x)B(x) = \sum_{k=0}^{n+m} c_kx^k\]

Trong đó, \(c_k = \sum_{j=0}^k a_j b_{k-j}\).

Độ phức tạp: \(O(nm)\), thường là \(O(n^2)\) nếu \(n \approx m\).

Tính giá trị đa thức tại một điểm

Để tính \(f(x_0)\), ta có thể sử dụng thuật toán Horner (sơ đồ Horner) để giảm số phép nhân:

\[\begin{aligned} f(x) &= a_nx^n + a_{n-1}x^{n-1} + \dots + a_1x + a_0 \\ &= (\dots((a_nx + a_{n-1})x + a_{n-2})x + \dots + a_1)x + a_0 \end{aligned}\]

Độ phức tạp: \(O(n)\).

2. Biểu diễn bằng tập hợp điểm-giá trị

Một đa thức bậc \(n\) có thể được xác định duy nhất bởi \(n+1\) cặp điểm-giá trị \(\{(x_0, A(x_0)), (x_1, A(x_1)), \dots, (x_n, A(x_n))\}\), với các \(x_i\) phân biệt.

Cộng đa thức

Để cộng hai đa thức \(A(x)\) và \(B(x)\) khi chúng được biểu diễn bằng tập hợp điểm-giá trị tại cùng một tập hợp các điểm \(x_i\):

\[C(x_i) = A(x_i) + B(x_i)\]

Độ phức tạp: \(O(n)\).

Nhân đa thức

Tương tự như phép cộng, phép nhân đa thức trong biểu diễn điểm-giá trị cũng rất hiệu quả:

\[C(x_i) = A(x_i)B(x_i)\]

Độ phức tạp: \(O(n)\).

Khôi phục đa thức từ tập hợp điểm-giá trị (Nội suy Lagrange)

Để tìm lại các hệ số của đa thức từ \(n+1\) cặp điểm-giá trị, ta có thể dùng công thức nội suy Lagrange:

\[A(x) = \sum_{j=0}^{n} A(x_j) \prod_{\substack{k=0 \\ k \neq j}}^{n} \frac{x - x_k}{x_j - x_k}\]

Độ phức tạp: \(O(n^2)\) để tính giá trị tại một điểm \(x\), và \(O(n^3)\) để tìm tất cả các hệ số.

Ví dụ cài đặt nội suy Lagrange:

double interpolate_lagrange(const std::vector<double>& points_x, const std::vector<double>& values_y, double target_x) {
    int n = points_x.size();
    double result = 0.0;
    for (int i = 0; i < n; ++i) {
        double term = values_y[i];
        for (int j = 0; j < n; ++j) {
            if (i == j) continue;
            term *= (target_x - points_x[j]) / (points_x[i] - points_x[j]);
        }
        result += term;
    }
    return result;
}

So sánh hai cách biểu diễn, ta thấy phép nhân trong biểu diễn điểm-giá trị có độ phức tạp thấp hơn nhiều. Vấn đề là làm thế nào để chuyển đổi giữa hai biểu diễn một cách hiệu quả, đặc biệt là chuyển từ hệ số sang điểm-giá trị và ngược lại. Đây chính là nơi FFT phát huy tác dụng.

Số phức

Các phép toán cơ bản

Số phức có dạng \(a + bi\), trong đó \(a, b\) là số thực và \(i\) là đơn vị ảo (\(i^2 = -1\)). Các phép toán cơ bản:

• Cộng: \((a+bi) + (c+di) = (a+c) + (b+d)i\)
• Trừ: \((a+bi) - (c+di) = (a-c) + (b-d)i\)
• Nhân: \((a+bi)(c+di) = (ac-bd) + (ad+bc)i\)
• Chia: \(\frac{a+bi}{c+di} = \frac{(a+bi)(c-di)}{(c+di)(c-di)} = \frac{ac+bd}{c^2+d^2} + \frac{bc-ad}{c^2+d^2}i\)

Công thức Euler

Công thức Euler thiết lập mối liên hệ sâu sắc giữa hàm mũ và các hàm lượng giác:

\[e^{i\theta} = \cos \theta + i\sin \theta\]

Công thức này cho phép biểu diễn số phức dưới dạng tọa độ cực trên mặt phẳng phức, rất hữu ích trong các phép quay và biến đổi.

Định lý cơ bản của Đại số

Một định lý quan trọng trong đại số phát biểu rằng mọi đa thức bậc \(n\) với hệ số phức (và \(n \ge 1\)) đều có đúng \(n\) nghiệm phức (kể cả bội). Điều này đảm bảo rằng phương trình \(x^n = 1\) sẽ có đúng \(n\) nghiệm phức.

Căn bậc n của đơn vị (Unit Roots)

Các nghiệm của phương trình \(x^n = 1\) trong mặt phẳng phức được gọi là các căn bậc \(n\) của đơn vị. Chúng có dạng \(\omega_n^k = e^{i\frac{2\pi k}{n}}\) với \(k = 0, 1, \dots, n-1\).

Các căn bậc \(n\) của đơn vị có một số tính chất quan trọng được sử dụng trong FFT:

Tính chất đối xứng (Folding property)

Nếu \(n\) là số chẵn:

\[\omega_{2n}^{k+n} = -\omega_{2n}^k\]

Chứng minh:

\[\begin{aligned} \omega_{2n}^{k+n} &= e^{i\frac{2\pi(k+n)}{2n}} \\ &= e^{i(\frac{2\pi k}{2n} + \pi)} \\ &= e^{i\frac{2\pi k}{2n}} \cdot e^{i\pi} \\ &= \omega_{2n}^k \cdot (-1) \\ &= -\omega_{2n}^k \end{aligned}\]

Tính chất giảm bậc (Reduction property)

Nếu \(d\) là số nguyên dương:

\[\omega_{dn}^{dk} = \omega_n^k\]

Chứng minh:

\[\begin{aligned} \omega_{dn}^{dk} &= e^{i\frac{2\pi (dk)}{dn}} \\ &= e^{i\frac{2\pi k}{n}} \\ &= \omega_n^k \end{aligned}\]

Tổng các căn đơn vị

Tổng của tất cả các căn bậc \(n\) của đơn vị bằng 0 (trừ trường hợp \(k \equiv 0 \pmod n\)):

\[\sum_{j=0}^{n-1} (\omega_n^k)^j = \begin{cases} n & \text{if } k \equiv 0 \pmod n \\ 0 & \text{if } k \not\equiv 0 \pmod n \end{cases}\]

Đối với trường hợp \(k \not\equiv 0 \pmod n\), ta có:

\[\sum_{j=0}^{n-1} (\omega_n^k)^j = \frac{(\omega_n^k)^n - 1}{\omega_n^k - 1} = \frac{(\omega_n^n)^k - 1}{\omega_n^k - 1} = \frac{1^k - 1}{\omega_n^k - 1} = 0\]

Biến đổi Fourier Nhanh (FFT)

Mục tiêu của FFT là chuyển đổi một đa thức từ biểu diễn hệ số sang biểu diễn điểm-giá trị một cách hiệu quả. Như đã phân tích, phép nhân đa thức sẽ nhanh hơn rất nhiều ở dạng điểm-giá trị. Cụ thể, FFT thực hiện phép biến đổi Fourier rời rạc (DFT) trong \(O(N \log N)\) thời gian.

Giả sử chúng ta có đa thức \(A(x) = \sum_{i=0}^{N-1} a_ix^i\). Để đơn giản hóa, ta giả sử \(N\) là lũy thừa của 2 (nếu không, ta thêm các hệ số 0 vào cuối).

FFT sử dụng chiến lược chia để trị. Ta chia đa thức \(A(x)\) thành hai đa thức con: một chứa các hệ số của các bậc chẵn và một chứa các hệ số của các bậc lẻ.

\[\begin{aligned} A(x) &= a_0 + a_1x + a_2x^2 + \dots + a_{N-1}x^{N-1} \\ &= (a_0 + a_2x^2 + \dots + a_{N-2}x^{N-2}) + (a_1x + a_3x^3 + \dots + a_{N-1}x^{N-1}) \\ &= (a_0 + a_2(x^2) + \dots + a_{N-2}(x^2)^{(N-2)/2}) + x(a_1 + a_3(x^2) + \dots + a_{N-1}(x^2)^{(N-2)/2}) \end{aligned}\]

Đặt \(A_{even}(y) = a_0 + a_2y + \dots + a_{N-2}y^{(N-2)/2}\) và \(A_{odd}(y) = a_1 + a_3y + \dots + a_{N-1}y^{(N-2)/2}\). Khi đó:

\[A(x) = A_{even}(x^2) + x \cdot A_{odd}(x^2)\]

Mỗi đa thức \(A_{even}\) và \(A_{odd}\) đều có bậc \(N/2 - 1\).

Để chuyển đổi sang biểu diễn điểm-giá trị, ta cần tính \(A(\omega_N^k)\) cho \(k = 0, \dots, N-1\). Sử dụng các tính chất của căn đơn vị, ta có:

• Với \(k < N/2\):

  \[\begin{aligned} A(\omega_N^k) &= A_{even}((\omega_N^k)^2) + \omega_N^k \cdot A_{odd}((\omega_N^k)^2) \\ &= A_{even}(\omega_N^{2k}) + \omega_N^k \cdot A_{odd}(\omega_N^{2k}) \\ &= A_{even}(\omega_{N/2}^k) + \omega_N^k \cdot A_{odd}(\omega_{N/2}^k) \end{aligned}\]

• Với \(k \ge N/2\) (đặt \(k' = k - N/2\), sao cho \(k' < N/2\)):

  \[\begin{aligned} A(\omega_N^k) &= A(\omega_N^{k'+N/2}) \\ &= A_{even}((\omega_N^{k'+N/2})^2) + \omega_N^{k'+N/2} \cdot A_{odd}((\omega_N^{k'+N/2})^2) \\ &= A_{even}(\omega_N^{2k'+N}) + \omega_N^{k'+N/2} \cdot A_{odd}(\omega_N^{2k'+N}) \\ &= A_{even}(\omega_N^{2k'}) + \omega_N^{k'+N/2} \cdot A_{odd}(\omega_N^{2k'}) \\ &= A_{even}(\omega_{N/2}^{k'}) - \omega_N^{k'} \cdot A_{odd}(\omega_{N/2}^{k'}) \quad (\text{do } \omega_N^{k'+N/2} = -\omega_N^{k'}) \end{aligned}\]

Hai công thức này cho phép ta tính \(A(\omega_N^k)\) cho tất cả các \(k\) bằng cách đệ quy tính \(A_{even}\) và \(A_{odd}\). Độ phức tạp của thuật toán đệ quy này là \(T(N) = 2T(N/2) + O(N)\), dẫn đến \(T(N) = O(N \log N)\).

Biến đổi Fourier Nhanh Ngược (Inverse FFT - IFFT)

Sau khi nhân hai đa thức ở dạng điểm-giá trị, ta cần chuyển kết quả trở lại dạng hệ số. Đây là lúc IFFT phát huy tác dụng.

Biến đổi Fourier rời rạc (DFT) có thể được biểu diễn dưới dạng nhân ma trận: \(\mathbf{y} = \mathbf{F} \mathbf{a}\), trong đó \(\mathbf{y}\) là vector giá trị đầu ra, \(\mathbf{a}\) là vector hệ số, và \(\mathbf{F}\) là ma trận Fourier với phần tử \(F_{jk} = \omega_N^{jk}\).

\[\begin{bmatrix}y_0 \\ y_1 \\ y_2 \\ \vdots \\ y_{N-1} \end{bmatrix}=\begin{bmatrix}1 & 1 & 1 & \cdots & 1 \\1 & \omega_N^1 & \omega_N^2 & \cdots & \omega_N^{N-1} \\1 & \omega_N^2 & \omega_N^4 & \cdots & \omega_N^{2(N-1)} \\\vdots & \vdots & \vdots & \ddots & \vdots \\1 & \omega_N^{N-1} & \omega_N^{2(N-1)} & \cdots & \omega_N^{(N-1)^2} \end{bmatrix}\begin{bmatrix} a_0 \\ a_1 \\ a_2 \\ \vdots \\ a_{N-1} \end{bmatrix}\]

Để tìm lại các hệ số \(\mathbf{a}\) từ \(\mathbf{y}\), ta cần tính \(\mathbf{a} = \mathbf{F}^{-1} \mathbf{y}\). May mắn thay, ma trận nghịch đảo của \(\mathbf{F}\) có dạng rất đẹp: \(F^{-1}_{jk} = \frac{1}{N}\omega_N^{-jk}\).

Điều này có nghĩa là IFFT về cơ bản giống như FFT, chỉ khác ở chỗ ta dùng \(\omega_N^{-1} = e^{-i\frac{2\pi}{N}}\) thay vì \(\omega_N^1\), và chia kết quả cuối cùng cho \(N\). Điều này cho phép chúng ta dùng cùng một hàm đệ quy FFT, chỉ cần truyền một tham số \(direction\) (chẳng hạn, \(1\) cho FFT và \(-1\) cho IFFT) để điều khiển dấu của \(\sin\) trong căn đơn vị.

Cài đặt đệ quy của FFT/IFFT:

#include <vector>
#include <complex>
#include <cmath>
#include <cstdio> // For scanf/printf

const double PI_CONST = acos(-1.0);

// Sử dụng std::complex để đơn giản hóa việc xử lý số phức
using Complex = std::complex<double>;

// Hàm thực hiện Biến đổi Fourier Nhanh (FFT) hoặc Biến đổi Fourier Nhanh Ngược (IFFT)
// poly_coeffs: mảng các hệ số (hoặc giá trị điểm) của đa thức
// transform_size: kích thước của phép biến đổi (phải là lũy thừa của 2)
// direction: 1 cho FFT (biến đổi thuận), -1 cho IFFT (biến đổi ngược)
void fft_recursive(Complex poly_coeffs[], int transform_size, int direction) {
    if (transform_size == 1) {
        return;
    }

    Complex poly_even[transform_size / 2];
    Complex poly_odd[transform_size / 2];

    for (int i = 0; i < transform_size / 2; ++i) {
        poly_even[i] = poly_coeffs[2 * i];
        poly_odd[i] = poly_coeffs[2 * i + 1];
    }

    fft_recursive(poly_even, transform_size / 2, direction);
    fft_recursive(poly_odd, transform_size / 2, direction);

    Complex root_of_unity_N_k = Complex(1, 0); // omega_N^k
    // Gốc đơn vị cơ bản: e^(i * 2 * pi / N) hoặc e^(-i * 2 * pi / N)
    Complex base_root_of_unity = Complex(cos(2.0 * PI_CONST / transform_size), direction * sin(2.0 * PI_CONST / transform_size));

    for (int k = 0; k < transform_size / 2; ++k) {
        // poly_coeffs[k] = A_even(omega_{N/2}^k) + omega_N^k * A_odd(omega_{N/2}^k)
        poly_coeffs[k] = poly_even[k] + root_of_unity_N_k * poly_odd[k];
        // poly_coeffs[k + transform_size/2] = A_even(omega_{N/2}^k) - omega_N^k * A_odd(omega_{N/2}^k)
        poly_coeffs[k + transform_size / 2] = poly_even[k] - root_of_unity_N_k * poly_odd[k];
        root_of_unity_N_k *= base_root_of_unity;
    }
}

int main() {
    int poly_degree_A, poly_degree_B;
    scanf("%d %d", &poly_degree_A, &poly_degree_B);

    int N_transform_size = 1;
    while (N_transform_size <= poly_degree_A + poly_degree_B) {
        N_transform_size <<= 1;
    }

    // Khởi tạo mảng số phức với kích thước N_transform_size
    // Gán phần thực từ các hệ số đa thức, phần ảo bằng 0
    std::vector<Complex> coeffs_A(N_transform_size);
    std::vector<Complex> coeffs_B(N_transform_size);

    for (int i = 0; i <= poly_degree_A; ++i) {
        double real_part;
        scanf("%lf", &real_part);
        coeffs_A[i] = Complex(real_part, 0);
    }
    for (int i = 0; i <= poly_degree_B; ++i) {
        double real_part;
        scanf("%lf", &real_part);
        coeffs_B[i] = Complex(real_part, 0);
    }

    // Thực hiện FFT cho cả hai đa thức
    fft_recursive(coeffs_A.data(), N_transform_size, 1);
    fft_recursive(coeffs_B.data(), N_transform_size, 1);

    // Nhân các giá trị điểm
    for (int i = 0; i < N_transform_size; ++i) {
        coeffs_A[i] *= coeffs_B[i];
    }

    // Thực hiện IFFT để chuyển về hệ số
    fft_recursive(coeffs_A.data(), N_transform_size, -1);

    // In kết quả
    for (int i = 0; i <= poly_degree_A + poly_degree_B; ++i) {
        // Chia cho N_transform_size (tính chất của IFFT)
        // Làm tròn đến số nguyên gần nhất
        printf("%d ", (int)(coeffs_A[i].real() / N_transform_size + 0.5));
    }
    printf("\n");

    return 0;
}

Các tối ưu hóa cho FFT lặp

Phiên bản FFT đệ quy có thể gặp vấn đề về overhead khi stack quá sâu. Ta có thể chuyển sang phiên bản lặp để cải thiện hiệu suất.

1. Hoán vị đảo bit (Bit-Reversal Permutation)

Để thực hiện FFT lặp, ta cần sắp xếp lại các hệ số đầu vào sao cho khi các phép toán "bướm" được thực hiện ở các tầng khác nhau, chúng luôn thao tác trên đúng các phần tử cần thiết. Thứ tự sắp xếp này chính là hoán vị đảo bit. Nếu ta có một mảng kích thước \(N\), phần tử ở chỉ số \(i\) sẽ được chuyển đến vị trí có chỉ số là số \(i\) sau khi đảo ngược các bit nhị phân của nó.

Ví dụ với \(N=8\) (3 bit):

• \(000_2 \to 000_2\) (0)
• \(001_2 \to 100_2\) (4)
• \(010_2 \to 010_2\) (2)
• \(011_2 \to 110_2\) (6)
• \(100_2 \to 001_2\) (1)
• \(101_2 \to 101_2\) (5)
• \(110_2 \to 011_2\) (3)
• \(111_2 \to 111_2\) (7)

Công thức tính chỉ số đảo bit \(\text{bit\_reversed\_indices}[i]\) từ \(\text{bit\_reversed\_indices}[i \gg 1]\) có thể là: \(\text{bit\_reversed\_indices}[i] = (\text{bit\_reversed\_indices}[i \gg 1] \gg 1) | ((i \& 1) \ll (\text{log}_2 N - 1))\).

2. Phép toán cánh bướm (Butterfly Operation)

Ở mỗi bước kết hợp trong FFT, ta tính hai giá trị \(V_1\) và \(V_2\) dựa trên các giá trị đã được tính toán ở bước trước. Cụ thể, trong công thức đệ quy:

\[\begin{aligned} V_1 &= A_{even}(\omega_{N/2}^k) + \omega_N^k \cdot A_{odd}(\omega_{N/2}^k) \\ V_2 &= A_{even}(\omega_{N/2}^k) - \omega_N^k \cdot A_{odd}(\omega_{N/2}^k) \end{aligned}\]

Ta thấy \(\omega_N^k \cdot A_{odd}(\omega_{N/2}^k)\) là một thành phần chung. Bằng cách lưu trữ giá trị này vào một biến tạm thời (\(term\_t\)), ta có thể giảm số phép nhân phức từ 2 xuống 1 cho mỗi cặp phép toán. Phép toán này, được gọi là phép toán cánh bướm, là nền tảng của thuật toán FFT lặp.

Cài đặt FFT/IFFT lặp với hoán vị đảo bit và phép toán cánh bướm:

#include <vector>
#include <complex>
#include <cmath>
#include <cstdio> // For scanf/printf
#include <algorithm> // For std::swap

const double PI_CONST = acos(-1.0);

using Complex = std::complex<double>;

// Mảng lưu trữ các chỉ số đảo bit
std::vector<int> bit_reversed_indices;

// Hàm tính toán các chỉ số đảo bit cho một kích thước biến đổi N
void precompute_bit_reversal(int N_size, int log_N) {
    bit_reversed_indices.resize(N_size);
    for (int i = 0; i < N_size; ++i) {
        bit_reversed_indices[i] = (bit_reversed_indices[i >> 1] >> 1) | ((i & 1) << (log_N - 1));
    }
}

// Hàm thực hiện FFT/IFFT lặp
void fft_iterative(std::vector<Complex>& poly_coeffs, int N_size, int direction) {
    // Sắp xếp lại các phần tử theo thứ tự đảo bit
    for (int i = 0; i < N_size; ++i) {
        if (i < bit_reversed_indices[i]) {
            std::swap(poly_coeffs[i], poly_coeffs[bit_reversed_indices[i]]);
        }
    }

    // Các tầng biến đổi (từ 1 đến log_N)
    for (int stage = 1; stage <= std::log2(N_size); ++stage) {
        int block_size = 1 << stage; // Kích thước của các khối đang xử lý
        // Gốc đơn vị cho tầng hiện tại (e^(i*2*pi/block_size) hoặc nghịch đảo của nó)
        Complex base_root_of_unity = Complex(cos(2.0 * PI_CONST / block_size), direction * sin(2.0 * PI_CONST / block_size));

        // Lặp qua từng khối
        for (int k_block = 0; k_block < N_size; k_block += block_size) {
            Complex current_twiddle_factor = Complex(1, 0); // omega_block_size^j
            // Lặp qua nửa đầu của mỗi khối
            for (int j = 0; j < block_size / 2; ++j) {
                Complex term_t = current_twiddle_factor * poly_coeffs[k_block + j + block_size / 2];
                Complex term_u = poly_coeffs[k_block + j];

                poly_coeffs[k_block + j] = term_u + term_t;
                poly_coeffs[k_block + j + block_size / 2] = term_u - term_t;
                current_twiddle_factor *= base_root_of_unity;
            }
        }
    }

    // Nếu là IFFT, chia kết quả cho N_size
    if (direction == -1) {
        for (int i = 0; i < N_size; ++i) {
            poly_coeffs[i] /= N_size;
        }
    }
}

int main() {
    int n, m;
    scanf("%d %d", &n, &m);

    int N_size = 1;
    int log_N = 0;
    while (N_size <= n + m) {
        N_size <<= 1;
        log_N++;
    }

    precompute_bit_reversal(N_size, log_N);

    std::vector<Complex> poly_A(N_size);
    std::vector<Complex> poly_B(N_size);

    for (int i = 0; i <= n; ++i) {
        double val;
        scanf("%lf", &val);
        poly_A[i] = Complex(val, 0);
    }
    for (int i = 0; i <= m; ++i) {
        double val;
        scanf("%lf", &val);
        poly_B[i] = Complex(val, 0);
    }

    fft_iterative(poly_A, N_size, 1); // FFT
    fft_iterative(poly_B, N_size, 1); // FFT

    for (int i = 0; i < N_size; ++i) {
        poly_A[i] *= poly_B[i]; // Element-wise multiplication in frequency domain
    }

    fft_iterative(poly_A, N_size, -1); // IFFT

    for (int i = 0; i <= n + m; ++i) {
        printf("%d ", (int)(poly_A[i].real() + 0.5)); // Round to nearest integer
    }
    printf("\n");

    return 0;
}