一些些数学的算法笔记

好好好，直接进入正题（

---

首先我们先要讲讲矩阵，矩阵你可以理解成 $n\times m$ 的一个二维数组，我们如下表示它：

$$\left[\begin{array}{cccc}{a}_{1,1}& a_{1,2}& \cdots & a_{1,m}\\ a_{2,1}& a_{2,2}& \cdots & a_{2,m}\\ ⋮& ⋮& \ddots & ⋮\\ a_{n,1}& a_{n,2}& \cdots & a_{n,m}\end{array}\right]$$

一、矩阵乘法

1 操作方法

矩阵最重要的一个运算就是矩阵乘法，矩阵乘法就是给定一个 $n\times m$ 的矩阵 $A$，和一个 $m\times l$ 的矩阵 $B$，得到一个 $n\times l$ 的矩阵 $C$，而 $C$ 的每一位 $c_{i,j}=\sum _{k=1}^m{a}_{i,k}\times b_{k,j}$。简单来说就是 $c_{i,j}$ 等于 $A$ 矩阵的第 $i$ 行，与 $B$ 矩阵的第 $j$ 列的元素顺序分别相乘再相加得到的数。

---

矩阵我们一般会把他封装到一个结构体里，接下来看看一道题吧。

LG-P3390 【模板】矩阵快速幂

给定 $n\times n$ 的矩阵 $A$，求 $A^k$。

首先我们要声明一件事，矩阵乘法满足结合律，但不满足交换律。前者可以直接爆算得出（或者用矩阵就是向量的变换得出），后者可以举例得出。我们知道任何满足结合律的运算都可以用快速幂求解，所以我们只要写出矩阵乘法的代码和快速幂的代码即可。

#include <bits/stdc++.h>

using namespace std;

const int kMaxN = 105, kM = 1e9 + 7;

int n;
long long k;

struct M { // 矩阵结构体
  int n, m;
  long long a[kMaxN][kMaxN];
  M(int h, int w) { n = h, m = w, memset(a, 0, sizeof(a)); }
  M() { memset(a, 0, sizeof(a)); } // 这两个都是初始化矩阵
  M operator*(M b) { // 重载乘法运算符，实现矩阵乘法
    M r(n, b.m);
    for (int k = 1; k <= m; k++) {
      for (int i = 1; i <= n; i++) {
        for (int j = 1; j <= b.m; j++) {
          r.a[i][j] = (r.a[i][j] + a[i][k] * b.a[k][j]) % kM;
        }
      }
    }
    return r;
  }
} A;

M fpow(M A, long long b) { // 矩阵快速幂
  M r(n, n);
  for (int i = 1; i <= r.n; i++) {
    r.a[i][i] = 1;
  }
  for (; b; b >>= 1, A = A * A) {
    if (b & 1) r = r * A;
  }
  return r;
}

int main() {
  ios::sync_with_stdio(0), cin.tie(0);
  cin >> n >> k, A = M(n, n);
  for (int i = 1; i <= n; i++) {
    for (int j = 1; j <= n; j++) {
      cin >> A.a[i][j];
    }
  }
  A = fpow(A, k);
  for (int i = 1; i <= n; i++) {
    for (int j = 1; j <= n; j++) {
      cout << A.a[i][j] << ' ';
    }
    cout << '\n';
  }
  return 0;
}

2 例题

介绍完矩阵乘法，你们是不是就要问：“这有什么用呢？”其实这个东西的用处很大，下面再讲一道题。

LG-P1962 斐波那契数列
给定 $n$，求斐波那契的第 $n$ 项模 ${10}^9+7$ 得到的值（$n<2^{63}$）。

我们设斐波那契数列第 $i$ 项的值为 $F_i$，那么 $F_i=F_{i-1}+F_{i-2}=1\times F_{i-1}+1\times F_{i-2}$，考虑构建矩阵 $[F_i{F}_{i+1}]$，若要将其转移成 $[F_{i+1}{F}_{i+2}]$ 即 $[1\times F_{i+1}1\times F_i+1\times F_{i+1}]$，我们可以使用矩阵乘法进行转移。

