高斯消元和矩阵快速幂

高斯消元

高斯消元是一种能在 \(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;
}

广义矩阵乘法

首先看到以下这个问题：

有 \(N\) 个饼干，第 \(i\) 个饼干的美味值为 \(a_i\)，你可以吃任意不相邻的一些饼干，问美味值之和最大能是多少。

这个问题看起来不能用普通的矩阵乘法解决，所以我们看到：

Max(Min)-plus 矩阵

首先令 \(dp_{i,0/1}\) 表示考虑前 \(i\) 个饼干，最后一个选 / 不选时美味值之和的最大值。然后我们能得到：

\[\begin{cases} dp_{i,0}=\max(dp_{i-1,0},dp_{i-1,1})\\ dp_{i,1}=dp_{i-1,0}+a_i \end{cases} \]

可以发现，这时上面出现了一个 \(\max\)，但下面没有。由于 \(\max\) 不能去掉，所以考虑把下面也变成 \(\max\) 的形式：

\[\begin{cases} dp_{i,0}=\max(dp_{i-1,0},dp_{i-1,1})\\ dp_{i,1}=\max(dp_{i-1,0}+a_i,dp_{i-1,1}-\infty) \end{cases} \]

这里我们把下面的 \(dp_{i-1,1}\) 减去了 \(\infty\)，表示不会从 \(dp_{i-1,1}\) 转移过来。

现在这个式子看起来就很舒服了，我们尝试把它也变成矩阵的形式：

\[\begin{bmatrix} dp_{i,0}\\dp_{i,1} \end{bmatrix}= \begin{bmatrix} dp_{i-1,0}\\dp_{i-1,1} \end{bmatrix} \begin{bmatrix} 0&0\\ a_i&-\infty \end{bmatrix} \]

我们再修改一下矩阵乘法的定义为：\(C=A\times B \Rightarrow C_{i,j}=\min \limits_{k=1}^{n}\{A_{k,j}+B_{i,k}\}\)。

这下我们就得到了 Max(Min)-plus 矩阵的乘法定义。

代码

struct Max_plus_Matrix {
  int n, m;
  ll 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] = -INF;
      }
    }
  }
  Max_plus_Matrix operator*(const Max_plus_Matrix &b) {
    Max_plus_Matrix ret;
    ret.Clear(n, b.m);
    for(int i = 1; i <= n; ++i) {
      for(int j = 1; j <= b.m; ++j) {
        for(int k = 1; k <= b.n; ++k) {
          ret.mat[i][j] = max(ret.mat[i][j], mat[k][i] + b.mat[j][k]);
        }
      }
    }
    return ret;
  }
  Max_plus_Matrix operator*=(const Max_plus_Matrix &b) {
    return *this = *this * b;
  }
};

posted @ 2024-07-02 16:17 Yaosicheng124 阅读(7) 评论(0) 编辑收藏举报

会员力量，点亮园子希望

刷新页面返回顶部

yaosicheng124

这个人很菜，什么也没有留下。

高斯消元和矩阵快速幂

高斯消元

代码

拓展

代码

矩阵快速幂

代码

应用

代码

广义矩阵乘法

Max(Min)-plus 矩阵

代码

公告