数学知识(二)
温馨提示:本章内容涉及线性代数知识较多!
一. 矩阵
1. 概念
矩阵\(\texttt{(Matrix)}\)是一个按照长方阵列排列的复数或实数集合,最早来自于方程组的系数及常数所构成的方阵。
说白了,一个 \(n\times m\) 的矩阵可以看作一个 \(n\times m\) 的二维数组,比如一个 \(2\times 2\) 的矩阵写作:
2. 矩阵的运算
(1) 加减法
矩阵的加减法就是把两个矩阵对位相加减。
运算规定:只有行数和列数都分别相等的矩阵才能进行加减运算!
比如:
那么令 \(C = A + B\),则:
若 \(C = A - B\),则:
具有的性质:交换律、结合律。
(2) 乘法
大名鼎鼎的矩阵乘法当然要复杂亿点点!
运算规定:只有行数和列数都分别相等的矩阵才能进行加减运算!
设 \(A\) 是一个 \(x\times y\) 的矩阵,\(B\) 是一个 \(y\times z\) 的矩阵,则 \(C = A \times B\) 是一个 \(x\times z\) 的矩阵,且满足 \(\forall i\in [1, x],\forall j\in [1, z]:\)
简单来说,\(A\) 矩阵的第 \(i\) 行乘上 \(B\) 矩阵的第 \(j\) 列就得到 \(C\) 矩阵的第 \(i\) 行第 \(j\) 列。
具有的性质:结合律。
单位矩阵:
只有对角线全是 \(1\),其他位置全是 \(0\) 的矩阵。
性质:任何矩阵 \(A\) 乘上单位矩阵都等于 \(A\) 自己。
3. 矩阵加速递推
终于迎来了矩阵乘法的第一个用途:矩阵加速递推!
在此之前,还得先了解一个叫做 “矩阵快速幂” 的东西。
其实本质上和整数的快速幂差不多,矩阵快速幂可以快速地将一个矩阵自乘多次,它将指数 \(k\) 拆分成 \(2\) 的次幂,从而能在 \(O(\log k)\) 的时间复杂度内解决问题。
P3390 【模板】矩阵快速幂
\(\texttt{Code:}\)
#include <cstring>
#include <iostream>
using namespace std;
const int N = 110, mod = 1e9 + 7;
typedef long long ll;
int n;
ll k;
int A[N][N];
//矩阵乘上一个矩阵
void Mul(int c[][N], int a[][N], int b[][N], int size) {
int tmp[N][N] = {0};
for(int i = 0; i < size; i++)
for(int j = 0; j < size; j++)
for(int k = 0; k < size; k++)
tmp[i][j] = (tmp[i][j] + (ll)a[i][k] * b[k][j]) % mod;
//注意 c 是一个指针,所以只能写成 sizeof tmp
memcpy(c, tmp, sizeof tmp);
}
int main() {
scanf("%d%lld", &n, &k);
for(int i = 0; i < n; i++)
for(int j = 0; j < n; j++)
scanf("%d", &A[i][j]);
int res[N][N];
//初始化一个单位矩阵
for(int i = 0; i < n; i++) res[i][i] = 1;
while(k) {
if(k & 1) Mul(res, res, A, n);
Mul(A, A, A, n);
k >>= 1;
}
for(int i = 0; i < n; i++) {
for(int j = 0; j < n; j++)
printf("%d ", res[i][j]);
puts("");
}
return 0;
}
知道了矩阵快速幂,那么它和加速递推有什么关系呢?
比如这道题:
P1962 斐波那契数列
让你求斐波那契数列的第 \(n\) 项,但是 \(n\le 2^{63}\)!
如果根据斐波那契数列的递推公式直接递推是 \(O(n)\) 的,无法通过此题。
不过可以发现,斐波那契数列的第 \(n\) 项只与第 \(n - 1\) 项和第 \(n - 2\) 项有关,所以只用保留最近的两个斐波那契数构造矩阵。
设斐波那契数列的第 \(i\) 项为 \(fib(i)\), \(F(n)\) 表示一个 \(1\times 2\) 的矩阵,\(F(n) = \begin{bmatrix}fib(n) & fib(n + 1)\end{bmatrix}\)。
现在我们希望能根据这个矩阵推出矩阵 \(F(n + 1) = \begin{bmatrix}fib(n + 1) & fib(n + 2)\end{bmatrix}\)。
设矩阵 \(A\) 为变换矩阵,则我们要构造一个矩阵 \(A\) 使 \(F(n)\times A = F(n + 1)\)。
再根据矩阵乘法的定义,我们构造的矩阵应该长这样:
斐波那契数列的每个相邻两项都能靠这种关系递推,所以 \(F(n) = F(0)\times A^n\),其中初始值 \(F(0) = \begin{bmatrix}0 & 1\end{bmatrix}\)。
由于矩阵乘法具有结合律,所以可以先通过矩阵快速幂算出 \(A^n\),再用 \(F(0)\) 左乘它,结果矩阵的第一列即为 \(fib(n)\)。
在矩阵乘法时顺便取模即可。
\(\texttt{Code:}\)
#include <cstring>
#include <iostream>
using namespace std;
const int N = 2, mod = 1e9 + 7;
typedef long long ll;
ll n;
//行向量乘上一个矩阵
void mul(int c[], int a[], int b[][N]) {
int tmp[N] = {0};
for(int i = 0; i < N; i++)
for(int j = 0; j < N; j++)
tmp[i] = (tmp[i] + (ll)a[j] * b[j][i]) % mod;
memcpy(c, tmp, sizeof tmp);
}
void Mul(int c[][N], int a[][N], int b[][N]) {
int tmp[N][N] = {0};
for(int i = 0; i < N; i++)
for(int j = 0; j < N; j++)
for(int k = 0; k < N; k++)
tmp[i][j] = (tmp[i][j] + (ll)a[i][k] * b[k][j]) % mod;
memcpy(c, tmp, sizeof tmp);
}
int main() {
scanf("%lld", &n);
int res[N] = {0, 1};
int A[N][N] = {
{0, 1},
{1, 1}
};
//矩阵快速幂
while(n) {
if(n & 1) mul(res, res, A);
Mul(A, A, A);
n >>= 1;
}
printf("%d\n", res[0]);
return 0;
}
习题
二. 高斯消元
1. 简介
高斯消元是一种用来求解线性方程组的方法,通俗说来就是我们小学学过的加减消元。
设一个线性方程组有 \(m\) 个未知数,\(n\) 个一次方程,那么它的所有系数可以写成一个 \(m\) 行 \(n\) 列的系数矩阵,再加上每个方程等号右侧的常数,可以写成一个 \(m\) 行 \(n + 1\) 列的增广矩阵。
比如当 \(m = n = 3\) 时:
\(\begin{cases} x_1 + 2x_2 - x_3 = -6\\ 2x_1 + x_2 - 3x_3 = -9\\ -x_1 - x_2 + 2x_3 = 7 \end{cases}\Rightarrow \left[\begin{array}{ccc | c} 1 & 2 & -1 & -6\\ 2 & 1 & -3 & -9\\ -1 & -1 & 2 & 7 \end{array}\right] \)
接下来需要用到一点线性代数的基础知识。
我们对这个增广矩阵可以进行一些操作,而原方程组的解不会发生改变。
- 用一个非零的数乘上某一行;
- 把其中一行的若干倍加到另一行上;
- 交换两行的位置。
以上三种操作称为初等行变换。
接下来我们的目标是用初等行变换,将原矩阵变为一个阶梯形矩阵,它的系数矩阵部分被称为上三角矩阵,比如:
\( \left[\begin{array}{ccc | c} 1 & 2 & -1 & -6\\ 2 & 1 & -3 & -9\\ -1 & -1 & 2 & 7 \end{array}\right]\Rightarrow \left[\begin{array}{ccc | c} 1 & 2 & -1 & -6\\ 0 & -3 & -1 & 3\\ -1 & -1 & 2 & 7 \end{array}\right]\Rightarrow \left[\begin{array}{ccc | c} 1 & 2 & -1 & -6\\ 2 & 1 & -3 & -9\\ 0 & 1 & 1 & 1 \end{array}\right] \\\qquad \\\Rightarrow\left[\begin{array}{ccc | c} 1 & 2 & -1 & -6\\ 0 & 1 & 1 & 1\\ 0 & -3 & -1 & 3 \end{array}\right]\Rightarrow \left[\begin{array}{ccc | c} 1 & 2 & -1 & -6\\ 0 & 1 & 1 & 1\\ 0 & 0 & 2 & 6 \end{array}\right]\Rightarrow \left[\begin{array}{ccc | c} 1 & 2 & -1 & -6\\ 0 & 1 & 1 & 1\\ 0 & 0 & 1 & 3 \end{array}\right] \)
最后就得到了一个上三角矩阵,即:
\(\left[\begin{array}{ccc | c} 1 & 2 & -1 & -6\\ 0 & 1 & 1 & 1\\ 0 & 0 & 1 & 3 \end{array}\right]\Rightarrow \left\{\begin{aligned} x_1 + 2x_2 + x_3 = -6 \\ x_2 + x_3 = 1 \\ x_3 = 3 \end{aligned}\right. \)
该矩阵还能继续化简:
\( \left[\begin{array}{ccc | c} 1 & 2 & -1 & -6\\ 0 & 1 & 1 & 1\\ 0 & 0 & 1 & 3 \end{array}\right]\Rightarrow \left[\begin{array}{ccc | c} 1 & 2 & 0 & -3\\ 0 & 1 & 0 & -2\\ 0 & 0 & 1 & 3 \end{array}\right]\Rightarrow \left[\begin{array}{ccc | c} 1 & 0 & 0 & 1\\ 0 & 1 & 0 & -2\\ 0 & 0 & 1 & 3 \end{array}\right] \)
最后得到的矩阵称为简化阶梯形矩阵,事实上,每一个标准的上三角矩阵都能化简成这样,它实际上已经表示出了每个未知数的解。
高斯消元法就是有目的地做这样一件事,它的基本步骤为:
-
找到当前列绝对值最大的一行。
-
用初等行变换 \(3\) 把这一行换到最上面(未确定阶梯型的行,并不是第一行)。
-
用初等行变换 \(1\) 将该行的第一个数变成 \(1\)(其余所有的数字依次跟着变化)。
-
用初等行变换 \(2\) 将下面所有行的当且列的值变成 \(0\)。
\(\texttt{Code:}\)
int gauss() {
int c, r; //c 表示列,r 表示行
for(c = 0, r = 0; c < n; c++) {
int t = r; //从这一行向下找
for(int j = r + 1; j < n; j++)
if(fabs(a[j][c]) > fabs(a[t][c]))
t = j;
if(fabs(a[t][c]) < eps) continue;
if(t != r) for(int i = c; i <= n; i++) swap(a[t][i], a[r][i]); //交换到最上面一行
for(int i = n; i >= c; i--) a[r][i] /= a[r][c]; //把这一行都消一遍
for(int i = r + 1; i < n; i++)
if(fabs(a[i][c]) > eps) //将下面所有非零行都消一遍
for(int j = n; j >= c; j--)
a[i][j] -= a[i][c] * a[r][j];
++r;
}
if(r < n) {
for(int i = r; i < n; i++)
//若消完后出现了非零等于零的情况则一定无解
if(fabs(a[i][n]) > eps)
return -1;
//否则有无穷多组解
return 0;
}
//有解时就化到最简,求出解
for(int i = n - 2; ~i; i--)
for(int j = i + 1; j < n; j++)
a[i][n] -= a[i][j] * a[j][n];
return 1;
}
P3389 【模板】高斯消元法
P2455 [SDOI2006] 线性方程组
2. 例题
P4035 [JSOI2008] 球形空间产生器
题目大意:
给定一个 \(n\) 维球上 \(n + 1\) 个点的坐标,求出该球的球心坐标。
思路:
先将题目条件转化一下:一个球的球心到球面上的任意点的距离都相等。利用这一性质,设球心坐标为 \((x_0, x_1,\cdots, x_{n - 1})\),给定的第 \(i\) 个点(下标从 \(0\) 开始)的坐标为 \((a_{i, 0}, a_{i, 1}\cdots, a_{i, n - 1})\),球的半径为 \(r\),那么就有:
但是这样是不能进行高斯消元的,因为 \(r\) 是未知数,而且方程中还有很多平方,这都是高斯消元不能解决的。
一共有 \(n + 1\) 个方程,\(n + 1\) 个未知数且没有无用方程,所以原方程一定能解出来,根据这一思想,我们先把每一个方程左式的所有平方展开,得:
发现所有方程的左式前面都是平方和,右边都是 \(r^2\),所以可以使用平方差公式,差分处理,具体地,我们用第 \(2\) 行减去第 \(1\) 行得到新方程组的第 \(1\) 行,用第 \(3\) 行减去第 \(2\) 行得到新方程组的第 \(2\) 行……以此类推,构造出一个 \(n\) 行,含有 \(n\) 个未知数的方程组(设 \(a_{i, 0}^2 + a_{i, 1}^2 + \cdots + a_{i, n - 1}^2) = sqsum_i\)):
这时候就可以快乐地高斯消元了。
\(\texttt{Code:}\)
#include <cmath>
#include <iostream>
using namespace std;
const int N = 15;
const double eps = 1e-9;
int n;
double a[N][N];
void gauss() {
int c, r;
for(c = 0, r = 0; c < n; c++) {
int t = r;
for(int i = r + 1; i < n; i++)
if(fabs(a[i][c]) > fabs(a[t][c]))
t = i;
if(fabs(a[t][c]) < eps) continue;
if(t != r) for(int i = c; i <= n; i++) swap(a[t][i], a[r][i]);
for(int i = n; i >= c; i--) a[r][i] /= a[r][c];
for(int i = r + 1; i < n; i++)
if(fabs(a[i][c]) > eps)
for(int j = n; j >= c; j--)
a[i][j] -= a[i][c] * a[r][j];
++r;
}
for(int i = n - 2; ~i; i--)
for(int j = i + 1; j < n; j++)
a[i][n] -= a[i][j] * a[j][n];
}
int main() {
scanf("%d", &n);
for(int i = 0; i <= n; i++)
for(int j = 0; j < n; j++)
scanf("%lf", &a[i][j]);
for(int i = 0; i < n; i++) {
double sqsum1 = 0, sqsum2 = 0;
for(int j = 0; j < n; j++) {
sqsum1 += a[i][j] * a[i][j];
sqsum2 += a[i + 1][j] * a[i + 1][j];
a[i][j] = 2 * (a[i + 1][j] - a[i][j]);
}
a[i][n] = sqsum2 - sqsum1;
}
gauss();
for(int i = 0; i < n; i++)
printf("%.3lf ", a[i][n]);
return 0;
}
P10315 [SHUPC 2024] 原神,启动!
题目大意:
有 \(n\) 个石碑,每个石碑有 \(0\sim m - 1\) 共 \(m\) 种状态,击打一个石碑会带动其他的石碑。若当前石碑的状态是 \(s\),则击打或被带动后的状态为 \((s + 1)\bmod m\)。
现给定这 \(n\) 个石碑的初始状态 \(s_i\)、每个石碑带动的石碑及末状态 \(t_i\),求每个石碑至少被击打几次。
思路:
首先把题面意思抽象出来,令 \(a_{i, j}\) 表示第 \(i\) 个石碑和第 \(j\) 个石碑的联系,若 \(a_{i, j} = 1\),则表示击打 \(j\) 会带动 \(i\),若 \(a_{i, j} = 0\) 表示无影响,特别的,因为击打自己就相当于带动自己,所以 \(a_{i, i} = 1\)。
再令 \(x_i\) 表示第 \(i\) 个石碑被击打的次数。
那么题面就变为了:
再移个项,得:
求出每个 \(x_i\) 的最小非负整数解即可。
同余只是纸老虎!直接转换成等号,高斯消元求解即可,只是需要把解映射到 \([0, m - 1]\)。
同时这道题还需要在无穷多组解时输出任意一组解,需要在消元时额外注意(其实应该只有我这种写法应该注意),不能直接回代求解。
当找到一个主行 \(r\) 时,不要只从 \(r + 1\) 消到 \(n\),而应该把 \(1\sim n\) 都消一遍,此时除每行的首变量(每个行向量中第一个系数非零的未知数)之外其他的都是自由元,直接将首变量赋值,自由元赋成 \(0\) 不管就行了。
这是我原来的高斯消元代码:
int gauss() {
int c, r;
for(c = 0, r = 0; c < n; c++) {
int t = r;
for(int i = r + 1; i < n; i++)
if(abs(a[i][c]) > abs(a[t][c]))
t = i;
if(!a[t][c]) continue;
if(t != r) for(int i = c; i <= n; i++) swap(a[t][i], a[r][i]);
for(int i = n; i >= c; i--) a[r][i] = (a[r][i] * (qpow(a[r][c], mod - 2) + mod) % mod) % mod;
for(int i = r + 1; i < n; i++) //原来是消第 r + 1 到 n 行
if(a[i][c])
for(int j = n; j >= c; j--)
a[i][j] = (mod + a[i][j] - a[i][c] * a[r][j] % mod) % mod;
++r;
}
if(r < n) {
for(int i = r; i < n; i++)
if(a[i][n] > 0)
return -1;
}
for(int i = n - 2; ~i; i--)
for(int j = i + 1; j < n; j++)
a[i][n] = (mod + a[i][n] - a[i][j] * a[j][n] % mod) % mod;
return 1;
}
它在这组数据时会出错:
3 3
2 2 3
2 1 3
1 1
0 0 0
2 1 2
Answer:
0 1 1
或:
1 0 1
My answer:
1 1 0
原矩阵:
这是因为用以上代码消出来的结果为:
而如果直接回代就会直接将 \(x_3\) 钦定为 \(0\),这是不对的,因为 \(x_2\) 才是自由元,而 \(x_3\) 有固定的解 \(1\)。
而采用全部重消一遍的方法就能保证所有首变量都只会在一个行向量中出现,这时候回代就完全不用考虑和其他首变量取值出现冲突的问题。
\(\texttt{Code:}\)
#include <cmath>
#include <iostream>
using namespace std;
const int N = 110;
typedef long long ll;
int n, mod;
ll a[N][N];
ll ans[N];
ll qpow(ll a, int b) {
ll ans = 1, base = a % mod;
while(b) {
if(b & 1) ans = ans * base % mod;
base = base * base % mod;
b >>= 1;
}
return ans;
}
void output() {
puts("---------");
for(int i = 0; i < n; i++) {
for(int j = 0; j <= n; j++)
printf("%d ", a[i][j]);
puts("");
}
puts("---------");
}
int gauss() {
int c, r;
for(c = 0, r = 0; c < n; c++) {
int t = r;
for(int i = r + 1; i < n; i++)
if(abs(a[i][c]) > abs(a[t][c]))
t = i;
if(!a[t][c]) continue;
if(t != r) for(int i = c; i <= n; i++) swap(a[t][i], a[r][i]);
for(int i = n; i >= c; i--) a[r][i] = (a[r][i] * (qpow(a[r][c], mod - 2) + mod) % mod) % mod;
for(int i = 0; i < n; i++) //全部重消一遍
if(a[i][c] && i != r)
for(int j = n; j >= c; j--)
a[i][j] = (mod + a[i][j] - a[i][c] * a[r][j] % mod) % mod; //注意取模时要加上模数以防负数
++r;
}
if(r < n) {
for(int i = r; i < n; i++)
if(a[i][n] > 0)
return -1;
//
for(int i = 0; i < n; i++) {
for(int j = 0; j < n; j++)
if(a[i][j]) { //寻找首变量
ans[j] = a[i][n]; //直接将首变量赋成方程右侧的值,其他变量都是自由元,直接赋成 0
break;
}
}
return 0;
}
for(int i = n - 2; ~i; i--)
for(int j = i + 1; j < n; j++)
a[i][n] = (mod + a[i][n] - a[i][j] * a[j][n] % mod) % mod;
for(int i = 0; i < n; i++)
ans[i] = a[i][n];
return 1;
}
int main() {
scanf("%d%d", &n, &mod);
ll x;
for(int i = 0; i < n; i++) {
int cnt;
scanf("%d", &cnt);
while(cnt--) {
scanf("%lld", &x);
a[x - 1][i] = 1;
}
a[i][i] = 1;
}
for(int i = 0; i < n; i++)
scanf("%lld", &a[i][n]);
for(int i = 0; i < n; i++) {
scanf("%lld", &x);
a[i][n] = x - a[i][n];
}
int type = gauss();
if(type >= 0)
for(int i = 0; i < n; i++)
printf("%lld ", ans[i]);
else puts("niuza");
return 0;
}
在此鸣谢大佬 Lu_xZ 的指导!
P10499 开关问题
题目大意:
有 \(n\) 个开关,\(0\) 表示关,\(1\) 表示开,每个开关还有带动的开关,若操作一个开关,那么它带动的开关也会相应变换。
现给出这 \(n\) 个开关的初始状态 \(s_i\) 和末状态 \(t_i\),询问有多少种方法能将初始态转变为末态(不考虑操作先后顺序且每个开关至多操作一次)。
思路:
高斯消元解异或方程组经典题。
先考虑将原题抽象成方程组。
令 \(x_i\) 表示第 \(i\) 个开关的操作次数,因为每个开关至多操作一次,所以 \(x_i = 0\) 或 \(x_i = 1\)。
令 \(a_{i, j}\) 表示第 \(i\) 个开关和第 \(j\) 个开关之间的联系,若 \(a_{i, j} = 1\),则表示操作 \(j\) 会带动 \(i\),若 \(a_{i, j} = 0\) 表示无影响,特别的,因为操作自己就相当于带动自己,所以 \(a_{i, i} = 1\)。
再根据操作效果:\(0\) 变 \(1\),\(1\) 变 \(0\),和异或一模一样。
所以可以列出以下方程组:
异或其实就是不进位加法,所以也可以用高斯消元来解,将加减法换为异或就行了。
这道题要求操作方案数,那么找自由元的数量就好了。因为若某个未知数是自由元,那么它取 \(0\) 和 \(1\) 都可以,于是贡献了两种方案,根据乘法原理,应该把自由元的数量这么多 \(2\) 乘起来,即 \(2^{cnt}\),\(cnt\) 为自由元的数量。
同时由于系数只能为 \(0\) 或 \(1\),所以一个行向量可以压缩为一个二进制整数或者用 bitset 来操作,这样就能一次异或一整行,时间复杂度降低为 \(O(\frac{n^3}{\omega})\),写起来也方便许多。
\(\texttt{Code:}\)
#include <cmath>
#include <bitset>
#include <iostream>
using namespace std;
const int N = 35;
int T;
int n;
bitset<N> a[N];
int ans;
int gauss() {
int c, r;
for(c = 0, r = 0; c < n; c++) {
int t = r;
for(int i = r + 1; i < n; i++)
if(a[i][c]) {
t = i;
break;
}
if(!a[t][c]) continue;
if(t != r) swap(a[t], a[r]);
for(int i = r + 1; i < n; i++)
if(a[i][c])
a[i] = a[i] ^ a[r];
++r;
}
if(r < n) {
for(int i = r; i < n; i++) {
if(a[i][n])
return -1;
}
ans = 1 << n - r;
return 0;
}
return 1;
}
int main() {
ios::sync_with_stdio(false);
cin.tie(nullptr);
cout.tie(nullptr);
cin >> T;
while(T--) {
cin >> n;
ans = 1;
int x;
for(int i = 0; i < n; i++)
a[i].reset(), a[i].set(i, 1);
for(int i = 0; i < n; i++) {
cin >> x;
a[i][n] = x;
}
for(int i = 0; i < n; i++) {
cin >> x;
a[i][n] = a[i][n] ^ x;
}
int y;
while(cin >> x >> y && x && y)
a[y - 1].set(x - 1, 1);
int type = gauss();
if(type >= 0) cout << ans << '\n';
else cout << "Oh,it's impossible~!!" << '\n';
}
return 0;
}