$$\left[\begin{array}{cc}{F}_i& F_{i+1}\end{array}\right]\times \left[\begin{array}{cc}0& 1\\ 1& 1\end{array}\right]=\left[\begin{array}{cc}{F}_{i+1}& F_{i+2}\end{array}\right]$$

所以我们可以通过多次乘系数矩阵，得到 $F_n$ 的值。由于矩阵乘法满足结合律，所以我们你可以运用矩阵快速幂快速计算 $F_n$。

$$\left[\begin{array}{cc}{F}_1& F_2\end{array}\right]\times {\left[\begin{array}{cc}0& 1\\ 1& 1\end{array}\right]}^{n-1}=\left[\begin{array}{cc}{F}_n& F_{n+1}\end{array}\right]$$

#include <bits/stdc++.h>

using namespace std;

const int kM = 1e9 + 7;

long long n;

struct M {
  int n, m;
  long long a[5][5];
  M(int h, int w) { n = h, m = w, memset(a, 0, sizeof(a)); }
  M() { memset(a, 0, sizeof(a)); }
  M operator*(M b) {
    M r(n, b.m);
    for (int k = 1; k <= m; k++) {
      for (int i = 1; i <= n; i++) {
        for (int j = 1; j <= b.m; j++) {
          r.a[i][j] = (r.a[i][j] + a[i][k] * b.a[k][j]) % kM;
        }
      }
    }
    return r;
  }
} A(1, 2), B(2, 2);

M fpow(M A, long long b) {
  M r(2, 2);
  for (int i = 1; i <= r.n; i++) {
    r.a[i][i] = 1;
  }
  for (; b; b >>= 1, A = A * A) {
    if (b & 1) r = r * A;
  }
  return r;
}

int main() {
  cin >> n;
  A.a[1][1] = A.a[1][2] = B.a[1][2] = B.a[2][1] = B.a[2][2] = 1;
  A = A * fpow(B, n - 1);
  cout << A.a[1][1];
  return 0;
}

看完这道题，我们也可以总结出一些规律，我们能够发现一些简单递推，如果范围较大，可以使用矩阵乘法求解，例如：

$$F_i=a_1\times F_{i-1}+a_2\times F_{i-2}+\cdots +a_m\times F_{i-m}(i<m)$$

像上述递推式可以如下求解：

$$\left[\begin{array}{cccc}{F}_{i-m}& F_{i-m+1}& \cdots & F_i\end{array}\right]\times B=\left[\begin{array}{cccc}{F}_{i-m+1}& F_{i-m+2}& \cdots & F_{i+1}\end{array}\right]$$

$B$ 为系数矩阵，自己求。

---

T234597 路径问题（困难）
有一个 $n\times m$ 的网格，初始时有个人在 $(1,1)$ 处，他每次只能向右走一步或者向下走一步，而且不能走出网格。
请问到达 $(n,m)$ 有多少种方案，请输出方案数模 ${10}^9+7$ 的结果。（$1\le n\le {10}^{18}$，$1\le m\le 100$）

设 $F_{i,j}$ 为到达 $(i,j)$ 的方案数，我们很容易推出转移方程：

