线性代数

线性代数在 OI 中也有非常多的应用,简单地介绍一下。

3b1b 线性代数的本质

矩阵加速

一、例题

P3216 [HNOI2011]数学作业

显然这样的数有递推式 fi=fi×10k+i(k 为 i 的位数)f_i= f_i \times10^k + i (k \ \text{为} \ i \ \text{的位数})i1018i\leq 10^{18},显然要用矩阵加速优化递推,注意要分段。

注意! 重载运算符的话一定要注意前后顺序!否则的话会出错。

矩阵行列式

一、表示

  • 一个 n×nn\times n 的正方形矩阵 AA 的行列式记为 A|A|detA\det A,一个 2×22\times 2 矩阵的行列式可表示如下:

det(abcd)=abcd=adbc\det \begin{pmatrix} a & b \\ c & d\end{pmatrix}= \begin{vmatrix} a & b \\ c & d \end{vmatrix} = ad-bc

  • 一个 n×nn\times n 矩阵的行列式等于其任意行(或列)的元素与对应的代数与余子式乘积之和,即:

det(A)=ai1Ai1++ainAin=j=1nai,j(1)i+jdet(Ai,j)\det(A)=a_{i1}A_{i1}+\cdots+a_{in}A_{in}=\sum_{j=1}^na_{i,j}(-1)^{i+j}\det(A_{i,j})

二、性质

  • 矩阵某两行(列)交换,行列式 ×(1)\times (-1)

  • 矩阵某一行(列)加上另一行(列),行列式不变。

  • 矩阵某一行(列) ×k\times k,行列式 ×k\times k

三、代码

P7112 【模板】行列式求值

#include <bits/stdc++.h>
using namespace std;
typedef long long ll;

inline ll read() {
    ll sum = 0, ff = 1;
    char ch = getchar();
    while(ch < '0' || ch > '9') {
        if(ch == '-') {
            ff = -1;
        }
        ch = getchar();
    }
    while(ch >= '0' && ch <= '9') {
    	sum = sum * 10 + ch - 48;
        ch = getchar();
    }
    return sum * ff;
}

void write(ll x) {
	if(x < 0)
	    putchar('-'), x = -x;
	if(x > 9)
		write(x / 10);
	putchar(x % 10 + '0');
}

const int N = 607;
int n, mod, a[N][N], tmp, ans = 1, opt = 1;

void sol() {
	for(int i = 1; i <= n; i++)
		for(int j = i + 1; j <= n; j++) {
			while(a[i][i]) {//辗转相减法 
				tmp = a[j][i] / a[i][i];
				for(int k = i; k <= n; k++)
					a[j][k] = (a[j][k] - (ll)a[i][k] * tmp % mod + mod) % mod;
				swap(a[i], a[j]), opt = -opt;
			}
			swap(a[i], a[j]), opt = -opt;
		}	
	for(int i = 1; i <= n; i++)
		ans = (ll)ans * a[i][i] % mod;
	ans *= opt;
}

int main() {
//  	freopen(".in", "r", stdin);
//  	freopen(".out", "w", stdout);
	n = read(), mod = read();
	for(int i = 1; i <= n; i++)
		for(int j = 1; j <= n; j++)
			a[i][j] = read();//, a[i][j] %= mod;
	sol();
	write((ans + mod) % mod);
    return 0;
}

高斯消元

一、简介

高斯消元是求解线性方程组、行列式计算等的经典算法,在很多方面都有涉及,因此有必要学一下。

二、增广矩阵

有方程组

2x+yz=82x+y-z=8

3xy+2z=11-3x-y+2z=-11

2x+y+2z=3-2x+y+2z=-3

写成矩阵的形式为

(2113122128113)\left( \begin{matrix} 2 & 1 & -1 \\ -3 & -1 & 2 \\ -2 & 1 & 2\end{matrix} \middle| \begin{matrix} 8 \\ -11 \\ -3 \end{matrix} \right)

这就是 增广矩阵,即方程组系数矩阵 AA 与常数列 bb 的并发生的新矩阵,表示为 (Ab)(A|b)

三、操作

从整体上来说,我们从上到下依次处理每一行,处理完第 ii 行后,让 ai,ia_{i,i}00,而 aj,i (j>i)a_{j,i} \ (j > i) 均为 00,最后让 ai,ia_{i,i}11,过程如下:

