矩陣運算
高斯消去法
題目
假設有一方程式
\[
\begin{cases}
2x + 3y + 5z \equiv 69 \pmod{101} \\
0x + 6y + 4z \equiv 59 \pmod{101} \\
1x + 7y + 2z \equiv 89 \pmod{101} \\
\end{cases}
\]
我們要如何解呢?
如果使用高斯消去法我們可以得到
\[
\begin{cases}
1x + 0y + 0z \equiv x \pmod{101} \\
0x + 1y + 0z \equiv y \pmod{101} \\
0x + 0y + 1z \equiv z \pmod{101} \\
\end{cases}
\]
此時常數項就是答案
接下來就一步驟一步驟做高斯消去法
計算
在計算之前要注意模運算結果必為正數，參見運算的符號
先將原本的轉換為矩陣：
\[
\left[\begin{array}{ccc|c}
2 & 3 & 5 & 69 \\
0 & 6 & 4 & 59 \\
1 & 7 & 2 & 89 \\
\end{array}\right]
\]
我們先將第一行第一個換成1，也就是第一個的全部都要除以2，但是要進行模運算，所以是乘上模逆元就達到除法的效果
2模101的逆元是51，所以全部乘上51再模101。
變成：
\[
\left[\begin{array}{ccc|c}
1 & 52 & 53 & 85 \\
0 & 6 & 4 & 59 \\
1 & 7 & 2 & 89 \\
\end{array}\right]
\]
接著要把其他的x都變成0，第二列已經為0，第三列則是加上-1倍的第一列。
第三列：
((7-52)%101+101)%101=56
((2-53)%101+101)%101=56
\[
\left[\begin{array}{ccc|c}
1 & 52 & 53 & 85 \\
0 & 6 & 4 & 59 \\
0 & 56 & 50 & 4 \\
\end{array}\right]
\]
接著要把第二列第二個改為1，全部乘上17(6模101的逆元為17)再模101，
\[
\left[\begin{array}{ccc|c}
1 & 52 & 53 & 85 \\
0 & 1 & 68 & 94 \\
0 & 56 & 50 & 4 \\
\end{array}\right]
\]
接著要把其他的x都變成0，第一列要加上-52倍的第二列，第三列要加上45倍的第二列。
第一列：
((53 - 68 * 52)%101+101)%101=52
((85 - 94 * 52)%101+101)%101=45
第三列：
((50 + 68 * 45)%101+101)%101=80
((4 + 94 * 45)%101+101)%101=93
\[
\left[\begin{array}{ccc|c}
1 & 0 & 52 & 45 \\
0 & 1 & 68 & 94 \\
0 & 0 & 80 & 93 \\
\end{array}\right]
\]
接著要把第三列第三個改為1，全部乘上24(80模101的逆元為24)再模101，
\[
\left[\begin{array}{ccc|c}
1 & 0 & 52 & 45 \\
0 & 1 & 68 & 94 \\
0 & 0 & 1 & 10 \\
\end{array}\right]
\]
接著要把其他的x都變成0，第一列要加上49倍的第三列，第二列要加上-68倍的第三列。
第一列：
((45 - 10 * 52)%101+101)%101=30
第二列：
((94 - 10 * 68)%101+101)%101=20
\[
\left[\begin{array}{ccc|c}
1 & 0 & 0 & 30 \\
0 & 1 & 0 & 20 \\
0 & 0 & 1 & 10 \\
\end{array}\right]
\]
例外狀況-目標為0
假設今天遇到的方程式如下
\[
\begin{cases}
0x + 6y + 4z \equiv 59 \pmod{101} \\
2x + 3y + 5z \equiv 69 \pmod{101} \\
1x + 7y + 2z \equiv 89 \pmod{101} \\
\end{cases}
\]
可以看到第一列的第一個為0，我們沒辦法變成1，所以要和其他行交換，直接往下找，找到不是0的交換即可。
在這邊我們將第一列與第二列交換，如下。
\[
\begin{cases}
2x + 3y + 5z \equiv 69 \pmod{101} \\
0x + 6y + 4z \equiv 59 \pmod{101} \\
1x + 7y + 2z \equiv 89 \pmod{101} \\
\end{cases}
\]
程式碼實作
我們要先判斷是否有例外狀況，再依照計算步驟計算。
程式碼
#include <bits/stdc++.h>
using namespace std;
#define int long long
#define nn "\n"
#define nnn " "
#define nnnn "\n----------------------\n"
#define N 105
int n,mod;
int v[N][N];
int fpow(int x,int y){
if(!y)return 1;
if(y&1)return fpow(x,y-1)*x%mod;
int z=fpow(x,y/2)%mod;
return z*z%mod;
}
void out(){
cout<<nnnn;
for(int i=1;i<=n;i++){
for(int j=1;j<=n+1;j++){
cout<<v[i][j]<<nnn;
}
cout<<nn;
}
}
signed main(){
ios::sync_with_stdio(0);
cin.tie(0);
istringstream cin("3 101 \
69 59 89 \
5 4 2 \
3 6 7 \
2 0 1");
cin>>n>>mod;
for(int i=1;i<=n+1;i++){//存入陣列
for(int j=1;j<=n;j++){
cin>>v[j][n-i+2];
}
}
out();
for(int i=1;i<=n;i++){
//例外狀況交換(目標為0)
if(v[i][i]==0){
for(int j=i;j<=n;j++){//往後找不是0的
if(v[j][i]!=0){//找到不是0的
for(int k=1;k<=n+1;k++){//交換
swap(v[i][k],v[j][k]);
}
break;
}
}
}
//第i個變成1
int e=fpow(v[i][i],mod-2);
for(int j=i;j<=n+1;j++){//前面都是0了，往後乘模逆元就好
v[i][j]=v[i][j]*e%mod;
}
//其他變成0
for(int j=1;j<=n;j++){
if(j==i)continue;//跳過自己
int tmp=v[j][i];//要先存起來不然會被改掉
for(int k=i;k<=n+1;k++){
v[j][k]= ((v[j][k] - (v[i][k]*tmp))%mod+mod)%mod;
}
}
out();
}
}
輸出
------------------------------
2 3 5 69
0 6 4 59
1 7 2 89
------------------------------
1 52 53 85
0 6 4 59
0 56 50 4
------------------------------
1 0 52 45
0 1 68 94
0 0 80 93
------------------------------
1 0 0 30
0 1 0 20
0 0 1 10
例題：2170 . 地圖編修 (Map)
題解
#include <bits/stdc++.h>
using namespace std;
#define int long long
#define nn "\n"
#define nnn " "
#define nnnn "\n----------------------\n"
#define N 105
int n,mod;
int v[N][N];
int fpow(int x,int y){
if(!y)return 1;
if(y&1)return fpow(x,y-1)*x%mod;
int z=fpow(x,y/2)%mod;
return z*z%mod;
}
void out(){
cout<<nnnn;
for(int i=1;i<=n;i++){
for(int j=1;j<=n+1;j++){
cout<<v[i][j]<<nnn;
}
cout<<nn;
}
cout<<nnnn;
}
signed main(){
ios::sync_with_stdio(0);
cin.tie(0);
//istringstream cin("3 101 \
69 59 89 \
5 4 2 \
3 6 7 \
2 0 1");
cin>>n>>mod;
for(int i=1;i<=n+1;i++){//存入陣列
for(int j=1;j<=n;j++){
cin>>v[j][n-i+2];
}
}
for(int i=1;i<=n;i++){
//例外狀況交換(目標為0)
if(v[i][i]==0){
for(int j=i;j<=n;j++){//往後找不是0的
if(v[j][i]!=0){//找到不是0的
for(int k=1;k<=n+1;k++){//交換
swap(v[i][k],v[j][k]);
}
break;
}
}
}
//第i個變成1
int e=fpow(v[i][i],mod-2);
for(int j=i;j<=n+1;j++){//前面都是0了，往後乘模逆元就好
v[i][j]=v[i][j]*e%mod;
}
//其他變成0
for(int j=1;j<=n;j++){
if(j==i)continue;//跳過自己
int tmp=v[j][i];//要先存起來不然會被改掉
for(int k=i;k<=n+1;k++){
v[j][k]= ((v[j][k] - (v[i][k]*tmp))%mod+mod)%mod;
}
}
}
for(int i=n;i>=1;i--)cout<<v[i][n+1]<<" ";
}