1. 項(xiàng)目概述從一道經(jīng)典算法題看高斯消元與模運(yùn)算的融合最近在翻看POJPKU Online Judge上的老題又遇到了那道經(jīng)典的2065SETI。這道題第一次做還是好多年前當(dāng)時(shí)對(duì)高斯消元解線性方程組的理解還停留在解實(shí)數(shù)域上的問題一看到要對(duì)系數(shù)取模就有點(diǎn)發(fā)怵?，F(xiàn)在回頭看它其實(shí)是一個(gè)絕佳的案例把數(shù)論里的模運(yùn)算和線性代數(shù)里的高斯消元法巧妙地結(jié)合在了一起。題目本身模擬了一個(gè)簡(jiǎn)化的SETI搜尋地外文明信號(hào)解碼場(chǎng)景給你一個(gè)基于質(zhì)數(shù)模數(shù)P的方程系統(tǒng)要求你解出每個(gè)變量可以理解為“字母”的值。這不僅僅是算法競(jìng)賽中的一道題其背后“在有限域上求解線性方程組”的思想在密碼學(xué)、編碼理論乃至通信系統(tǒng)的糾錯(cuò)解碼中都有實(shí)實(shí)在在的應(yīng)用。今天我就結(jié)合這道題把高斯消元解模方程模數(shù)為質(zhì)數(shù)的完整思路、實(shí)現(xiàn)細(xì)節(jié)、以及我踩過的那些坑系統(tǒng)地梳理一遍。無論你是正在備賽的選手還是對(duì)算法如何應(yīng)用于實(shí)際問題感興趣的開發(fā)者相信這篇都能給你帶來一些直接的參考。2. 問題背景與數(shù)學(xué)模型建立2.1 題意解析與問題轉(zhuǎn)化POJ 2065的題目描述大致是這樣的我們接收到一個(gè)長(zhǎng)度為n的字符串str它實(shí)際上是由一個(gè)函數(shù)f(k)生成。函數(shù)定義為對(duì)于第i個(gè)位置i從0開始有f(i) str[i]如果str[i]是*則對(duì)應(yīng)值為0否則為str[i] - a 1。同時(shí)存在一個(gè)關(guān)于變量a[0], a[1], ..., a[n-1]的線性方程組對(duì)于每一個(gè)i0 i n滿足sum_{j0}^{n-1} (a[j] * (i1)^j) ≡ f(i) (mod P)其中P是一個(gè)給定的質(zhì)數(shù)。我們需要求解的就是這個(gè)n元一次模線性方程組的解即每個(gè)a[j]在模P意義下的值。舉個(gè)例子如果字符串是abcP31那么n3。f(0)對(duì)應(yīng)字符a值為1f(1)對(duì)應(yīng)字符b值為2f(2)對(duì)應(yīng)字符c值為3。方程組就是a0 * 1^0 a1 * 1^1 a2 * 1^2 ≡ 1 (mod 31) a0 * 2^0 a1 * 2^1 a2 * 2^2 ≡ 2 (mod 31) a0 * 3^0 a1 * 3^1 a2 * 3^2 ≡ 3 (mod 31)化簡(jiǎn)后1*a0 1*a1 1*a2 ≡ 1 1*a0 2*a1 4*a2 ≡ 2 1*a0 3*a1 9*a2 ≡ 3這樣問題就清晰地轉(zhuǎn)化為了求解一個(gè)模PP為質(zhì)數(shù)意義下的n階線性方程組A * X ≡ B (mod P)。其中系數(shù)矩陣A的第i行第j列元素為(i1)^j mod P常數(shù)向量B的第i個(gè)元素為f(i)未知向量X就是我們要找的a[0..n-1]。2.2 為什么高斯消元法依然適用在實(shí)數(shù)域上我們使用高斯消元法解線性方程組核心是通過行初等變換交換兩行、某行乘以非零常數(shù)、一行加上另一行的倍數(shù)將增廣矩陣化為行階梯形或最簡(jiǎn)行階梯形然后回代求解。這里有一個(gè)關(guān)鍵這些行變換必須是“可逆”的且不能改變方程組的解集。在模PP為質(zhì)數(shù)的有限域上這些性質(zhì)依然成立但運(yùn)算規(guī)則變成了模運(yùn)算。這就要求我們?cè)诿恳徊讲僮髦卸急仨毧紤]模P下的逆元。交換兩行顯然可行不涉及乘除。某行乘以一個(gè)非零常數(shù)k要求k在模P下存在逆元k^{-1}這樣變換才是可逆的。因?yàn)镻是質(zhì)數(shù)所以任何1到P-1的整數(shù)在模P下都有逆元由費(fèi)馬小定理或擴(kuò)展歐幾里得算法求得。因此只要k不是0模P意義下我們就可以進(jìn)行這個(gè)操作。將第i行加上第j行的k倍這相當(dāng)于對(duì)第i行做了一個(gè)線性組合。只要這個(gè)操作本身是定義良好的即k在模P下有意義它就是可逆的其逆操作是減去第j行的k倍因此也不會(huì)改變解集。所以只要模數(shù)P是質(zhì)數(shù)我們就能保證在消元過程中用來做主元的非零系數(shù)稱為主元都存在模逆元從而可以像實(shí)數(shù)域上一樣進(jìn)行歸一化將主元變?yōu)?和消去其他行對(duì)應(yīng)列的元素。這就是解模線性方程組的高斯消元法常被稱為高斯-約當(dāng)消元法能夠成立的理論基礎(chǔ)。如果P不是質(zhì)數(shù)那么不是所有非零數(shù)都有逆元消元過程可能會(huì)失敗或需要更復(fù)雜的處理比如模合數(shù)下的解方程通常分解質(zhì)因數(shù)后用中國(guó)剩余定理組合本題保證了P是質(zhì)數(shù)大大簡(jiǎn)化了問題。3. 算法核心模P意義下的高斯-約當(dāng)消元法詳解3.1 算法流程與實(shí)數(shù)域消元的異同整體流程和實(shí)數(shù)域上的高斯-約當(dāng)消元法幾乎一致目標(biāo)都是將增廣矩陣化為行最簡(jiǎn)形式每一行只有一個(gè)主元1且該1所在的列其他元素全為0。不同點(diǎn)全部體現(xiàn)在具體的運(yùn)算上所有加減乘除都必須模P。假設(shè)我們有n個(gè)方程n個(gè)未知數(shù)增廣矩陣aug大小為n x (n1)其中前n列是系數(shù)矩陣A最后一列是常數(shù)向量B。算法步驟如下初始化設(shè)當(dāng)前列col 0當(dāng)前行row 0。尋找主元遍歷第col列從第row行開始找到一個(gè)aug[r][col] % P ! 0的行r。如果找不到說明這一列是自由變量在模P意義下全為0理論上應(yīng)該處理自由變量但根據(jù)POJ 2065的題目描述該方程組保證有唯一解所以我們可以認(rèn)為總能找到主元。找到后交換第r行和第row行。主元?dú)w一化令pivot aug[row][col]。因?yàn)閜ivot非零且P是質(zhì)數(shù)所以pivot在模P下存在逆元inv_pivot。我們將第row行的每一個(gè)元素aug[row][j]都乘以inv_pivot然后對(duì)P取模。這樣aug[row][col]就變成了1。注意這里必須先求逆元再乘而不是直接除以pivot。因?yàn)槟＿\(yùn)算下沒有直接的除法除法需要用乘以逆元來實(shí)現(xiàn)。inv_pivot可以通過擴(kuò)展歐幾里得算法或費(fèi)馬小定理pow(pivot, P-2, P)快速計(jì)算。消去其他行對(duì)于所有非當(dāng)前行ii從0到n-1且i ! row計(jì)算倍數(shù)factor aug[i][col]因?yàn)榇藭r(shí)aug[row][col]已經(jīng)是1了。然后將第i行的每一個(gè)元素aug[i][j]減去factor * aug[row][j]并對(duì)結(jié)果取模P。這一步會(huì)使得第col列上除了第row行是1其他行都變?yōu)?。注意這里的減法和乘法都要進(jìn)行模P運(yùn)算防止中間結(jié)果溢出尤其是在用C/C等語言實(shí)現(xiàn)時(shí)。通常我們會(huì)在每一步運(yùn)算后立即取模。移動(dòng)指針完成上述操作后row,col準(zhǔn)備處理下一列。循環(huán)重復(fù)步驟2-5直到row或col達(dá)到n。提取解此時(shí)增廣矩陣已經(jīng)被化為行最簡(jiǎn)形式。方程組的解X[j]就直接存儲(chǔ)在aug[j][n]即變換后的常數(shù)項(xiàng)列中。因?yàn)榈趈行的主元在第j列值為1所以a[j] ≡ aug[j][n] (mod P)。3.2 關(guān)鍵工具模逆元的計(jì)算這是整個(gè)算法實(shí)現(xiàn)中的基石。有兩種常見方法擴(kuò)展歐幾里得算法 (Extended Euclidean Algorithm)求解a * x P * y 1的整數(shù)解x這個(gè)x模P就是a的逆元。這是通用方法即使P不是質(zhì)數(shù)只要gcd(a, P)1就能用。// 返回 a 在模 mod 下的逆元假設(shè) gcd(a, mod) 1 long long inv(long long a, long long mod) { long long x, y; exgcd(a, mod, x, y); // 解 ax mod*y 1 return (x % mod mod) % mod; // 調(diào)整到 0~mod-1 范圍內(nèi) }費(fèi)馬小定理 (Fermat‘s Little Theorem)當(dāng)P為質(zhì)數(shù)且a不是P的倍數(shù)時(shí)有a^(P-1) ≡ 1 (mod P)。因此a的逆元就是a^(P-2) mod P?？梢杂每焖賰缬?jì)算。// 快速冪求模逆元僅當(dāng) mod 為質(zhì)數(shù)時(shí)可用 long long inv(long long a, long long mod) { return pow_mod(a, mod - 2, mod); // pow_mod 是快速冪取模函數(shù) }在POJ 2065的場(chǎng)景下P是質(zhì)數(shù)兩種方法都可以。由于P不大題目說小于30000且我們需要頻繁求逆元每選一個(gè)主元就要算一次使用擴(kuò)展歐幾里得通常更穩(wěn)定且常數(shù)可能更小?？焖賰缧枰鰞邕\(yùn)算如果P很大比如1e97用快速冪配合快速冪函數(shù)也很高效。在實(shí)際編碼中我更喜歡用擴(kuò)展歐幾里得感覺更“基礎(chǔ)”一些。3.3 代碼實(shí)現(xiàn)框架與注釋這里給出一個(gè)用C風(fēng)格描述的核心消元函數(shù)框架包含了上述所有要點(diǎn)#include iostream #include vector #include cmath using namespace std; typedef long long ll; // 擴(kuò)展歐幾里得算法 ll exgcd(ll a, ll b, ll x, ll y) { if (b 0) { x 1; y 0; return a; } ll d exgcd(b, a % b, y, x); y - (a / b) * x; return d; } // 求a在模mod下的逆元假設(shè)gcd(a, mod)1 ll inv(ll a, ll mod) { ll x, y; exgcd(a, mod, x, y); return (x % mod mod) % mod; } // 高斯消元解模線性方程組 (模數(shù)P為質(zhì)數(shù)) // aug為增廣矩陣n為方程數(shù)/未知數(shù)個(gè)數(shù)P為模數(shù) // 返回解向量X如果無解或多解本題保證唯一解可調(diào)整返回值類型 vectorll gauss_mod(vectorvectorll aug, ll P) { int n aug.size(); // 行數(shù) vectorll X(n, 0); for (int col 0, row 0; col n row n; col) { // 1. 選主元找第col列中第row行以下第一個(gè)非零元 int pivot row; while (pivot n aug[pivot][col] % P 0) { pivot; } if (pivot n) { continue; // 這一列全為0理論上可能是自由變量本題忽略 } // 2. 交換行 swap(aug[row], aug[pivot]); // 3. 主元?dú)w一化 ll val aug[row][col] % P; ll inv_val inv(val, P); // 核心計(jì)算模逆元 for (int j col; j n; j) { // 注意要處理到常數(shù)項(xiàng)列 aug[row][j] (aug[row][j] % P) * inv_val % P; } // 4. 消去其他行 for (int i 0; i n; i) { if (i ! row aug[i][col] ! 0) { ll factor aug[i][col] % P; for (int j col; j n; j) { // 模運(yùn)算下的減法消元: aug[i][j] - factor * aug[row][j] aug[i][j] (aug[i][j] - factor * aug[row][j]) % P; // 保證結(jié)果非負(fù) if (aug[i][j] 0) aug[i][j] P; } } } row; } // 5. 提取解 for (int i 0; i n; i) { // 行最簡(jiǎn)形式下第i行的主元在第i列解就在常數(shù)項(xiàng)列 X[i] aug[i][n] % P; if (X[i] 0) X[i] P; // 調(diào)整為非負(fù) } return X; }4. 針對(duì)POJ 2065的完整解題實(shí)現(xiàn)與優(yōu)化4.1 輸入處理與系數(shù)矩陣構(gòu)建有了上面的通用消元函數(shù)解決POJ 2065就剩下構(gòu)建增廣矩陣了。根據(jù)題目描述讀入測(cè)試用例數(shù)T。對(duì)每個(gè)用例讀入模數(shù)P和字符串str。n str.length()。構(gòu)建n x (n1)的增廣矩陣aug。aug[i][j] (i1)^j mod P。這里(i1)^j可能很大需要用快速冪取模計(jì)算。aug[i][n] f(i) mod P。f(i)根據(jù)str[i]計(jì)算*為0否則為str[i] - a 1。這里有一個(gè)性能優(yōu)化點(diǎn)計(jì)算(i1)^j mod P。如果對(duì)每個(gè)i, j都單獨(dú)用快速冪計(jì)算復(fù)雜度是O(n^3 log P)對(duì)于n最大為70的數(shù)據(jù)范圍POJ典型范圍雖然可以接受但不夠優(yōu)雅。我們可以利用遞推對(duì)于固定的i(i1)^0 1。(i1)^j (i1)^(j-1) * (i1) mod P。 這樣對(duì)于每一行i我們可以在O(n)時(shí)間內(nèi)計(jì)算出所有j對(duì)應(yīng)的系數(shù)整體構(gòu)建矩陣的復(fù)雜度就降到了O(n^2)。4.2 完整AC代碼參考與逐行分析下面是我在POJ上通過的代碼結(jié)合了上述所有討論#include iostream #include string #include vector #include cmath #include algorithm using namespace std; typedef long long ll; ll P; // 模數(shù)全局變量方便傳遞 // 快速冪取模用于計(jì)算系數(shù) (i1)^j % P ll pow_mod(ll a, ll b, ll mod) { ll res 1; a % mod; while (b) { if (b 1) res (res * a) % mod; a (a * a) % mod; b 1; } return res; } // 擴(kuò)展歐幾里得求逆元 ll exgcd(ll a, ll b, ll x, ll y) { if (!b) { x 1, y 0; return a; } ll d exgcd(b, a % b, y, x); y - a / b * x; return d; } ll inv(ll a, ll mod) { ll x, y; exgcd(a, mod, x, y); return (x % mod mod) % mod; } // 模意義下的高斯消元 vectorll gauss(vectorvectorll aug) { int n aug.size(); for (int col 0, row 0; col n row n; col) { // 選主元 int pivot row; while (pivot n aug[pivot][col] % P 0) pivot; if (pivot n) continue; swap(aug[row], aug[pivot]); // 歸一化 ll div inv(aug[row][col], P); for (int j col; j n; j) { aug[row][j] aug[row][j] * div % P; } // 消元 for (int i 0; i n; i) { if (i ! row aug[i][col]) { ll factor aug[i][col]; for (int j col; j n; j) { aug[i][j] (aug[i][j] - factor * aug[row][j]) % P; if (aug[i][j] 0) aug[i][j] P; // 立即調(diào)整非負(fù) } } } row; } // 回代實(shí)際上消元完成后解已在常數(shù)項(xiàng)列 vectorll res(n); for (int i 0; i n; i) { res[i] aug[i][n] % P; } return res; } int main() { int T; cin T; while (T--) { cin P; string s; cin s; int n s.size(); vectorvectorll aug(n, vectorll(n 1, 0)); // 構(gòu)建增廣矩陣 for (int i 0; i n; i) { // 計(jì)算常數(shù)項(xiàng) f(i) ll b (s[i] *) ? 0 : (s[i] - a 1); aug[i][n] b % P; // 計(jì)算系數(shù) (i1)^j % P使用遞推優(yōu)化 ll base (i 1) % P; ll cur 1; // (i1)^0 for (int j 0; j n; j) { aug[i][j] cur; cur (cur * base) % P; // 遞推計(jì)算下一個(gè)冪 } } vectorll ans gauss(aug); for (int i 0; i n; i) { cout ans[i] (i n - 1 ? \n : ); } } return 0; }逐行分析關(guān)鍵點(diǎn)pow_mod函數(shù)雖然我們?cè)跇?gòu)建矩陣時(shí)用了更優(yōu)的遞推但這個(gè)快速冪函數(shù)是基礎(chǔ)工具保留以備不時(shí)之需。inv函數(shù)使用擴(kuò)展歐幾里得求逆元這是模運(yùn)算消元的核心。gauss函數(shù)嚴(yán)格遵循了之前描述的流程。注意if (pivot n) continue;這一句在本題有唯一解的保證下遇到全零列直接跳過處理下一列是安全的。如果是一般情況這里需要更復(fù)雜的處理來判斷無解或多解。立即取模與調(diào)整非負(fù)在消元循環(huán)aug[i][j] (aug[i][j] - factor * aug[row][j]) % P;之后立刻跟一個(gè)if (aug[i][j] 0) aug[i][j] P;。這是防止負(fù)數(shù)模運(yùn)算產(chǎn)生意外結(jié)果的關(guān)鍵習(xí)慣。在C/C中-1 % 5的結(jié)果是-1而不是4我們必須手動(dòng)調(diào)整到[0, P-1]的范圍。主循環(huán)中的遞推在構(gòu)建矩陣的雙重循環(huán)里內(nèi)層循環(huán)用cur變量遞推計(jì)算(i1)^j避免了重復(fù)計(jì)算快速冪是一個(gè)有效的常數(shù)優(yōu)化。輸出格式注意行末空格的處理POJ比較嚴(yán)格。5. 常見陷阱、調(diào)試技巧與擴(kuò)展思考5.1 實(shí)戰(zhàn)中踩過的坑負(fù)數(shù)取模問題這是最大的坑沒有之一。C/C中%運(yùn)算符對(duì)負(fù)數(shù)的處理是“商向零取整”導(dǎo)致-1 % 5 -1。而在數(shù)論和本題中我們需要的是“最小非負(fù)剩余”即-1 ≡ 4 (mod 5)。所以在任何可能產(chǎn)生負(fù)數(shù)的運(yùn)算特別是減法后必須立即判斷并加模數(shù)調(diào)整到非負(fù)范圍。我的代碼中在消元后和提取解后都做了這個(gè)操作。中間結(jié)果溢出即使模數(shù)P只有30000但(i1)^j在取模前可能非常大比如70^69直接計(jì)算會(huì)溢出long long。因此必須在運(yùn)算過程中步步取模。我們遞推計(jì)算系數(shù)時(shí)cur (cur * base) % P;以及消元時(shí)的乘法和減法都立即跟了% P。逆元不存在雖然題目保證P是質(zhì)數(shù)但你的代碼是否處理了aug[row][col]為0的情況在選主元階段如果找到的pivot行該列元素為0歸一化時(shí)計(jì)算inv(0, P)就會(huì)出錯(cuò)0沒有逆元。所以選主元的循環(huán)判斷條件必須是aug[pivot][col] % P ! 0而不是簡(jiǎn)單的aug[pivot][col] ! 0因?yàn)橐粋€(gè)數(shù)可能是P的倍數(shù)模P后為0但原值不為0。浮點(diǎn)數(shù)誤區(qū)絕對(duì)不要試圖用double類型來解模方程有些初學(xué)者想先按實(shí)數(shù)解再取模。這是完全錯(cuò)誤的因?yàn)楦↑c(diǎn)數(shù)的精度誤差和取模運(yùn)算的離散性會(huì)導(dǎo)致結(jié)果毫無意義。必須全程使用整數(shù)運(yùn)算。時(shí)間復(fù)雜度估算高斯消元是O(n^3)的n70時(shí)計(jì)算量大約在70^3343,000次運(yùn)算級(jí)別加上模運(yùn)算完全在合理范圍內(nèi)。但如果n達(dá)到幾百就需要考慮優(yōu)化如使用bitset優(yōu)化異或方程組或迭代法不過這是后話了。5.2 調(diào)試技巧與測(cè)試數(shù)據(jù)當(dāng)你覺得代碼邏輯沒錯(cuò)但WAWrong Answer時(shí)可以按以下步驟排查小數(shù)據(jù)測(cè)試自己構(gòu)造n1,2,3的案例。比如P31, sa那么方程就是a0 ≡ 1 (mod 31)解應(yīng)為1。再比如P3, sab方程組是1*a0 1*a1 ≡ 1 (mod 3) 1*a0 2*a1 ≡ 2 (mod 3)手算一下解應(yīng)該是a00, a11。用你的程序跑一下看對(duì)不對(duì)。打印中間矩陣在消元的關(guān)鍵步驟如選主元后、歸一化后、消元后打印出增廣矩陣。對(duì)比手算過程看數(shù)值是否正確特別是模P后的值是否在[0, P-1]范圍內(nèi)。驗(yàn)證解求出解向量X后再代回原方程組sum(a[j]*(i1)^j) mod P看是否等于f(i)。這是一個(gè)非常有效的驗(yàn)證手段可以寫一個(gè)簡(jiǎn)單的驗(yàn)證函數(shù)。邊界測(cè)試測(cè)試P2的情況最小的質(zhì)數(shù)以及字符串全為*常數(shù)項(xiàng)全為0的情況看程序是否能正確處理。這里提供一個(gè)簡(jiǎn)單的測(cè)試用例輸入 2 31 abc 3 ab 輸出 1 2 3 0 1第一個(gè)用例就是我們最開始舉的例子。第二個(gè)用例上面分析過。5.3 從POJ 2065延伸出去的思考解出這道題不僅僅是AC了一道OJ題。它給你裝備了一套處理“模質(zhì)數(shù)線性方程組”的完整工具箱。這套工具能用在很多地方密碼學(xué)一些公鑰密碼算法如RSA的某些變種或線性同余生成器的分析中可能會(huì)涉及模方程組。編碼理論在糾錯(cuò)碼如Reed-Solomon碼的解碼過程中需要求解有限域上的線性方程組來恢復(fù)原始信息。算法競(jìng)賽這是處理一類“模運(yùn)算線性方程組”問題的標(biāo)準(zhǔn)板子。比如一些涉及組合數(shù)取模、多項(xiàng)式插值取模的問題最終可能歸結(jié)為此類方程。理解有限域通過親手實(shí)現(xiàn)你能更深刻地理解“在有限域上加減乘除除即乘逆元構(gòu)成一個(gè)封閉的代數(shù)系統(tǒng)”這一概念這與實(shí)數(shù)域或復(fù)數(shù)域有很大不同。更進(jìn)一步你可以思考如果模數(shù)P不是質(zhì)數(shù)而是合數(shù)該怎么辦通常的思路是將P分解質(zhì)因數(shù)對(duì)每個(gè)質(zhì)因子冪p^k求解方程組然后用中國(guó)剩余定理CRT將解組合起來。這要復(fù)雜得多因?yàn)槟^k下不是所有非零元都有逆元消元時(shí)需要更小心地處理。如何判斷模方程組無解或有多個(gè)解這需要在高斯消元完成后檢查系數(shù)矩陣的秩和增廣矩陣的秩。具體來說如果消元后出現(xiàn)0 non-zero (mod P)的行則無解如果系數(shù)矩陣的秩小于未知數(shù)個(gè)數(shù)則有無窮多解自由變量。如何優(yōu)化大規(guī)模模方程組的求解O(n^3)的復(fù)雜度對(duì)于n500可能就力不從心了。在實(shí)際工程中可能會(huì)利用系數(shù)矩陣的特殊性如稀疏性、對(duì)稱性、正定性等使用迭代法如共軛梯度法在有限域上的變種或?qū)Ｓ盟惴?。把這些都搞明白你對(duì)線性代數(shù)和數(shù)論在計(jì)算機(jī)中的應(yīng)用就算真正入門了。最后代碼實(shí)現(xiàn)時(shí)養(yǎng)成好習(xí)慣模運(yùn)算后立即調(diào)整非負(fù)、小心處理逆元、用遞推避免重復(fù)計(jì)算。這些細(xì)節(jié)決定了你的程序是優(yōu)雅地AC還是在WA和TLE中掙扎。