(2113122128113)(31201/31/305/32/3112/313/3)(31205/32/3001/51113/31/5)(11/32/3012/500111/313/51)\left( \begin{matrix} 2 & 1 & -1 \\ -3 & -1 & 2 \\ -2 & 1 & 2\end{matrix} \middle| \begin{matrix} 8 \\ -11 \\ -3 \end{matrix} \right) \Rightarrow \left( \begin{matrix} -3 & -1 & 2 \\ 0 & 1/3 & 1/3 \\ 0 & 5/3 & 2/3\end{matrix} \middle| \begin{matrix} -11 \\ 2/3 \\ 13/3 \end{matrix} \right) \Rightarrow \left( \begin{matrix} -3 & -1 & 2 \\ 0 & 5/3 & 2/3 \\ 0 & 0 & 1/5\end{matrix} \middle| \begin{matrix} -11 \\ 13/3 \\ -1/5 \end{matrix} \right) \Rightarrow \left( \begin{matrix} 1 & 1/3 & -2/3 \\ 0 & 1 & 2/5 \\ 0 & 0 & 1\end{matrix} \middle| \begin{matrix} 11/3 \\ 13/5 \\ -1 \end{matrix} \right)

然后从下往上依次回带,就能得到所有的 xx 值。

  • 将增广矩阵中所有 aj,i (1in,j>i)a_{j,i} \ (1 \leq i \leq n,j>i) 全部变成 00

    根据加减消元法,用第 ii 行的系数将所有 aj,i (j>i)a_{j,i} \ (j > i) 都变成 00。具体来说,将 ai,ia_{i,i} 变成 aj,ia_{j,i},第 ii 行按照比例同时变化。然后对于第 jj 行的每一个数 aj,ka_{j,k} 都减去 ai,ka_{i,k}。相当于每个 aj,ka_{j,k} 减去 ai,jai,i=aj,i×ai,k\frac{a_{i,j}}{a_{i,i}} = a_{j,i} \times a_{i,k}

  • ai,ia_{i,i} 变成 11

    这个很简单,第 ii 行的每个数都除以 ai,ia_{i,i} 即可。

  • 从下到上回带求出每个 xx

四、注意点

  • 假设正在处理第 ii 行,则首先要先找一个 r>ir>i 且绝对值更大的 ar,ia_{r,i},然后交换第 ii 行和第 rr 行。这样可以减少误差。

  • 若当前 ai,i=0a_{i,i}=0,则高斯消元无法继续,无法确定 xix_i,判断无解。

五、代码

P3389 【模板】高斯消元法

#include <bits/stdc++.h>
using namespace std;
typedef long long ll;

inline ll read() {
    ll sum = 0, ff = 1;
    char ch = getchar();
    while(ch < '0' || ch > '9') {
        if(ch == '-') {
            ff = -1;
        }
        ch = getchar();
    }
    while(ch >= '0' && ch <= '9') {
    	sum = sum * 10 + ch - 48;
        ch = getchar();
    }
    return sum * ff;
}

void write(ll x) {
	if(x < 0)
	    putchar('-'), x = -x;
	if(x > 9)
		write(x / 10);
	putchar(x % 10 + '0');
}

const double eps = 1e-8;
const int N = 1e2 + 7;
double tmp, a[N][N], ans[N];
int n, pos;

inline double Abs(double x) {
	return x < 0 ? -x : x;
}

inline int cmp(double x, double y) {
	if(Abs(x - y) < eps)	
		return 0;
	return Abs(x) > Abs(y) ? 1 : -1;
}

int main() {
//  	freopen(".in", "r", stdin);
//  	freopen(".out", "w", stdout);
	n = read();
	for(int i = 1; i <= n; i++)
		for(int j = 1; j <= n + 1; j++)
			scanf("%lf", &a[i][j]);
	for(int i = 1; i <= n; i++) {
		//找到 pos 使得 pos > i 且 a[pos][i] 的绝对值 > a[i][i] 的绝对值 
		pos = i;
		for(int j = i + 1; j <= n; j++)
			if(cmp(a[j][i], a[pos][i]) == 1)
				pos = j;
		if(cmp(a[pos][i], 0) == 0)//无解 
			return puts("No Solution"), 0;
		for(int j = 1; j <= n + 1; j++)
			swap(a[pos][j], a[i][j]);//交换
		//使 a[i][i] = 1 
		tmp = a[i][i];
		for(int j = i; j <= n + 1; j++)
			a[i][j] /= tmp;
		//使所有 f[j][i] (j > i) = 0
		for(int j = i + 1; j <= n; j++) {
			tmp = a[j][i];
			for(int k = i; k <= n + 1; k++)
				a[j][k] -= a[i][k] * tmp;
			}
		} 
	//从下到上回带	
	for(int i = n; i >= 1; i--) {
		ans[i] = a[i][n + 1];
		for(int j = n; j >= i + 1; j--)
			ans[i] -= ans[j] * a[i][j];
	}
	for(int i = 1; i <= n; i++)
		printf("%.2f\n", ans[i]);
    return 0;
}

六、例题

P2455 [SDOI2006]线性方程组

用到 高斯-约旦消元法,难点在于判断无解和无数解。

七、参考文献

C++高斯消元详解

算法竞赛入门经典 训练指导