Trong lĩnh vực đại số máy tính và toán học rời rạc, việc xử lý các phép toán trên đa thức với hiệu suất cao là rất quan trọng. Các thuật toán như biến đổi Fourier nhanh (FFT) và biến đổi số lý thuyết (NTT) là nền tảng cho nhiều phép toán đa thức phức tạp hơn. Bài viết này sẽ đi sâu vào một số kỹ thuật nâng cao sử dụng các công cụ này.
Phép Nhân Đa Thức Chia Để Trị (Divide and Conquer NTT)
Kỹ thuật này áp dụng ý tưởng của chia để trị, tương tự như CDQ divide and conquer, nhưng không yêu cầu kiến thức chuyên sâu về CDQ. Nguyên tắc cốt lõi là: ở mỗi bước, chúng ta chỉ xem xét đóng góp của nửa bên trái vào nửa bên phải. Đầu tiên, giải quyết đệ quy nửa bên trái, sau đó thêm đóng góp của nửa bên trái vào nửa bên phải, và cuối cùng, giải quyết đệ quy nửa bên phải. Bằng cách tích lũy các đóng góp này, chúng ta tính toán giá trị tại mỗi vị trí.
Một tình huống điển hình là khi \(P_i\) phụ thuộc vào tích chập của \(Q_j\) (với \(j \in [0, i]\)) và \(P_k\) (với \(k \in [0, i]\)). Kỹ thuật này yêu cầu một trong hai đa thức tham gia tích chập phải đã biết trước. Ví dụ, nếu chúng ta cần tính \(f_i = \sum_{j=0}^i g_j f_{i-j}\), trong đó \(g\) là đa thức đã biết, chúng ta có thể áp dụng chia để trị. Lưu ý rằng các chỉ số thường bắt đầu từ 0.
Các bước cơ bản trong quá trình chia để trị:
void process_polynomial_convolution(int left_idx, int right_idx) {
if (left_idx == right_idx) {
// Xử lý trường hợp cơ sở
return;
}
int mid_idx = (left_idx + right_idx) / 2;
process_polynomial_convolution(left_idx, mid_idx); // Giải quyết nửa trái
// Tính toán đóng góp từ nửa trái (P[left_idx...mid_idx]) đến nửa phải (P[mid_idx+1...right_idx])
// Đa thức A chứa các hệ số đã biết của Q.
// Đa thức B chứa các hệ số đã được tính của P từ nửa trái.
int current_len_q = mid_idx - left_idx + 1; // Số lượng hệ số của Q được dùng
for (int i = 0; i < current_len_q; ++i) {
temp_q_coeffs[i] = known_q_polynomial[i]; // Ví dụ: g[i]
}
int current_len_p = mid_idx - left_idx + 1; // Số lượng hệ số của P đã tính
for (int i = 0; i < current_len_p; ++i) {
temp_p_coeffs[i] = result_p_polynomial[left_idx + i]; // Ví dụ: f[L+j]
}
// Thực hiện phép nhân đa thức (ví dụ: bằng NTT)
// Sau khi nhân, cộng kết quả vào các hệ số tương ứng của P trong nửa phải.
// Ví dụ: result_p_polynomial[mid_idx + 1 + k] += (convolution_result[k])
process_polynomial_convolution(mid_idx + 1, right_idx); // Giải quyết nửa phải
}
Điều quan trọng là phải áp dụng chia để trị cho đa thức mà chúng ta đang tìm kiếm (ví dụ, \(f\)), để các đóng góp từ các đoạn nhỏ hơn đã được tính toán trước khi chúng được sử dụng.
Mở rộng: Tích chập giữa hai hàm phụ thuộc lẫn nhau
Đối với các bài toán có dạng \(f = \sum f \cdot g + \dots\) và \(g = \sum f \cdot g + \dots\), một phương pháp là giải bằng hàm sinh (generation function), nhưng nó có thể phức tạp. Một phương pháp khác là thực hiện chia để trị cho cả hai hàm cùng một lúc. Khi \(left\_idx = 0\), chúng ta tính đóng góp của \(f,g\) trong đoạn \(0 \dots mid\_idx\). Khi \(left\_idx > 0\), chúng ta phải tính cả \(f(left\_idx \dots mid\_idx) \times g(0 \dots len)\) và \(g(left\_idx \dots mid\_idx) \times f(0 \dots len)\), vì cả hai đều chưa được tính đóng góp.
NTT Với Modulo Bất Kỳ (MTT - Modulo Polynomial Multiplication)
Đôi khi, modulo \(P\) không có dạng \(a \cdot 2^k + 1\), khiến NTT tiêu chuẩn không thể áp dụng trực tiếp. Nếu \(P\) nhỏ (khoảng \(10^4\)), FFT số thực có thể được sử dụng. Tuy nhiên, với \(P\) lớn (khoảng \(10^9\)), chúng ta cần một kỹ thuật đặc biệt.
Một phương pháp là phân tách các hệ số đa thức thành dạng \(ax + b\), trong đó \(x\) là một lũy thừa của 2 đã chọn (thường là \(2^{15}\)). Khi đó, \(a, b\) sẽ không vượt quá \(2^{15}\) và có thể được xử lý riêng. Tích của hai số \((a_1x + b_1)(a_2x + b_2)\) sẽ là \(a_1a_2x^2 + (a_1b_2 + a_2b_1)x + b_1b_2\).
Để giảm hằng số, chúng ta có thể sử dụng các đa thức phức: \(P_c = (a, a')\), \(Q_c = (b, b')\), và \(R_c = (a, -a')\). Đặt \(F = P_c Q_c = (ab - a'b', ab' + a'b)\) và \(F' = R_c Q_c = (ab + a'b', ab' - a'b)\). Sau đó, ta có thể kết hợp chúng: \(F + F' = (2ab, 2ab')\) và \(F - F' = (-2a'b', 2a'b)\). Bằng cách này, chúng ta có thể thu được các thành phần cần thiết chỉ với 5 lần FFT.
Đoạn mã chính mô tả phương pháp này:
// Struct cho số phức hoặc cặp giá trị
struct ComplexPair { double real, imag; };
// Hàm FFT (hoặc NTT nếu có modulo phù hợp)
void perform_transform(ComplexPair* arr, int n, int type);
void multiply_polynomials_arbitrary_mod(int n_terms_a, int n_terms_b, long long modulus) {
// Giả sử n_terms_a và n_terms_b là số lượng hệ số (ví dụ: n+1 và m+1)
// Limi là kích thước biến đổi FFT/NTT, là lũy thừa của 2 lớn hơn tổng số hệ số
int total_degree = n_terms_a + n_terms_b - 2;
int transform_size = 1;
while (transform_size <= total_degree) transform_size <<= 1;
ComplexPair poly_A[transform_size], poly_B[transform_size], poly_A_conj[transform_size];
ComplexPair transformed_Q[transform_size]; // Để lưu trữ Q sau biến đổi
int split_point = 1 << 15; // Ví dụ: 2^15
long long mod_val_1 = 1LL << 30; // 2^30
long long mod_val_2 = 1LL << 15; // 2^15
for (int i = 0; i < n_terms_a; ++i) {
int val; // Đọc hệ số A[i]
// Ví dụ: read(val);
poly_A[i].real = val / split_point;
poly_A[i].imag = val % split_point;
poly_A_conj[i].real = val / split_point;
poly_A_conj[i].imag = -(val % split_point); // Liên hợp
}
for (int i = n_terms_a; i < transform_size; ++i) {
poly_A[i] = {0,0}; poly_A_conj[i] = {0,0};
}
for (int i = 0; i < n_terms_b; ++i) {
int val; // Đọc hệ số B[i]
// Ví dụ: read(val);
transformed_Q[i].real = val / split_point;
transformed_Q[i].imag = val % split_point;
}
for (int i = n_terms_b; i < transform_size; ++i) {
transformed_Q[i] = {0,0};
}
perform_transform(poly_A, transform_size, 1);
perform_transform(poly_A_conj, transform_size, 1);
perform_transform(transformed_Q, transform_size, 1);
ComplexPair result_P_Q[transform_size];
ComplexPair result_R_Q[transform_size];
for (int i = 0; i < transform_size; ++i) {
// Chia trước để tránh tràn số sau biến đổi ngược
// Đây là cách đơn giản, trong thực tế có thể cần chia sau hoặc xử lý khác
transformed_Q[i].real /= transform_size;
transformed_Q[i].imag /= transform_size;
// Tính P_c * Q_c và R_c * Q_c trong miền tần số
result_P_Q[i].real = poly_A[i].real * transformed_Q[i].real - poly_A[i].imag * transformed_Q[i].imag;
result_P_Q[i].imag = poly_A[i].real * transformed_Q[i].imag + poly_A[i].imag * transformed_Q[i].real;
result_R_Q[i].real = poly_A_conj[i].real * transformed_Q[i].real - poly_A_conj[i].imag * transformed_Q[i].imag;
result_R_Q[i].imag = poly_A_conj[i].real * transformed_Q[i].imag + poly_A_conj[i].imag * transformed_Q[i].real;
}
perform_transform(result_P_Q, transform_size, -1);
perform_transform(result_R_Q, transform_size, -1);
for (int i = 0; i <= total_degree; ++i) {
long long term_a1b1 = (long long)round((result_P_Q[i].real + result_R_Q[i].real) / 2.0);
long long term_a1b2 = (long long)round((result_P_Q[i].imag + result_R_Q[i].imag) / 2.0);
long long term_a2b1 = (long long)round((result_P_Q[i].imag - result_R_Q[i].imag) / 2.0);
long long term_a2b2 = (long long)round((result_R_Q[i].real - result_P_Q[i].real) / 2.0);
long long current_ans = (term_a1b1 % modulus * (mod_val_1 % modulus)) % modulus;
current_ans = (current_ans + (term_a2b1 + term_a1b2) % modulus * (mod_val_2 % modulus)) % modulus;
current_ans = (current_ans + term_a2b2 % modulus) % modulus;
current_ans = (current_ans + modulus) % modulus;
// In hoặc lưu current_ans
}
}
Nhân Đa Thức Lũy Thừa Giảm Dần
Một đa thức lũy thừa giảm dần có dạng: \[F(x) = \sum_{i} a_i x^{\underline{i}}\] trong đó \(x^{\underline{i}} = x(x-1)\dots(x-i+1)\) là lũy thừa giảm dần. Các hệ số của đa thức này có thể được chuyển đổi giữa hàm sinh thông thường (OGF) và hàm sinh mũ (EGF) của các giá trị điểm. Cụ thể, nếu \(f_i = \sum_j a_j i^{\underline{j}}\), thì:
\[\begin{aligned} F_{EGF}(x) &= \sum_{i} \frac{f_i}{i!} x^i \\ &= \sum_i \frac{x^i}{i!} \sum_j a_j i^{\underline{j}} \\ &= \sum_i x^i \sum_j a_j \frac{1}{(i-j)!} \\ &= \sum_j a_j x^j \sum_i \frac{x^i}{i!} \\ &= A_{OGF}(x) \cdot e^x \end{aligned}\]
Như vậy, để chuyển đổi từ các hệ số OGF của \(A(x)\) sang EGF của các giá trị điểm \(F(x)\), chúng ta nhân OGF với \(e^x\). Ngược lại, để chuyển đổi từ EGF của các giá trị điểm sang các hệ số OGF, chúng ta nhân với \(e^{-x}\). Lưu ý rằng \(e^{-x} = \sum_i \frac{(-1)^i}{i!} x^i\).
Vì các giá trị điểm có thể được nhân trong thời gian \(O(N)\), chúng ta có thể chuyển đổi đa thức lũy thừa giảm dần sang dạng điểm-giá trị, nhân chúng, rồi chuyển ngược lại. Độ phức tạp là \(O(N \log N)\).
// Giả sử C = A * B (trong miền lũy thừa giảm dần)
// A_coeffs, B_coeffs là các hệ số của A và B (OGF)
// factorial_inv[i] = 1/i! mod P
// factorial[i] = i! mod P
// 1. Chuyển đổi A và B sang EGF
// poly_exp_x_coeffs[i] = 1/i!
// poly_exp_neg_x_coeffs[i] = (i % 2 == 1 ? P - 1 : 1) * factorial_inv[i] % P;
// Tạo đa thức tạm thời cho e^x và e^-x
vector<long long> exp_x_poly(max_degree, 0);
vector<long long> exp_neg_x_poly(max_degree, 0);
for (int i = 0; i < max_degree; ++i) {
exp_x_poly[i] = factorial_inv[i];
exp_neg_x_poly[i] = (i % 2 == 1 ? MOD - 1 : 1) * factorial_inv[i] % MOD;
}
vector<long long> A_egf(max_degree), B_egf(max_degree);
// Hàm multiply_polynomials_ntt thực hiện nhân hai đa thức bằng NTT
multiply_polynomials_ntt(A_coeffs, exp_x_poly, A_egf, max_degree);
multiply_polynomials_ntt(B_coeffs, exp_x_poly, B_egf, max_degree);
// 2. Nhân các EGF điểm-giá trị (hệ số)
vector<long long> C_egf_product(max_degree);
for (int i = 0; i < max_degree; ++i) {
C_egf_product[i] = A_egf[i] * B_egf[i] % MOD;
// Chú ý: Ở đây có thể cần nhân với factorial[i] % MOD nếu các hệ số A_egf[i] và B_egf[i]
// đã được chia cho i! trong định nghĩa, để đưa về giá trị điểm.
// Nếu A_egf[i] = f_i/i!, thì C_egf_product[i] = (f_i/i!) * (g_i/i!) = (f_i g_i) / (i! * i!)
// Để có (f_i g_i)/i!, ta cần nhân với i!
C_egf_product[i] = C_egf_product[i] * factorial[i] % MOD;
}
// 3. Chuyển đổi C_egf_product trở lại OGF
vector<long long> C_coeffs(max_degree);
multiply_polynomials_ntt(C_egf_product, exp_neg_x_poly, C_coeffs, max_degree);
Các Thuật Toán Đa Thức "Công Nghiệp"
Phần này tập trung vào các phép toán đa thức cơ bản, thường được triển khai bằng phương pháp Newton hoặc đệ quy chia để trị.
Nghịch Đảo Đa Thức
Cho đa thức \(A(x)\), chúng ta muốn tìm \(B(x)\) sao cho \(A(x)B(x) \equiv 1 \pmod{x^N}\). Đây là một bài toán kinh điển giải quyết bằng phương pháp Newton.
Áp dụng lặp Newton cho hàm \(F(B) = AB - 1 = 0\):
\[B_{k+1} = B_k - \frac{F(B_k)}{F'(B_k)} = B_k - \frac{A B_k - 1}{A} = B_k - B_k(A B_k - 1) = B_k(2 - A B_k)\]
Quá trình này đệ quy: để tìm \(B(x) \pmod{x^N}\), chúng ta trước hết tìm \(B(x) \pmod{x^{\lceil N/2 \rceil}}\).
long long poly_coeffs_a[MAX_N], poly_inv_b[MAX_N];
long long temp_c_coeffs[MAX_N]; // Đa thức tạm
int current_ntt_size;
int current_bit_length; // log2(current_ntt_size)
int bit_reverse_indices[MAX_N]; // Mảng đảo bit
void compute_bit_reverse_indices(int size) {
// Tính toán mảng đảo bit cho NTT
if (current_ntt_size != size) { // Cập nhật khi kích thước thay đổi
current_ntt_size = size;
current_bit_length = 0;
while ((1 << current_bit_length) < size) current_bit_length++;
for (int i = 0; i < size; ++i) {
bit_reverse_indices[i] = (bit_reverse_indices[i >> 1] >> 1) | ((i & 1) << (current_bit_length - 1));
}
}
}
// Hàm NTT (viết tắt)
void ntt_transform(long long* arr, int size, int type);
// type = 1 cho forward, type = -1 cho inverse, đã tích hợp inverse của size
void polynomial_inverse(int degree, long long* input_a, long long* output_b) {
if (degree == 1) {
output_b[0] = power(input_a[0], MOD - 2); // Nghịch đảo modulo của hệ số hằng
return;
}
// Đệ quy để tìm nghịch đảo cho nửa độ
polynomial_inverse((degree + 1) / 2, input_a, output_b);
// Tính toán kích thước NTT cho bước hiện tại
int transform_size = 1;
while (transform_size < degree * 2) transform_size <<= 1;
compute_bit_reverse_indices(transform_size);
// Sao chép input_a vào temp_c_coeffs và điền 0
for (int i = 0; i < degree; ++i) temp_c_coeffs[i] = input_a[i];
for (int i = degree; i < transform_size; ++i) temp_c_coeffs[i] = 0;
// Biến đổi NTT
ntt_transform(output_b, transform_size, 1);
ntt_transform(temp_c_coeffs, transform_size, 1);
// Áp dụng công thức Newton
for (int i = 0; i < transform_size; ++i) {
output_b[i] = (output_b[i] * (2 - temp_c_coeffs[i] * output_b[i] % MOD + MOD)) % MOD;
}
// Biến đổi NTT ngược
ntt_transform(output_b, transform_size, -1);
// Đặt 0 cho các hệ số ngoài phạm vi degree
for (int i = degree; i < transform_size; ++i) {
output_b[i] = 0;
}
}
Một phiên bản lặp (iterative) có thể hiệu quả hơn:
void polynomial_inverse_iterative(long long* input_a, long long* output_b, int degree) {
output_b[0] = power(input_a[0], MOD - 2); // Khởi tạo b[0]
// len là độ hiện tại của nghịch đảo chúng ta đang tính (len * 2 sẽ là độ tiếp theo)
for (int len = 1; len < degree * 2; len <<= 1) {
int transform_size = len << 1;
compute_bit_reverse_indices(transform_size); // Cập nhật mảng đảo bit
// Sao chép các hệ số cần thiết
for (int i = 0; i < len; ++i) {
temp_c_coeffs[i] = input_a[i];
// temp_d_coeffs[i] = output_b[i]; // Nếu output_b bị sửa đổi trong ntt_transform
}
for (int i = len; i < transform_size; ++i) {
temp_c_coeffs[i] = 0;
// temp_d_coeffs[i] = 0;
}
ntt_transform(temp_c_coeffs, transform_size, 1);
ntt_transform(output_b, transform_size, 1); // Sử dụng output_b trực tiếp hoặc temp_d_coeffs
for (int i = 0; i < transform_size; ++i) {
output_b[i] = (output_b[i] * (2 - temp_c_coeffs[i] * output_b[i] % MOD + MOD)) % MOD;
}
ntt_transform(output_b, transform_size, -1);
// Đặt 0 cho các hệ số không cần thiết (từ độ 'len' trở đi)
for (int i = len; i < transform_size; ++i) {
output_b[i] = 0;
}
}
}
Căn Bậc Hai Đa Thức
Cho đa thức \(A(x)\) với \(A(0) = 1\), tìm \(B(x)\) sao cho \(B(x)^2 \equiv A(x) \pmod{x^N}\). Tương tự như nghịch đảo, ta sử dụng phương pháp Newton. Đặt \(F(B) = B^2 - A = 0\). Lặp Newton cho ra:
\[B_{k+1} = B_k - \frac{B_k^2 - A}{2B_k} = \frac{B_k^2 + A}{2B_k} = \frac{1}{2} \left(B_k + \frac{A}{B_k}\right)\]
Điều kiện \(A(0) = 1\) đảm bảo \(B(0) = 1\) (hoặc \(-1\)), giúp xác định điều kiện cơ sở đệ quy. Nếu \(A(0)\) không phải là 1, ta cần sử dụng các kỹ thuật số dư bậc hai.
long long poly_a[MAX_N], poly_sqrt_b[MAX_N];
long long temp_inv_b[MAX_N]; // Để lưu nghịch đảo của B
long long temp_prod_ac[MAX_N]; // Để lưu A * (1/B)
const long long INV2 = power(2, MOD - 2); // (MOD + 1) / 2
void polynomial_sqrt(int degree, long long* input_a, long long* output_b) {
if (degree == 1) {
output_b[0] = 1; // Vì A[0] = 1, B[0]^2 = 1, chọn B[0] = 1
return;
}
polynomial_sqrt((degree + 1) / 2, input_a, output_b); // Đệ quy
int transform_size = 1;
while (transform_size < degree * 2) transform_size <<= 1;
compute_bit_reverse_indices(transform_size);
// Tính nghịch đảo của output_b
polynomial_inverse_iterative(output_b, temp_inv_b, degree);
// Sao chép input_a vào temp_prod_ac
for (int i = 0; i < degree; ++i) temp_prod_ac[i] = input_a[i];
for (int i = degree; i < transform_size; ++i) temp_prod_ac[i] = 0;
// Nhân input_a với nghịch đảo của output_b (A * (1/B))
ntt_transform(temp_prod_ac, transform_size, 1);
ntt_transform(temp_inv_b, transform_size, 1);
for (int i = 0; i < transform_size; ++i) {
temp_prod_ac[i] = temp_prod_ac[i] * temp_inv_b[i] % MOD;
}
ntt_transform(temp_prod_ac, transform_size, -1);
// Áp dụng công thức B = 1/2 * (B + A/B)
for (int i = 0; i < degree; ++i) {
output_b[i] = (output_b[i] + temp_prod_ac[i]) % MOD * INV2 % MOD;
}
for (int i = degree; i < transform_size; ++i) { // Xóa các hệ số không cần thiết
output_b[i] = 0;
}
}
Lưu ý: Luôn nhớ khởi tạo hoặc xóa sạch các mảng tạm thời trước khi sử dụng để tránh lỗi từ các phép toán NTT trước đó.
Phép Chia Đa Thức
Cho hai đa thức \(D(x)\) và \(A(x)\), tìm đa thức thương \(B(x)\) và đa thức dư \(C(x)\) sao cho \(D(x) = A(x)B(x) + C(x)\), với \(\deg(C) < \deg(A)\). Số bậc của \(D(x)\) là \(N\), và \(A(x)\) là \(M\) (\(M \le N\)).
Chúng ta định nghĩa phép đảo ngược đa thức \(F^R(x) = x^k F(1/x)\), trong đó \(k = \deg(F)\). Phép toán này đảo ngược thứ tự các hệ số.
Áp dụng phép đảo ngược cho phương trình: \(D(x) = A(x)B(x) + C(x)\). Giả sử \(\deg(D) = n-1\), \(\deg(A) = m-1\), \(\deg(B) = n-m\), \(\deg(C) = m-2\).
Thay \(x\) bằng \(1/x\) và nhân với \(x^{n-1}\) (bậc cao nhất):
\[D^R(x) = A^R(x) B^R(x) + C^R(x) x^{n-m+1}\]
Lấy modulo \(x^{n-m+1}\):
\[D^R(x) \equiv A^R(x) B^R(x) \pmod{x^{n-m+1}}\]
Từ đây, chúng ta có thể tìm \(B^R(x) \pmod{x^{n-m+1}}\) bằng cách nhân \(D^R(x)\) với nghịch đảo của \(A^R(x)\) (tính modulo \(x^{n-m+1}\)). Sau khi có \(B^R(x)\), ta đảo ngược nó để có \(B(x)\). Cuối cùng, tính \(C(x) = D(x) - A(x)B(x)\).
long long poly_D[MAX_N], poly_A[MAX_N];
long long poly_D_orig[MAX_N], poly_A_orig[MAX_N]; // Sao lưu bản gốc
long long poly_A_rev_inv[MAX_N]; // Nghịch đảo của A_rev
long long poly_B[MAX_N], poly_C[MAX_N]; // Thương và số dư
long long poly_B_rev_temp[MAX_N]; // B_rev tạm thời
long long poly_AB_product[MAX_N]; // A * B
// Hàm đảo ngược đa thức
void reverse_polynomial(long long* poly, int degree) {
for (int i = 0; i < degree / 2; ++i) {
swap(poly[i], poly[degree - 1 - i]);
}
}
int main() {
int N_degree, M_degree; // Bậc của D và A
// Ví dụ: read(N_degree); read(M_degree);
for (int i = 0; i <= N_degree; ++i) { // Lưu D và sao lưu
// Ví dụ: read(poly_D[i]);
poly_D_orig[i] = poly_D[i];
}
for (int i = 0; i <= M_degree; ++i) { // Lưu A và sao lưu
// Ví dụ: read(poly_A[i]);
poly_A_orig[i] = poly_A[i];
}
int N_terms = N_degree + 1; // Số lượng hệ số
int M_terms = M_degree + 1;
// Đảo ngược A và D
reverse_polynomial(poly_A, M_terms);
reverse_polynomial(poly_D, N_terms);
int B_terms_count = N_terms - M_terms + 1; // Số hệ số của B
// Tính nghịch đảo của A_rev modulo x^(B_terms_count)
polynomial_inverse_iterative(poly_A, poly_A_rev_inv, B_terms_count);
// Nhân D_rev với (A_rev)^-1 để tìm B_rev
// Thực hiện nhân đa thức bằng NTT
multiply_polynomials_ntt(poly_D, poly_A_rev_inv, poly_B_rev_temp, B_terms_count, B_terms_count, B_terms_count);
// Đảo ngược B_rev để có B
reverse_polynomial(poly_B_rev_temp, B_terms_count);
for (int i = 0; i < B_terms_count; ++i) {
poly_B[i] = poly_B_rev_temp[i];
}
// Tính C(x) = D(x) - A(x)B(x)
// Nhân A_orig và B để có A*B
multiply_polynomials_ntt(poly_A_orig, poly_B, poly_AB_product, M_terms, B_terms_count, N_terms);
for (int i = 0; i < M_terms - 1; ++i) { // Độ của C < độ của A
poly_C[i] = (poly_D_orig[i] - poly_AB_product[i] + MOD) % MOD;
}
// In B và C
// ...
return 0;
}
Logarit Đa Thức (ln)
Cho đa thức \(A(x)\) với \(A(0) = 1\), tìm \(\ln(A(x)) \pmod{x^N}\). Ta sử dụng tính chất \(\frac{d}{dx} \ln(A(x)) = \frac{A'(x)}{A(x)}\). Vậy, \(\ln(A(x)) = \int \frac{A'(x)}{A(x)} dx\). Các bước thực hiện:
- Tính đạo hàm của \(A(x)\): \(A'(x)\).
- Tính nghịch đảo của \(A(x)\): \(A^{-1}(x)\).
- Nhân \(A'(x)\) với \(A^{-1}(x)\).
- Lấy tích phân của kết quả để có \(\ln(A(x))\).
long long poly_input_A[MAX_N], poly_output_ln_B[MAX_N];
long long temp_A_prime[MAX_N]; // Đạo hàm của A
long long temp_A_inv[MAX_N]; // Nghịch đảo của A
long long temp_product_poly[MAX_N]; // Tích A' * A^-1
// Hàm tính đạo hàm đa thức
void polynomial_derivative(long long* input_poly, long long* output_poly, int degree) {
for (int i = 0; i < degree - 1; ++i) {
output_poly[i] = input_poly[i+1] * (i+1) % MOD;
}
output_poly[degree - 1] = 0; // Bậc cao nhất sau đạo hàm sẽ là 0
}
// Hàm tính tích phân đa thức
void polynomial_integral(long long* input_poly, long long* output_poly, int degree) {
output_poly[0] = 0; // Hệ số hằng của ln(A(x)) là 0 vì A(0)=1
for (int i = 0; i < degree - 1; ++i) {
output_poly[i+1] = input_poly[i] * power(i+1, MOD - 2) % MOD;
}
}
void polynomial_ln(long long* input_a, long long* output_b, int degree) {
// 1. Tính A'(x)
polynomial_derivative(input_a, temp_A_prime, degree);
// 2. Tính A^-1(x)
polynomial_inverse_iterative(input_a, temp_A_inv, degree);
// 3. Nhân A'(x) * A^-1(x)
int transform_size = 1;
while (transform_size < degree * 2) transform_size <<= 1;
compute_bit_reverse_indices(transform_size);
ntt_transform(temp_A_prime, transform_size, 1);
ntt_transform(temp_A_inv, transform_size, 1);
for (int i = 0; i < transform_size; ++i) {
temp_product_poly[i] = temp_A_prime[i] * temp_A_inv[i] % MOD;
}
ntt_transform(temp_product_poly, transform_size, -1);
// Đảm bảo các hệ số ngoài phạm vi degree bị xóa
for (int i = degree -1; i < transform_size; ++i) temp_product_poly[i] = 0;
// 4. Lấy tích phân
polynomial_integral(temp_product_poly, output_b, degree);
}
Mũ Đa Thức (exp)
Cho đa thức \(A(x)\) với \(A(0) = 0\), tìm \(B(x)\) sao cho \(B(x) = e^{A(x)} \pmod{x^N}\). Áp dụng phương pháp Newton cho \(F(B) = \ln(B) - A = 0\).
Lặp Newton: \[B_{k+1} = B_k - \frac{\ln(B_k) - A}{\frac{1}{B_k}} = B_k(1 - \ln(B_k) + A)\]
Điều kiện \(A(0) = 0\) đảm bảo \(B(0) = e^0 = 1\), là điều kiện cơ sở cho đệ quy. Biểu thức \(1 - \ln(B_k) + A\) được thực hiện bằng cách cộng/trừ các hệ số tương ứng. Lưu ý rằng "1" chỉ ảnh hưởng đến hệ số hằng số (hệ số bậc 0).
Lũy Thừa Đa Thức Nhanh (Polynomial Power)
Cho đa thức \(A(x)\) và số nguyên \(k\), tính \(A(x)^k \pmod{x^N}\). Phương pháp đơn giản là lũy thừa nhị phân (binary exponentiation) với phép nhân đa thức, cho độ phức tạp \(O(N \log N \log k)\). Sử dụng \(\ln\) và \(\exp\), ta có thể đạt \(O(N \log N)\):
\[A(x)^k = e^{k \ln(A(x))}\]
Các bước thực hiện:
- Nếu \(A(0) = 0\) và \(k > 0\): Kết quả là \(0\) nếu \(N > 1\) (trừ khi \(k\) đủ nhỏ để \(A(x)^k\) có bậc 0 là \(0\) và các bậc cao hơn cũng là \(0\)).
- Nếu \(A(0) \ne 0\):
- Trích hệ số hằng số \(c = A(0)\). Đặt \(A_0(x) = A(x)/c\), sao cho \(A_0(0) = 1\).
- Tính \(\ln(A_0(x))\).
- Nhân \(k\) với \(\ln(A_0(x))\).
- Tính \(e^{k \ln(A_0(x))}\).
- Nhân kết quả với \(c^k \pmod P\).
Đoạn mã dưới đây là một khung tổng thể tích hợp các hàm \(\text{NTT}\), \(\text{inv}\), \(\text{dao}\) (đạo hàm), \(\text{ji}\) (tích phân), \(\ln\), \(\exp\), và sử dụng chúng để tính \((e^x + 1)^n\).
namespace PolynomialAdvancedOperations {
const int PRIME_MOD = 998244353;
const int GENERATOR = 3; // Nguyên hàm căn bậc nguyên thủy
long long power(long long base, long long exp) {
long long res = 1;
base %= PRIME_MOD;
while (exp > 0) {
if (exp % 2 == 1) res = (res * base) % PRIME_MOD;
base = (base * base) % PRIME_MOD;
exp /= 2;
}
return res;
}
long long INV_GENERATOR = power(GENERATOR, PRIME_MOD - 2);
long long INV_2 = power(2, PRIME_MOD - 2);
int bit_rev_arr[MAX_N];
int current_bit_len;
void ntt_transform(long long* arr, int size, int type) {
if ((1 << current_bit_len) != size) {
current_bit_len = 0;
while ((1 << current_bit_len) < size) current_bit_len++;
for (int i = 0; i < size; ++i) {
bit_rev_arr[i] = (bit_rev_arr[i >> 1] >> 1) | ((i & 1) << (current_bit_len - 1));
}
}
for (int i = 0; i < size; ++i) {
if (i < bit_rev_arr[i]) swap(arr[i], arr[bit_rev_arr[i]]);
}
for (int len = 1; len < size; len <<= 1) {
long long root = power(type == 1 ? GENERATOR : INV_GENERATOR, (PRIME_MOD - 1) / (len << 1));
for (int i = 0; i < size; i += (len << 1)) {
long long w = 1;
for (int j = 0; j < len; ++j) {
long long u = arr[i + j];
long long v = (arr[i + len + j] * w) % PRIME_MOD;
arr[i + j] = (u + v) % PRIME_MOD;
arr[i + len + j] = (u - v + PRIME_MOD) % PRIME_MOD;
w = (w * root) % PRIME_MOD;
}
}
}
if (type == -1) {
long long inv_size = power(size, PRIME_MOD - 2);
for (int i = 0; i < size; ++i) {
arr[i] = (arr[i] * inv_size) % PRIME_MOD;
}
}
}
long long temp_inv_A[MAX_N], temp_inv_B[MAX_N];
void get_polynomial_inverse(long long* poly_A, long long* poly_B_inv, int degree) {
poly_B_inv[0] = power(poly_A[0], PRIME_MOD - 2);
for (int current_len = 1; (current_len >> 1) < degree; current_len <<= 1) {
int transform_size = current_len << 1;
for (int i = 0; i < current_len; ++i) {
temp_inv_A[i] = poly_A[i];
temp_inv_B[i] = poly_B_inv[i];
}
for (int i = current_len; i < transform_size; ++i) {
temp_inv_A[i] = temp_inv_B[i] = 0;
}
ntt_transform(temp_inv_A, transform_size, 1);
ntt_transform(temp_inv_B, transform_size, 1);
for (int i = 0; i < transform_size; ++i) {
poly_B_inv[i] = (temp_inv_B[i] * (2 - temp_inv_A[i] * temp_inv_B[i] % PRIME_MOD + PRIME_MOD)) % PRIME_MOD;
}
ntt_transform(poly_B_inv, transform_size, -1);
for (int i = current_len; i < transform_size; ++i) poly_B_inv[i] = 0;
}
}
void polynomial_derivative(long long* input_poly, int degree) {
for (int i = 0; i < degree - 1; ++i) input_poly[i] = (input_poly[i+1] * (i+1)) % PRIME_MOD;
input_poly[degree - 1] = 0;
}
// (fact_inv[x] = 1/x! mod P)
void polynomial_integral(long long* input_poly, int degree, const vector<long long>& fact_inv) {
for (int i = degree - 1; i > 0; --i) input_poly[i] = (input_poly[i-1] * power(i, PRIME_MOD - 2)) % PRIME_MOD; // or use fact_inv[i] * factorial[i-1]
input_poly[0] = 0;
}
long long ln_temp_A[MAX_N], ln_temp_B[MAX_N];
void get_polynomial_ln(long long* poly_A_in, long long* poly_B_ln, int degree, const vector<long long>& fact_inv) {
for (int i = 0; i < degree; ++i) ln_temp_A[i] = poly_A_in[i];
get_polynomial_inverse(ln_temp_A, ln_temp_B, degree); // B_inv = 1/A
polynomial_derivative(ln_temp_A, degree); // A'
int transform_size = 1; while(transform_size < degree * 2) transform_size <<= 1;
ntt_transform(ln_temp_A, transform_size, 1); // A' transformed
ntt_transform(ln_temp_B, transform_size, 1); // A_inv transformed
for (int i = 0; i < transform_size; ++i) ln_temp_A[i] = (ln_temp_A[i] * ln_temp_B[i]) % PRIME_MOD; // A' * A_inv
ntt_transform(ln_temp_A, transform_size, -1);
for (int i = degree -1; i < transform_size; ++i) ln_temp_A[i] = 0;
polynomial_integral(ln_temp_A, degree, fact_inv); // Integral of A' * A_inv
for (int i = 0; i < degree; ++i) poly_B_ln[i] = ln_temp_A[i];
for (int i = 0; i < transform_size; ++i) ln_temp_A[i] = ln_temp_B[i] = 0; // Clear temps
}
long long exp_temp_A[MAX_N], exp_temp_B[MAX_N];
void get_polynomial_exp(long long* poly_A_in, long long* poly_B_exp, int degree, const vector<long long>& fact_inv) {
poly_B_exp[0] = 1; // B[0] = e^(A[0]) = e^0 = 1
for (int current_len = 1; (current_len >> 1) < degree; current_len <<= 1) {
int transform_size = current_len << 1;
for (int i = 0; i < current_len; ++i) exp_temp_A[i] = poly_B_exp[i]; // Copy B_exp for ln
for (int i = current_len; i < transform_size; ++i) exp_temp_A[i] = 0;
get_polynomial_ln(exp_temp_A, exp_temp_B, current_len, fact_inv); // exp_temp_B = ln(B_exp)
// B_new = B_exp * (1 - ln(B_exp) + A)
// Handle (1 + A - ln(B_exp)) directly on A_in
exp_temp_A[0] = (poly_A_in[0] + 1 - exp_temp_B[0] + PRIME_MOD) % PRIME_MOD;
for (int i = 1; i < current_len; ++i) exp_temp_A[i] = (poly_A_in[i] - exp_temp_B[i] + PRIME_MOD) % PRIME_MOD;
for (int i = current_len; i < transform_size; ++i) exp_temp_A[i] = 0;
// Copy current B_exp to exp_temp_B for multiplication
for (int i = 0; i < current_len; ++i) exp_temp_B[i] = poly_B_exp[i];
for (int i = current_len; i < transform_size; ++i) exp_temp_B[i] = 0;
ntt_transform(exp_temp_A, transform_size, 1);
ntt_transform(exp_temp_B, transform_size, 1);
for (int i = 0; i < transform_size; ++i) exp_temp_A[i] = (exp_temp_A[i] * exp_temp_B[i]) % PRIME_MOD;
ntt_transform(exp_temp_A, transform_size, -1);
for (int i = 0; i < transform_size; ++i) poly_B_exp[i] = exp_temp_A[i];
for (int i = current_len; i < transform_size; ++i) poly_B_exp[i] = 0; // Clear unused parts
}
}
// Main polynomial power function for A(x)^k
void polynomial_power(long long* poly_A_in, long long* poly_res, int degree_N, long long exponent_K, const vector<long long>& fact_inv) {
if (exponent_K == 0) {
poly_res[0] = 1;
for (int i = 1; i < degree_N; ++i) poly_res[i] = 0;
return;
}
int first_nonzero = 0;
while (first_nonzero < degree_N && poly_A_in[first_nonzero] == 0) {
first_nonzero++;
}
if (first_nonzero == degree_N) { // A(x) is zero polynomial
if (exponent_K > 0) {
for (int i = 0; i < degree_N; ++i) poly_res[i] = 0;
} else { // 0^0 is usually 1, 0^negative undefined
poly_res[0] = 1; // Convention for 0^0
for (int i = 1; i < degree_N; ++i) poly_res[i] = 0;
}
return;
}
if (first_nonzero * exponent_K >= degree_N) { // Result will be 0 up to degree_N
for (int i = 0; i < degree_N; ++i) poly_res[i] = 0;
return;
}
long long c_const = poly_A_in[first_nonzero];
long long inv_c = power(c_const, PRIME_MOD - 2);
long long c_to_k = power(c_const, exponent_K);
long long temp_poly_normalized[MAX_N];
for (int i = 0; i < degree_N - first_nonzero; ++i) {
temp_poly_normalized[i] = (poly_A_in[i + first_nonzero] * inv_c) % PRIME_MOD;
}
for (int i = degree_N - first_nonzero; i < degree_N; ++i) temp_poly_normalized[i] = 0;
long long ln_result[MAX_N];
get_polynomial_ln(temp_poly_normalized, ln_result, degree_N - first_nonzero, fact_inv);
for (int i = 0; i < degree_N - first_nonzero; ++i) {
ln_result[i] = (ln_result[i] * exponent_K) % PRIME_MOD;
}
for (int i = degree_N - first_nonzero; i < degree_N; ++i) ln_result[i] = 0;
long long exp_result[MAX_N];
get_polynomial_exp(ln_result, exp_result, degree_N - first_nonzero, fact_inv);
for (int i = 0; i < degree_N; ++i) poly_res[i] = 0; // Clear result poly
for (int i = 0; i < degree_N - first_nonzero; ++i) {
if (i + first_nonzero * exponent_K < degree_N) {
poly_res[i + first_nonzero * exponent_K] = (exp_result[i] * c_to_k) % PRIME_MOD;
}
}
}
}
Các Thuật Toán Đa Thức "Brute Force"
Khi không thể sử dụng FFT/NTT (ví dụ, do giới hạn về bộ nhớ hoặc các ràng buộc khác), có thể dùng các thuật toán đa thức \(O(N^2)\) bằng cách quan sát hệ số thứ \(m\) trong một quan hệ. Các phương pháp này thường dựa trên công thức truy hồi.
Nghịch Đảo
Cho \(A(x)B(x) = 1\). Xét hệ số thứ \(m\):
\[\sum_{i=0}^m a_i b_{m-i} = [m=0]\] \[a_0 b_m + \sum_{i=1}^m a_i b_{m-i} = [m=0]\] Khi \(m > 0\), \(a_0 b_m = -\sum_{i=1}^m a_i b_{m-i}\). Từ đó, \(b_m = -\frac{1}{a_0}\sum_{i=1}^m a_i b_{m-i}\). Với \(b_0 = \frac{1}{a_0}\), chúng ta có thể tính \(b_m\) theo các \(b_i\) đã biết. Độ phức tạp là \(O(N^2)\).
Logarit/Mũ
Đối với \(\ln A = B \implies A'/A = B' \implies A' = B'A\). Xét hệ số thứ \(m\) của \(A'\) và \(B'A\):
\[(m+1)a_{m+1} = \sum_{i=0}^m (i+1)b_{i+1} a_{m-i}\] Từ đây, có thể suy ra công thức truy hồi cho \(b_{m+1}\). Độ phức tạp \(O(N^2)\).
Lũy Thừa Nhanh
Đối với \(A^k = B \implies k \ln A = \ln B \implies k A'/A = B'/B \implies k A'B = B'A\). Xét hệ số thứ \(m\):
\[k \sum_{i=0}^m b_i a'_{m-i} = \sum_{i=0}^m (i+1) b_{i+1} a_{m-i}\] Tương tự, có thể suy ra công thức truy hồi cho \(b_{m+1}\). Độ phức tạp \(O(N^2)\).