高斯消元
高斯消元是一种能在 \(O(N^3)\) 的时间内求解 \(N\) 元一次方程组的算法。
其思路大致如下：
- 使第一个未知数只有第一个式子中系数非 \(0\)。
- 使第二个未知数只有第二个式子中系数非 \(0\)。
- \(\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \vdots\)
- 使第 \(N\) 个未知数只有第 \(N\) 个式子中系数非 \(0\)。
在这时第 \(i\) 个未知数的解就为等号右边的数除以自己的系数。如果有式子无法满足，则无解或有无穷多解。
但如果是这样呢？
\[\begin{cases} x_2=0\\ x_1=0 \end{cases} \]这样这种算法也会认为它不合法。所以当我们在让第 \(i\) 个未知数满足条件之前要先从后面的式子中找到一个满足条件的并交换即可。
代码
#include<bits/stdc++.h>
using namespace std;
const int MAXN = 105;
const long double EPS = 1e-6;
int n;
vector<long double> a[MAXN];
int main() {
ios::sync_with_stdio(false), cin.tie(0), cout.tie(0);
cin >> n;
for(int i = 1; i <= n; ++i) {
a[i].resize(n + 2);
for(int j = 1; j <= n + 1; ++j) {
cin >> a[i][j];
}
}
for(int i = 1; i <= n; ++i) {
for(int j = i; j <= n; ++j) {
if(fabs(a[j][i]) > EPS) {
swap(a[i], a[j]);
break;
}
}
if(fabs(a[i][i]) <= EPS) {
cout << "No Solution";
return 0;
}
for(int j = 1; j <= n; ++j) {
if(i != j) {
long double res = 1.0l * a[j][i] / a[i][i];
for(int k = 1; k <= n + 1; ++k) {
a[j][k] -= 1.0l * a[i][k] * res;
}
}
}
}
for(int i = 1; i <= n; ++i) {
cout << 'x' << i << '=' << fixed << setprecision(2) << 1.0l * a[i][n + 1] / a[i][i] << "\n";
}
return 0;
}
拓展
当我们需要区分无解和无穷解时该怎么办呢？
我们就需要在交换式子时枚举所有的式子，如果这个式子中不应为 \(0\) 的未知数却为 \(0\) 或该式子在当前式子之后，那么就交换，并且当当前式子对应未知数为 \(0\) 时就不管他，直接 continue;
等消元完成后，如果发现式子左边为 \(0\) 而右边不为 \(0\)，那么该方程组无解。如果发现式子左边为 \(0\) 而右边也为 \(0\)，那么该方程组有无限组解。
注意：只要有一个式子发现无解就无解，即使既有无限解的式子又有无解的式子也是这样。
代码
#include<bits/stdc++.h>
using namespace std;
const int MAXN = 105;
const long double EPS = 1e-6;
int n;
vector<long double> a[MAXN];
int main() {
ios::sync_with_stdio(false), cin.tie(0), cout.tie(0);
cin >> n;
for(int i = 1; i <= n; ++i) {
a[i].resize(n + 2);
for(int j = 1; j <= n + 1; ++j) {
cin >> a[i][j];
}
}
for(int i = 1; i <= n; ++i) {
for(int j = 1; j <= n; ++j) {
if((j >= i || fabs(a[j][j]) <= EPS) && fabs(a[j][i]) > EPS) {
swap(a[i], a[j]);
break;
}
}
if(fabs(a[i][i]) <= EPS) {
continue;
}
for(int j = 1; j <= n; ++j) {
if(i != j) {
long double res = 1.0l * a[j][i] / a[i][i];
for(int k = 1; k <= n + 1; ++k) {
a[j][k] -= 1.0l * a[i][k] * res;
}
}
}
}
for(int i = 1; i <= n; ++i) {
bool op = 1;
for(int j = 1; j <= n; ++j) {
op &= (fabs(a[i][j]) <= EPS);
}
if(op && fabs(a[i][n + 1]) > EPS) {
cout << "No Solution";
return 0;
}
}
for(int i = 1; i <= n; ++i) {
bool op = 1;
for(int j = 1; j <= n; ++j) {
op &= (fabs(a[i][j]) <= EPS);
}
if(op && fabs(a[i][n + 1]) <= EPS) {
cout << "Infinite Solution";
return 0;
}
}
for(int i = 1; i <= n; ++i) {
cout << 'x' << i << '=' << fixed << setprecision(2) << 1.0l * a[i][n + 1] / a[i][i] << "\n";
}
return 0;
}
矩阵快速幂
首先介绍矩阵是什么以及四则运算。
矩阵就是一个 \(N\times M\) 的二维数组。
加减法：两个矩阵要做加减法必须满足两个矩阵的行和列相等，只需对对应的元素进行加减即可。如 \(\begin{bmatrix}1\ \ 2\\3\ \ 4\end{bmatrix}+\begin{bmatrix}5\ \ 6\\7\ \ 8\end{bmatrix}=\begin{bmatrix}6\ \ 8\\10\ 12\end{bmatrix}\)
乘法：设矩阵 \(A\) 大小为 \(N\times M\)，矩阵 \(B\) 大小为 \(K\times N\)，\(C=AB\)，那么 \(C\) 的大小为 \(K\times M\) 且 \(C_{i,j}=\sum \limits_{k=1}^N A_{k,j}\cdot B_{i,k}\)。如 \(\begin{bmatrix}1\ \ 2\\3\ \ 4\end{bmatrix} \begin{bmatrix}5\ \ 6\\7\ \ 8\end{bmatrix}=\begin{bmatrix}23\ \ 34\\31\ \ 46\end{bmatrix}\)。
幂：设矩阵 \(A\) 大小为 \(N\times N\)，那么 \(A^x=\underbrace{AAA\dots A}_ x\)。如 \(\begin{bmatrix}1\ \ 2\\3\ \ 4\end{bmatrix}^3=\begin{bmatrix}37\ \ 54\\81\ \ 118\end{bmatrix}\)。这里我们定义 \(A^0=\begin{bmatrix}1\ \ 0\ \ 0\ \ \cdots\ \ 0\\0\ \ 1\ \ 0\ \ \cdots\ \ 0\\0\ \ 0\ \ 1\ \ \cdots\ \ 0\\\vdots\ \ \vdots\ \ \vdots\ \ \ddots \ \ \vdots\\0\ \ 0\ \ 0\ \ \cdots\ \ 1\end{bmatrix}\)，这个矩阵也是 \(N\times N\) 的。
代码
const int MAXN = 101;
struct Matrix {
int n, m, mat[MAXN][MAXN];
void clear(int a, int b) {
n = a, m = b;
for(int i = 1; i <= n; ++i) {
for(int j = 1; j <= m; ++j) {
mat[i][j] = 0;
}
}
}
void Set(int a) {
n = m = a;
for(int i = 1; i <= n; ++i) {
for(int j = 1; j <= m; ++j) {
mat[i][j] = (i == j);
}
}
}
Matrix operator*(const Matrix &a) {
Matrix ans;
ans.clear(a.n, m);
for(int i = 1; i <= a.n; ++i) {
for(int j = 1; j <= m; ++j) {
for(int k = 1; k <= n; ++k) {
ans.mat[i][j] += mat[k][j] * a.mat[i][k];
}
}
}
return ans;
}
Matrix operator*=(const Matrix &a) {
return (*this = *this * a);
}
};
Matrix Pow(const Matrix &a, int b) {
Matrix res;
res.Set(a.n);
for(; b; a *= a, b >>= 1) {
if(b & 1) {
res *= a;
}
}
return res;
}
应用
矩阵快速幂可以用来加速递推。比如：
求出斐波那契数列的第 \(N\) 项 \(\bmod (10^9+7)\) 的结果。\((1\le N \le 10^{18})\)
我们可以从这个角度来看：
设斐波那契数列的第 \(i\) 项为 \(f_i\)。如果我们知道了 \(f_{x-1},f_{x}\)，我们怎么求 \(f_x,f_{x+1}\) 呢？
很容易可以得到：
\[\begin{cases} f_x=f_x\\ f_{x+1}=f_{x-1}+f_x \end{cases} \]转换一下：
\[\begin{cases} f_x\ \ \ =0f_{x-1}+1f_x\\ f_{x+1}=1f_{x-1}+1f_x \end{cases} \]我们可以再将式子转变为：
\[\begin{bmatrix}f_{x-1}\\f_x\end{bmatrix}\begin{bmatrix}0\ \ 1\\1\ \ 1\end{bmatrix}=\begin{bmatrix}f_x\\f_{x+1}\end{bmatrix} \]所以我们可以得到斐波那契数列的第 \(N\) 项为矩阵 \(\begin{bmatrix}1\\1\end{bmatrix}\begin{bmatrix}0\ \ 1\\1\ \ 1\end{bmatrix}^{N-1}\) 的 \((1,1)\) 位置的元素。
代码
#include<bits/stdc++.h>
using namespace std;
using ll = long long;
const int MAXN = 3, MOD = int(1e9) + 7;
struct Matrix {
int n, m, mat[MAXN][MAXN];
void clear(int a, int b) {
n = a, m = b;
for(int i = 1; i <= n; ++i) {
for(int j = 1; j <= m; ++j) {
mat[i][j] = 0;
}
}
}
void Set(int a) {
n = m = a;
for(int i = 1; i <= n; ++i) {
for(int j = 1; j <= m; ++j) {
mat[i][j] = (i == j);
}
}
}
Matrix operator*(const Matrix &a) {
Matrix ans;
ans.clear(a.n, m);
for(int i = 1; i <= a.n; ++i) {
for(int j = 1; j <= m; ++j) {
for(int k = 1; k <= n; ++k) {
ans.mat[i][j] = (ans.mat[i][j] + 1ll * mat[k][j] * a.mat[i][k] % MOD) % MOD;
}
}
}
return ans;
}
Matrix operator*=(const Matrix &a) {
return (*this = *this * a);
}
}a, Mat;
ll n;
Matrix Pow(Matrix &a, int b) {
Matrix res;
res.Set(a.n);
for(; b; a *= a, b >>= 1) {
if(b & 1) {
res *= a;
}
}
return res;
}
int main() {
ios::sync_with_stdio(false), cin.tie(0), cout.tie(0);
cin >> n;
a.clear(2, 1), a.mat[1][1] = a.mat[2][1] = 1;
Mat.clear(2, 2), Mat.mat[1][2] = Mat.mat[2][1] = Mat.mat[2][2] = 1;
cout << (a * Pow(Mat, n - 1)).mat[1][1];
return 0;
}