$$F_{i,j}=\{\begin{array}{l}1,(i=1,j=1)\\ F_{i,j-1},(i=1,j>1)\\ F_{i-1,j},(i>1,j=1)\\ F_{i,j-1}+F_{i-1,j},(i>1,j>1)\end{array}$$

由于这样的转移方程牵扯到前一行的状态，所以我们将每一行的 $F$ 作为矩阵表示出来，做相应的矩阵乘法即可。

---

P1397 [NOI2013] 矩阵游戏

$F[i,j]$ 满足下面的递推式:

$$\begin{array}{rl}{F}_{1,1}& =1\\ F_{i,j}& =a\times F_{i,j-1}+b,& j\ne 1\\ F_{i,1}& =c\times F_{i-1,m}+d,& i\ne 1\end{array}$$

求 $F_{i,j}$ 除以 ${10}^9+7$ 的余数 （$1\le n,m\le {10}^{1000000}$；$1\le a,b,c,d\le {10}^9$）

对于上述递推式，我们先考虑每一行，每一行的变化，就是将上一个数 $\times a+b$。考虑初始矩阵为 $[1F_{i,j}]$，若要将其变为 $[1a\times F_{i,j}+b]$，转移方程为：

$$\left[\begin{array}{cc}1& F_{i,j}\end{array}\right]\times \left[\begin{array}{cc}1& b\\ 0& a\end{array}\right]=\left[\begin{array}{cc}1& F_{i,j+1}\end{array}\right]$$

那么考虑列的变化的话就是将系数矩阵的 $b$ 改成 $d$，$a$ 改成 $c$，即可。

设

$$A=\left[\begin{array}{cc}1& 1\end{array}\right],B=\left[\begin{array}{cc}1& b\\ 0& a\end{array}\right],C=\left[\begin{array}{cc}1& d\\ 0& c\end{array}\right]$$

$A$ 为初始矩阵，$B,C$ 为系数矩阵。

那么

$$\left[\begin{array}{cc}1& F_{n,m}\end{array}\right]=A\times {(B^{m-1}\times C)}^{n-1}\times B^{m-1}$$

由于 $n,m$ 太大，$n-1,m-1$ 需要用到高精度，快速幂也需要用到 $10$ 进制快速幂即可。

#include <bits/stdc++.h>

using namespace std;

const int kM = 1e9 + 7;

string n, m;

string J(string s) {
  for (int i = s.size() - 1; ~i; i--) {
    if (s[i] == '0') {
      s[i] = '9';
    } else {
      s[i]--;
      break;
    }
  }
  return s;
}

struct M {
  int n, m;
  long long a[5][5];
  M(int h, int w) { n = h, m = w, memset(a, 0, sizeof(a)); }
  M() { memset(a, 0, sizeof(a)); }
  M operator*(M b) {
    M r(n, b.m);
    for (int k = 1; k <= m; k++) {
      for (int i = 1; i <= n; i++) {
        for (int j = 1; j <= b.m; j++) {
          r.a[i][j] = (r.a[i][j] + a[i][k] * b.a[k][j]) % kM;
        }
      }
    }
    return r;
  }
};

M fpow(M A, string b) {
  M r(A.n, A.m), T;
  for (int i = 1; i <= r.n; i++) {
    r.a[i][i] = 1;
  }
  for (int i = b.size() - 1; ~i; i--) {
    T = M(A.n, A.m);
    if (b[i] == '9') {
      T = A * A, T = T * T, T = T * T, r = r * T * A;
    } else {
      for (int j = 1; j <= b[i] - '0'; j++) {
        r = r * A;
      }
    }
    A = A * A, T = A, A = A * A, A = A * A, A = A * T;
  }
  return r;
}

int main() {
  ios::sync_with_stdio(0), cin.tie(0);
  M A(1, 2), B(2, 2), C(2, 2);
  cin >> n >> m >> B.a[2][2] >> B.a[1][2] >> C.a[2][2] >> C.a[1][2];
  A.a[1][1] = A.a[1][2] = B.a[1][1] = C.a[1][1] = 1;
  A = A * fpow(fpow(B, J(m)) * C, J(n)) * fpow(B, J(m));
  cout << A.a[1][2];
  return 0;
}

---

二、高斯消元

1 操作方法

学习完矩阵乘法，我们可以发现一个东西，任何一个 $n$ 元一次方程都可以写成矩阵乘法的形式：

$$\{\begin{array}{l}{a}_{1,1}{x}_1+a_{1,2}{x}_2+\cdots +a_{1,n}=b_1\\ a_{2,1}{x}_1+a_{2,2}{x}_2+\cdots +a_{2,n}=b_2\\ \cdots \\ a_{n,1}{x}_1+a_{n,2}{x}_2+\cdots +a_{n,n}=b_n\end{array}$$

可以写成（上述方程为有 $n$ 个方程的 $n$ 原一次方程）：

$$\left[\begin{array}{cccc}{a}_{1,1}& a_{1,2}& \cdots & a_{1,n}\\ a_{2,1}& a_{2,2}& \cdots & a_{2,n}\\ ⋮& ⋮& \ddots & ⋮\\ a_{n,1}& a_{n,2}& \cdots & a_{n,n}\end{array}\right]\times \left[\begin{array}{c}{x}_1\\ x_2\\ ⋮\\ x_n\end{array}\right]=\left[\begin{array}{c}{b}_1\\ b_2\\ ⋮\\ b_n\end{array}\right]$$

是不是很神奇，虽然上面的内容和我们下面要讲的东西没有太大关系。（

我们在数学中时常遇到解这种线性方程组，但是元数不高，方程数不多，如果我们要解上述方程该怎么办呢？接下来我们用矩阵（分块矩阵）的思想解决这个问题。

我们设

$$A=\left[\begin{array}{cccc}{a}_{1,1}& a_{1,2}& \cdots & a_{1,n}\\ a_{2,1}& a_{2,2}& \cdots & a_{2,n}\\ ⋮& ⋮& \ddots & ⋮\\ a_{n,1}& a_{n,2}& \cdots & a_{n,n}\end{array}\right],X=\left[\begin{array}{c}{x}_1\\ x_2\\ ⋮\\ x_n\end{array}\right],B=\left[\begin{array}{c}{b}_1\\ b_2\\ ⋮\\ b_n\end{array}\right]$$

考虑将 $A$ 和 $B$ 塞进一个矩阵得到 $[AB]$，接下来引入矩阵的一种变换，它不会改变矩阵的本质，叫做初等行变换。初等行变化有三种变换方式：

• 将矩阵的任意一行中的每一个数乘常数 $c(c\ne 0)$。
• 将矩阵的任意两行交换位置。
• 将矩阵的任意一行中的每一个数乘常数 $c(c\ne 0)$，在将这一行的每一个数加到另一行对应的数上。

我们发现这三种操作就是方程的性质。

接下来我们将之前的矩阵 $[AB]$ 进行初等行变换，考虑我们之前学过的解多元一次方程组是如何解的，首先要消元！

我们考虑先消第一个元，如果 $a_{1,1}=0$ 我们就把下面的一行换上来即可，如果下面全都是 $0$，那么就是无解或者无数多组解（后面会提到如何判定这两种情况），这样操作完后，我们就可以将第一行下面的所有行的第一列都消成 $0$，然后看第二行，找到第二列，跟第一行做初等行变换消元，然后继续如此操作，继续消，以此类推。不出意外最后我们可能应该也许会消成一个阶梯矩阵，即：

$$\left[\begin{array}{ccccc}{a}_{1,1}^{{}^{\prime }}& a_{1,2}^{{}^{\prime }}& \cdots & a_{1,n}^{{}^{\prime }}& b_1^{{}^{\prime }}\\ 0& a_{2,2}^{{}^{\prime }}& \cdots & a_{2,n}^{{}^{\prime }}& b_2^{{}^{\prime }}\\ ⋮& ⋮& \ddots & ⋮& ⋮\\ 0& 0& \cdots & a_{n,n}^{{}^{\prime }}& b_n^{{}^{\prime }}\end{array}\right]$$

得到阶梯矩阵，我们就可以从第 $n$ 层开始消，将上一行得到的值代入然后解出所有的未知数。

或者说可以用倒数第二行与倒数第一行进行消元，倒数第三行与倒数第一行和倒数第二行进行消元以此类推，最终可以变为：

$$\left[\begin{array}{ccccc}{a}_{1,1}^{{}^{″}}& 0& \cdots & 0& b_1^{{}^{″}}\\ 0& a_{2,2}^{{}^{″}}& \cdots & 0& b_2^{{}^{″}}\\ ⋮& ⋮& \ddots & ⋮& ⋮\\ 0& 0& \cdots & a_{n,n}^{{}^{″}}& b_n^{{}^{″}}\end{array}\right]$$

这样就可以解得每一个 $x$ 啦。