Để có trải nghiệm tốt nhất khi đọc bài viết này, bạn đọc nên nắm vững kiến thức về:
Ảnh lấy từ báo Danviet
Bài viết sẽ mở đầu bằng bài toán "Hàn Tín điểm binh", được chép lại trong quyển hạ sách "Tôn Tử toán kinh". Đây cũng là lý do cách giải bài toán lại được gọi là "Định lý Thặng dư Trung Hoa".
Hàn Tín là một danh tướng nổi tiếng của Trung Quốc thời Hán Sở tranh hùng. Tên tuổi ông gắn bó với nhiều câu chuyện thú vị, trong đó có một bài toán với tên gọi "Hàn Tín điểm binh".
Tương truyền, để kiểm đếm số binh sĩ dưới quyền, Hàn Tín đã nghĩ ra một cách. Ông cho quân lính xếp hàng ba (mỗi hàng ba người), hàng năm, hàng bảy rồi ghi các số lẻ tương ứng và suy ra số lính.
Giả sử số người lẻ ra khi xếp hàng ba, hàng năm, hàng bảy lần lượt là , , . Gọi số lính cần tìm là . Khi đó, sử dụng phép toán modulo hiện đại, ta biểu diễn được các điều kiện của :
Lời giải đơn giản nhất là ... mò các bội của
Trở về lập trình thi đấu, đôi khi ta sẽ gặp những bài toán kiểu như trên, nhưng với số điều kiện và dữ liệu đầu vào khủng khiếp hơn nhiều. Lúc này, ta không thể "mò" như cách làm trên, và cũng không thể trông đợi vào một công thức "trên trời rơi xuống" để tìm ra nghiệm. Sử dụng các công cụ hiện đại, ta sẽ có cách để giải bài toán này.
Đề bài: LQDOJ - Hàn Tín điểm binh
Giải hệ gồm
trong đó
Dễ thấy nếu
Để thuận tiện, trong toàn bộ bài viết này tác giả sẽ sử dụng các ký hiệu sau:
Ta có thể xây dựng nghiệm của phương trình theo hướng sau: Tạo ra một tổng gồm nhiều số hạng, sao cho các số hạng này thoả mãn cả hai điều kiện:
Nói cách khác, nghiệm
Xét phương trình
Với ý tưởng xây dựng nghiệm như đã nói ở trên,
Để
Sau khi nhân xong, sẽ không có gì xảy ra nếu
Do
Bằng cách làm hoàn toàn tương tự như trên cho các phương trình khác trong hệ, ta sẽ xây dựng các số hạng
Nhắc lại, nếu một số
Đây chính là công thức nghiệm được đưa ra bởi định lý Thặng dư Trung Hoa.
Ngoài đưa ra công thức nghiệm, định lý thặng dư Trung hoa cũng khẳng định công thức nghiệm trên là duy nhất:
Định lý Thặng dư Trung Hoa (Chinese Remainder Theorem, CRT): Hệ phương trình (*) có họ nghiệm duy nhất là:
Chứng minh sự tồn tại: Do
Xét
Vậy
Chứng minh sự duy nhất: Giả sử
Do các số
Vậy (*) có nghiệm duy nhất theo modulo
Với định lý Thặng dư Trung Hoa trong tay, để giải hệ phương trình đồng dư ta chỉ cần tính các giá trị
Lưu ý rằng, không phải lúc nào
const int MAXN = 8;
const pair <long long, long long>
INVALID_ROOT = {LLONG_MAX, LLONG_MAX};
// CTDL biểu diễn đồng dư thức, gồm số dư và modulo
struct Congruence {
long long rem;
long long mod;
Congruence(long long __rem, long long __mod) {
rem = __rem;
mod = __mod;
if (mod) {
rem %= mod;
if (rem < 0) {
rem += mod;
}
}
}
};
// Hàm trả về ƯCLN của a và b, biến đổi x, y thoả mãn ax + by = gcd(a, b)
long long extEuclid(long long a, long long b,
long long &x, long long& y) {
if (b == 0) {
x = 1;
y = 0;
return a;
}
long long q = a / b;
long long r = a - b * q;
long long x1 = 0, y1 = 0;
long long d = extEuclid(b, r, x1, y1);
x = y1;
y = x1 - q * y1;
return d;
}
// Tính nghịch đảo modulo bằng cách giải phương trình Diophantus
long long invMod(long long num, long long mod) {
long long inv = 0;
long long y = 0;
extEuclid(num, mod, inv, y);
inv %= mod;
if (inv < 0) {
inv += mod;
}
return inv;
}
// Nghiệm của phương trình là một đồng dư thức mà, đúng không?
Congruence solveCongruenceEqSet(vector <Congruence>& eqSet) {
Congruence sol(0, 1);
for (int i = 0; i < (int)eqSet.size(); i++) {
sol.mod *= eqSet[i].mod;
}
vector <long long> p;
vector <long long> pInv;
for (int i = 0; i < (int)eqSet.size(); i++) {
p.push_back(sol.mod / eqSet[i].mod);
pInv.push_back(invMod(p[i], eqSet[i].mod));
}
for (int i = 0; i < (int)eqSet.size(); i++) {
long long prod = 1;
prod = prod * p[i] % sol.mod;
prod = prod * pInv[i] % sol.mod;
prod = prod * eqSet[i].rem % sol.mod;
sol.rem += prod;
}
sol.rem %= sol.mod;
return sol;
}
Với các tham số đã nêu trên, thì độ phức tạp không gian của thuật toán là
Về thời gian, với mỗi phương trình ta mất
Tất cả nhưng gì ta vừa làm đều chỉ đúng với các modulo nguyên tố cùng nhau đôi một. Khi các modulo không còn nguyên tố cùng nhau, giá trị
Ta đã biết đồng dư thức có hai tính chất sau với các số nguyên
Do đó, ta có thể tách một phương trình thành nhiều phương trình sao cho tích các modulo của chúng bằng đúng modulo của phương trình ban đầu. Cụ thể, giả sử
Thực hiện tương tự với các phương trình còn lại, ta sẽ thu được một hệ các phương trình có modulo là luỹ thừa của các số nguyên tố. Lúc này, việc kiểm tra tính mâu thuẫn của các điều kiện là tương đối dễ dàng: nếu có
Ngược lại, ta có thể kết luận ngay là hệ vô nghiệm.
Nếu không có cặp phương trình nào xung đột với nhau, hệ phương trình chắc chắn có nghiệm. Lúc này, dễ thấy các phương trình có modulo thuộc cùng một cơ số, chẳng hạn
Chẳng hạn, xét hệ phương trình:
Ta phân tích các modulo ra được:
Xét phương trình thứ nhất và phương trình thứ ba, với
Đoạn chương trình dưới đây cài đặt các ý tưởng trên, sử dụng phép phân tích bằng sàng nguyên tố và lưu lại ước nhỏ nhất của các số.
const int MAXX = 1e6 + 1;
int minPDiv[MAXX];
// Sàng nguyên tố (nhớ gọi hàm này trước khi làm bất cứ việc gì)
void sieve() {
for (int i = 2; i < MAXX; i ++) {
if (minPDiv[i]) {
continue;
}
minPDiv[i] = i;
for (int j = i + i; j < MAXX; j += i) {
minPDiv[j] = (minPDiv[j] ? minPDiv[j] : i);
}
}
}
// Phân tích một số ra TSNT dưới dạng tích của các luỹ thừa
vector<int> poweredPrimeDivisors(int x) {
vector<int> ret;
int pp = 1;
int curPDiv = minPDiv[x];
while (x != 1) {
if (minPDiv[x] != curPDiv) {
ret.push_back(pp);
curPDiv = minPDiv[x];
pp = curPDiv;
} else {
pp *= minPDiv[x];
}
x /= minPDiv[x];
}
ret.push_back(pp);
return ret;
}
// Hàm sắp xếp các phương trình theo cơ số của modulo,
// dùng trong sort
bool compCongruenceByModuloBase(
Congruence eq1, Congruence eq2) {
if (minPDiv[eq1.mod] == minPDiv[eq2.mod]) {
return eq1.mod < eq2.mod;
}
return minPDiv[eq1.mod] < minPDiv[eq2.mod];
}
// Tách các phương trình thành các phương trình nhỏ hơn với modulo
// nguyên tố cùng nhau
bool factor(vector<Congruence>& eqSet) {
// Phân tích các modulo ra thừa số nguyên tố
vector<Congruence> pEqSet;
for (auto eq : eqSet) {
vector<int> ppDivs = poweredPrimeDivisors(eq.mod);
for (auto d : ppDivs) {
pEqSet.push_back(Congruence(eq.rem, d));
}
}
// Sắp xếp các phương trình theo cơ số của modulo
sort(pEqSet.begin(), pEqSet.end(),
compCongruenceByModuloBase);
// Kiểm tra hệ phương trình vô nghiệm
for (int i = 1; i < (int)pEqSet.size(); i++) {
auto eq1 = pEqSet[i - 1];
auto eq2 = pEqSet[i];
if (minPDiv[eq1.mod] == minPDiv[eq2.mod]) {
if (eq2.rem % eq1.mod != eq1.rem) {
return 0;
}
}
}
// Loại bỏ các modulo nhỏ
eqSet.clear();
for (int i = pEqSet.size() - 1; i >= 0; i--) {
if (eqSet.empty() || minPDiv[pEqSet[i].mod]
!= minPDiv[eqSet.back().mod]) {
eqSet.push_back(pEqSet[i]);
}
}
return 1;
}
Congruence solveCongruenceEqSet(vector<Congruence>& eqSet) {
if (!factor(eqSet)) {
return Congruence(0, 0);
}
Congruence sol(0, 1);
for (int i = 0; i < (int)eqSet.size(); i++) {
sol.mod *= eqSet[i].mod;
}
vector <long long> p;
vector <long long> pInv;
for (int i = 0; i < (int)eqSet.size(); i++) {
p.push_back(sol.mod / eqSet[i].mod);
pInv.push_back(invMod(p[i], eqSet[i].mod));
}
for (int i = 0; i < (int)eqSet.size(); i++) {
long long prod = 1;
prod = prod * p[i] % sol.mod;
prod = prod * pInv[i] % sol.mod;
prod = prod * eqSet[i].rem % sol.mod;
sol.rem += prod;
}
sol.rem %= sol.mod;
return sol;
}
Nói về độ phức tạp, ta xét các công việc con:
std::sort
) là Nếu không nhớ công thức nghiệm của hệ phương trình, ta vẫn có thể giải được hệ phương trình đồng dư chỉ bằng các phương trình Diophantus. Đặc biệt, cách làm này không hề bị ảnh hưởng bởi việc các modulo có nguyên tố cùng nhau đôi một hay không.
Với
Theo định lý Thặng dư Trung Hoa và những gì ta đã làm ở trường hợp modulo không nguyên tố cùng nhau đôi một, (**) có thể không có nghiệm có một họ nghiệm duy nhất theo modulo
Dựa vào hai phương trình thành phần của (**), ta có:
(1) là một phương trình Diophantus tuyến tính với hai ẩn là
là một nghiệm của (**).
Với
Để ý dạng của (***) chẳng khác gì (**) cả. Ta có thể áp dụng mọi biện pháp giải (**) để giải hệ này. Ta cứ tiếp tục như vậy tới khi toàn bộ phương trình đã được giải.
const pair <long long, long long>
INVALID_ROOT = {LLONG_MAX, LLONG_MAX};
Congruence solveInduction(vector <Congruence>& eqSet, int solved) {
Congruence sol(0, 1), lastSol(0, 1);
if (solved == 1) {
lastSol = eqSet[0];
} else {
lastSol = solveInduction(eqSet, solved - 1);
}
if (!lastSol.mod) {
return Congruence(0, 0);
}
// Tìm cặp số x, y để ax + by = gcd(a, b)
pair<long long, long long> p = {0, 0};
long long d = extEuclid(lastSol.mod, eqSet[solved].mod,
p.first, p.second);
// Điều kiện vô nghiệm của phương trình Diophantus
if ((eqSet[solved].rem - lastSol.rem) % d) {
return Congruence(0, 0);
}
// Đưa ra nghiệm đúng của phương trình Diophantus và
// tìm nghiệm của hệ đến phương trình hiện tại
sol.rem = p.first * (eqSet[solved].rem - lastSol.rem) / d
% (eqSet[solved].mod / d) * lastSol.mod
+ lastSol.rem;
sol.mod = lastSol.mod * eqSet[solved].mod / d;
sol.rem %= sol.mod;
if (sol.rem < 0) {
sol.rem += sol.mod;
}
return sol;
}
Congruence solveCongruenceEqSet(vector <Congruence>& eqSet) {
if (eqSet.size() == 1) {
return eqSet.front();
}
return solveInduction(eqSet, eqSet.size() - 1);
}
Về độ phức tạp, ta đã thực hiện giải
Nhìn chung, cách làm này có vẻ tự nhiên hơn so với việc sử dụng công thức trực tiếp. Tuy nhiên, khi cài đặt, cần phải đặc biệt chú ý xử lý các phép toán lấy số dư, đặc biệt khi giải phương trình Diophantus để tìm nghiệm.