¶ Định lý Thặng dư Trung Hoa (CRT)
Tác giả:
- Nguyễn Đức Kiên - Đại học Công nghệ - ĐHQGHN
Reviewer:
- Lê Minh Hoàng - Đại học Khoa học Tự nhiên, ĐHQG-HCM
- Nguyễn Minh Nhật - Trường THPT chuyên Khoa học Tự nhiên, ĐHQGHN
- Đoàn Gia Huy - Trường THPT chuyên Khoa học Tự nhiên, ĐHQGHN
¶ Kiến thức cần biết
Để 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ề:
¶ Mở đầu
Ả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 , cộng thêm rồi kiểm tra các điều kiện còn lại. Trình Đại Vỹ thời Minh có nhắc đến một lời giải có thể tóm tắt lại bằng công thức với là một số nguyên phù hợp với thực tế.
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 toán Giải hệ phương trình đồng dư
Đề bài: LQDOJ - Hàn Tín điểm binh
Giải hệ gồm phương trình đồng dư ẩn như sau:
trong đó là các số nguyên bất kỳ sao cho với mọi .
Dễ thấy nếu là một nghiệm của (*) thì các giá trị cũng là nghiệm của (*). Vì vậy, thông thường các bài toán chỉ yêu cầu ta tìm nghiệm dương nhỏ nhất thoả mãn phương trình.
¶ Định lý thặng dư Trung Hoa
¶ Các biến phụ được sử dụng
Để 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:
- là tích của các với :
- là bội chung nhỏ nhất của các với : .
- với . Nói cách khác, là tích của các với
- là một giá trị nghịch đảo modulo của với :
¶ Ý tưởng
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:
- Mỗi số hạng chỉ thoả mãn một phương trình tương ứng với nó
- Mỗi số hạng đều chia hết cho tất cả các modulo của các phương trình khác. Như vậy nó sẽ không gây ảnh hưởng đến số dư của các phương trình còn lại
Nói cách khác, nghiệm của ta sẽ có dạng , trong đó mỗi đều thoả mãn:
¶ Xây dựng các số hạng thoả mãn điều kiện thứ nhất
Xét phương trình
Với ý tưởng xây dựng nghiệm như đã nói ở trên, cần phải thoả mãn (1). Dễ thấy thoả mãn (1), nên ta sẽ bắt đầu với .
¶ Xây dựng các số hạng thoả mãn điều kiện thứ hai
Để thoả mãn điều kiện thứ hai, cách đơn giản nhất là nhân hết các modulo khác vào nghiệm. Như vậy ta thu được
Sau khi nhân xong, sẽ không có gì xảy ra nếu . Nhưng nếu , số dư của khi chia cho cũng bị nhân lên một lượng tương ứng và không còn là nữa. Lúc này, ta cần tìm được một số nguyên sao cho:
Do , số như vậy chắc chắn tồn tại; số đó chính là . Việc nhân thêm vào không làm thay đổi tính chia hết của với các modulo khác . Tới bước này, ta đã có .
¶ Xây dựng nghiệm của hệ
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 còn lại. Cộng các số hạng đó lại, ta sẽ được một nghiệm thoả mãn phương trình là
Nhắc lại, nếu một số là nghiệm của (*) thì mọi số cũng sẽ là nghiệm của (*). Thêm vào đó, do các modulo nguyên tố cùng nhau đôi một nên . Như vậy, (*) sẽ có nghiệm là:
Đây chính là công thức nghiệm được đưa ra bởi định lý Thặng dư Trung Hoa.
¶ Nội dung đị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 (nhấn để hiện)
Chứng minh sự tồn tại: Do là các số đôi một nguyên tố cùng nhau nên dễ thấy với mọi thì . Do và nguyên tố cùng nhau nên tồn tại nghịch đảo modulo của p_i, chính là . Vì nên
Xét . Ta thấy với mọi thì do là tích của tất cả các với khác . Vì vậy, với mọi ta có
Vậy là một nghiệm của (*).
Chứng minh sự duy nhất: Giả sử và là hai số nguyên thoả mãn (*). Với mọi , ta đều có . Lấy bội chung nhỏ nhất của tất cả các đồng dư thức dạng trên ta được
Do các số nguyên tố cùng nhau đôi một nên . Như vậy
Vậy (*) có nghiệm duy nhất theo modulo .
¶ Cài đặt
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à xong.
Lưu ý rằng, không phải lúc nào cũng là một số nguyên tố, nên ta không thể áp dụng định lý Fermat nhỏ để tìm . Thay vào đó, ta phải giải phương trình Diophantus tuyến tính sau với hai ẩn bằng thuật toán Euclid mở rộng:
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;
}
¶ Độ phức tạp
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 (vì ) cho mỗi lần giải phương trình Diophantus và cho các việc còn lại, do vậy ta sẽ mất để giải hệ phương trình.
¶ Trường hợp modulo không nguyên tố cùng nhau
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ị không chắc chắn tồn tại nữa, và cũng không còn khẳng định được phương trình có nghiệm.
Ta đã biết đồng dư thức có hai tính chất sau với các số nguyên bất kỳ sao cho :
- thì , với là một ước bất kỳ của
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ử được phân tích thành với là các số nguyên tố. Khi đó, phương trình sẽ tương đương với hệ các phương trình sau:
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ó , trong đó là một số nguyên tố và , thì hệ gồm hai phương trình trên sẽ chỉ có nghiệm nếu
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 kết hợp lại sẽ tương đương với phương trình có số mũ lớn nhất (xem phần ví dụ để hiểu thêm). Ta chỉ cần giữ lại phương trình này. Sau khi thực hiện bước này, ta đã có một hệ gồm các phương trình có modulo nguyên tố cùng nhau. Việc cần làm lúc này là áp dụng định lý Thặng dư Trung Hoa để giải hệ.
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 , ta có . Như vậy không có mâu thuẫn nào xảy ra. Ta loại đi phương trình thứ nhất do đã có modulo . Áp dụng định lý Thặng dư Trung Hoa cho hệ này, ta có nghiệm của phương trình là .
Đ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:
- Để phân tích các modulo ra thừa số nguyên tố, sẽ có rất nhiều cách khác nhau. Nếu đặt , cách làm trên có độ phức tạp thời gian là . Còn nếu ta chỉ chuẩn bị các ước nguyên tố nhỏ nhất cho số nguyên đầu tiên, ta sẽ mất . Độ phức tạp không gian cũng vậy, hoàn toàn tuỳ thuộc vào cách chọn thuật toán phân tích. Về các thuật toán phân tích, có thể tham khảo tại đây.
- Sau khi phân tích, mỗi phương trình sẽ "nở" ra thành phương trình khác. Chúng ta có phương trình. Lúc này, sẽ xuất hiện các phương trình có modulo với cơ số giống nhau (cùng giá trị ) và có thể sẽ xuất hiện mâu thuẫn. Cách tốt nhất để kiểm tra và loại trừ các cặp phương trình mâu thuẫn là sắp xếp lại các phương trình có cùng cơ số theo độ lớn của số mũ và xét từng cặp kề nhau. Độ phức tạp của bước này cũng tuỳ vào cách bạn chọn thuật sắp xếp nhưng theo cách "phổ thông" nhất (
std::sort) là thời gian và không gian. - Sau khi sắp xếp và thu gọn hệ xong, trong trường hợp tệ nhất ta vẫn có phương trình (tuy nhiên thực tế số phương trình còn lại là nhỏ hơn rất nhiều). Như vậy, độ phức tạp của công việc này sẽ là thời gian và không gian.
¶ Cách tiếp cận khác
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 , ta có hệ phương trình:
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 . Như vậy, nghiệm ta cần tìm có dạng
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à và . Nếu nó có nghiệm,
là một nghiệm của (**).
Với lớn, ta giải từng cặp phương trình một bằng cách quy nạp. Chẳng hạn, giả sử ta đã giải xong phương trình đầu tiên và tìm được nghiệm . Khi đó ta sẽ giải hệ
Để ý 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 phương trình Diophantus. Làm như vậy sẽ tốn thời gian và bộ nhớ.
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.
¶ Chú ý thêm
- Định lý thặng dư Trung Hoa được ứng dụng trong các thuật toán sau:
- Tính tổ hợp moldulo với không nguyên tố bằng cách phân tích modulo ra thừa số nguyên tố, tính modulo cho các thừa số nguyên tố rồi kết hợp lại.
- Thuật toán Garner: cho phép biểu diễn và tính toán trên các số rất lớn
¶ Bài tập áp dụng
- Kattis - Chinese Remainder Theorem (non-relatively prime moduli): bài toán chỉ có hai phương trình với modulo không nguyên tố
- LightOJ - Monkey tradition
- Hackerank - Number of sequences
- Hackerank - nCr: áp dụng định lý Thặng dư Trung Hoa để tính tổ hợp modulo
- Codeforces - Billards
- UVa Online Judge - Code Feat
- POI 2016 - Arkanoid
- Codeforces - Remainders Game
- Google Code Jam - Golf Gopher
¶ Tài liệu tham khảo
- CP Algorithms - Chinese Remainder Theorem
- USACO Guide - Extended Euclidean Algorithm
- Codeforces - [Tutorial] Chinese Remainder Theorem (Tác giả: Valiors)
- Tạp chí Pi (cho phần lịch sử của thuật toán)