寒假讲课：线性代数

寒假讲课：线性代数

蓬蒿人

推荐3Blue1Brown《线性代数的本质》和漫士沉思录《无痛线代》系列。

在高中数学的基础上，有时间的话推荐阅读《线性代数导论（Introduction to Linear Algebra）》《线性代数应该这样学（Linear Algebra Done Right）》。前者是科普（？）和入门性质，偏工科应用；后者更加抽象（本义），适合第二遍学的时候看。后者的作者在网上免费发布电子版，中文版也是完全免费的。

向量和矩阵

向量

先让我们约定几个符号：

• 有序对：集合论一般用尖括号\(\langle1,4\rangle\)包裹有序对，但为了与一般的坐标形式符合（而且这里不涉及无序对），我们还是采用圆括号\((1,4)\)。
• 元组：有序的\(n\in\mathbb{N}\)个数称作\(n\)元组（tuple），同样使用圆括号包裹。
• 笛卡尔积：设\(S,T\)是集合，定义\(S\times T=\{(x,y):x\in S\land y\in T\}\)。特别地，记\(S\times S\)为\(S^2\)。
• 设\(S,T\)是集合，定义\(T^S\)为所有从\(S\)映射到\(T\)的函数，即\(T^S\triangleq\{f:(f:S\to T)\}\)。

几何意义

凡学过高中物理的同学应当都理解“矢量”的概念，例如力、速度、加速度、动量都是矢量。矢量（vector）就是向量（vector）。向量在几何上通常视作一个“箭头”。

通常我们说的“向量”都是自由向量，即认为起点不同但大小、方向相同的向量是相等的，与“有向线段”相区分。

定义

既有大小（magnitude）/长度（length）又有方向（direction）的量称为向量。向量的长度也称为模（module），用符号\(\Vert\cdot\Vert\)表示。特别地，称长度为\(0\)的向量为零向量（zero vector），记作\(\vec{0}\)。称长度为\(1\)的向量为单位向量（unit vector）。

称两个向量是平行的当且仅当它们方向相同或相反，称两个向量是垂直的当且仅当它们的方向垂直，即夹角为\(\dfrac{1}{2}\pi\)。特别地，规定\(\vec{0}\)与任何向量平行。

基本运算：

• 向量可以依照平行四边形法则相加。
• 向量可以进行数乘，定义为缩放和反向。
• 两个向量的数量/点/内积（inner product）定义为模长的乘积再乘上夹角的余弦，即\(\vec{u}\cdot\vec{v}\triangleq\Vert\vec{u}\Vert\Vert\vec{v}\Vert\cos\langle\vec{u},\vec{v}\rangle\)。

特别地，定义两个空间向量\(\vec{u},\vec{v}\)的外/向量/叉积\(\vec{u}\times\vec{v}\)为这样一个向量，它的模长为\(\Vert\vec{u}\Vert\Vert\vec{v}\Vert\sin\langle\vec{u},\vec{v}\rangle\)，方向按照右手定则确定，即把右手做成虚握的“点赞”姿势👍，手掌沿\(\vec{u}\)方向、四根手指向\(\vec{v}\)方向弯曲，则竖起的大拇指即沿\(\vec{u}\times\vec{v}\)方向，且与\(\vec{u},\vec{v}\)均垂直。

注意叉积的别名“外积”在线性代数语境下具有歧义（列向量与其转置进行矩阵乘法的结果也称为外积），注意避免使用。

性质

运算律：

• 向量加法满足交换律、结合律。
• 数乘对标量加法和向量加法都满足分配律：\((\lambda+\mu)(\vec{u}+\vec{v})=\lambda\vec{u}+\lambda\vec{v}+\mu\vec{u}+\mu\vec{v}\)。
• 内积满足交换律和对向量加法的分配律：\(\vec{u}\cdot(\vec{v}+\vec{w})=\vec{u}\cdot\vec{v}+\vec{u}\cdot\vec{w}\)。

平面向量基本定理：若\(O,A,B\)不共线，则\(C\in\text{平面}ABO\iff\exists\lambda,\mu\in\mathbb{R}\ (\overrightarrow{OC}=\lambda\overrightarrow{OA}+\mu\overrightarrow{OB})\)。

空间向量基本定理：若\(O,A,B,C\)不共面，则三维空间内任意一点\(D\)都满足\(\exists\lambda,\mu,\nu\in\mathbb{R}\ (\overrightarrow{OD}=\lambda\overrightarrow{OA}+\mu\overrightarrow{OB}+\nu\overrightarrow{OC})\)。

共线向量基本定理：设\(\vec{v}\neq\vec{0}\)，则\(\vec{u}\parallel\vec{v}\iff\exists\lambda\in\mathbb{R}\ (\vec{u}=\lambda\vec{v})\)。

推论：设\(O,A,B,C\)是空间内三点，则\(C\in\text{直线}AB\iff\exists\lambda,\mu\in\mathbb{R}\ (\lambda+\mu=1\land\overrightarrow{OC}=\lambda\overrightarrow{OA}+\mu\overrightarrow{OB})\)。

推论：设\(O,A,B,C,D\)是空间内三点，\(O,A,B,C\)不共面，则\(D\in\text{平面}ABC\iff\exists\lambda,\mu,\nu\in\mathbb{R}\ (\lambda+\mu+\nu=1\land\overrightarrow{OD}=\lambda\overrightarrow{OA}+\mu\overrightarrow{OB}+\nu\overrightarrow{OC})\)。

向量垂直的判定：若\(\vec{u},\vec{v}\neq\vec{0}\)，则\(\vec{u}\perp\vec{v}\iff\vec{u}\cdot\vec{v}=0\)。

将若干个向量数乘并相加的结果称为它们的线性组合（linear combination）。对于一组向量，如果其中任何一个向量都不能由其他的向量线性组合得到，则称这组向量互相线性无关（linearly independent），否则称为线性相关（linearly dependent）。线性无关的一组向量称为一组基(底)（basis）。利用基底，可以将向量运算分解为各分量上的标量运算。

代数意义

垂直的基底是最讨人喜欢的，因为交叉项的内积为\(0\)，避免了烦人的夹角余弦。而如果基底都是单位向量就更爽了。

以平面向量为例：设\(\vec{i},\vec{j}\)是一组正交单位基，\(\vec{v}_1=x_1\vec{i}+y_1\vec{j},\ \vec{v}_2=x_2\vec{i}+y_2\vec{j}\)，则

\[\begin{aligned} \vec{v}_1\cdot\vec{v}_2&=(x_1\vec{i}+y_1\vec{j})\cdot(x_2\vec{i}+y_2\vec{j})\\ &=x_1x_2\vec{i}^2+(x_1y_2+x_2y_1)(\vec{i}\cdot\vec{j})+y_1y_2\vec{j}^2\\ &=x_1x_2+y_1y_2. \end{aligned} \]

在笛卡尔（Descartes）坐标系中，总是把向量的起点固定在原点，则两个向量不同当且仅当其终点坐标不同，于是我们可以用终点坐标来唯一描述向量，这就是向量的坐标表示。分别取\(+x,+y,+z\)方向的单位向量为\(\vec{i},\vec{j},\vec{k}\)，则终点在\((x,y,z)\)处的向量的基底表示刚好就是\(x\vec{i}+y\vec{j}+z\vec{k}\)。

于是向量加法、数乘和内积都有了代数意义，即如果\(\vec{v}_1=(x_1,y_1,z_1),\ \vec{v}_2=(x_2,y_2,z_2),\ \lambda\in\mathbb{R}\)，则

\[\vec{v}_1+\vec{v}_2=(x_1+x_2,y_1+y_2,z_1+z_2),\\ \lambda\vec{v_1}=(\lambda x_1,\lambda y_1,\lambda z_1),\\ \vec{v}_1\cdot\vec{v}_2=x_1x_2+y_1y_2+z_1z_2,\\ \vec{v}_1\times\vec{v}_2=(y_1z_2-y_2z_1,z_1x_2-z_2x_1,x_1y_2-x_2y_1). \]

在线性代数中，我们通常将坐标写成一列，例如

\[(1,1,4)=\begin{bmatrix} 1\\ 1\\ 4 \end{bmatrix}, \]

称为列向量。如果写成一行\(\begin{bmatrix}1 & 1 & 4\end{bmatrix}\)，则称为行向量。在不加额外说明的情况下，“向量”默认指列向量，但常写作\((1,1,4)\)，因为竖着写不方便，无论是手写还是\(\LaTeX\)。

在笛卡尔坐标系下：直线对应\(\mathbb{R}^1\)，在其中描述一个向量，只需用一个实数，此时向量与标量不需要区分；平面对应\(\mathbb{R}^2\)，在其中描述一个向量需要两个实数；立体对应\(\mathbb{R}^3\)，在其中描述向量需要三个实数。超过三维的空间是人类难以想象的，但这启发我们用\(n\in\mathbb{N}^*\)个实数来描述\(n\)维空间\(\mathbb{R}^n\)中的向量，即我们可以用一种纯代数的方式定义\(n\)维空间的向量，而无需理解高维空间的样貌。实际上，这种定义方式的向量是更让人感觉“严谨”的，因为它不要求我们对立体几何有任何理解（无需欧氏几何的五条公理），可以严格建立在仅以ZFC公理系统为基石的集合论上。

---

不严谨地说，\(\mathbb{R}^n\ (n\in\mathbb{N}^*)\)中的元素，即\(n\)元有序实数元组，称为\(n\)维向量。

于是可以重新定义各运算：设\(\vec{u}=(u_i)_{i=1}^n,\ \vec{v}=(v_i)_{i=1}^n,\ \lambda\in\mathbb{R}\)。

• 向量加法：\(\vec{u}+\vec{v}=(u_i+v_i)_{i=1}^n\)。
• 向量数乘：\(\lambda\vec{v}=(\lambda v_i)_{i=1}^n\)。
• 内积：\(\vec{u}\cdot\vec{v}=\sum_{i=1}^nu_iv_i\)。特别地，将\(\vec{v}\cdot\vec{v}\)记作\(\vec{v}^2\)。
• 模：\(\Vert\vec{v}\Vert=\sqrt{\vec{v}^2}=\sqrt{\sum_{i=1}^nv_i^2}\)。
• 夹角：\(\langle\vec{u},\vec{v}\rangle=\arccos\dfrac{\vec{u}\cdot\vec{v}}{\Vert\vec{u}\Vert\Vert\vec{v}\Vert}\)。

由于复数允许实数的所有运算，所以我们完全可以将上面描述向量的方式推广到\(\mathbb{C}^n\)。

至少对于工科学生而言，对向量的理解到这一步完全够用了。

最广义的向量

经过进一步抽象后，向量可以是满足一定条件的任何东西，不一定是箭头，也不一定是一列数。

数学家发现，无论是几何箭头、数组，还是多项式、连续函数、甚至是矩阵，它们都共享同一套运算结构（加法和数乘）。为了统一研究这些对象，数学家提出了向量空间的概念。

域

一个集合\(\mathbb{F}\)如果包含至少两个不同的元素（分别记作\(0,1\)），且对二元运算\(+,\times\)（分别称作加法、乘法）满足

• 交换律：\(\forall x,y\in\mathbb{F}\ (x+y=y+x\land x\times y=y\times x)\)。
• 结合律：\(\forall x,y,z\in\mathbb{F}\ ((x+y)+z=x+(y+z)\land(x\times y)\times z=x\times(y\times z))\)。
• 存在恒等元：\(\forall x\in\mathbb{F}\ (x+0=x\land x\times1=x)\)，即\(0,1\)分别是\(+,\times\)的单位/幺/恒等元（identity）。
• 存在逆元：对任何\(x\in\mathbb{F}\)存在唯一的\(y\in\mathbb{F}\)使\(x+y=0\)；对任何\(x\in\mathbb{F}\)满足\(x\neq0\)，存在唯一的\(y\in\mathbb{F}\)使\(x\times y=1\)。
• 乘法对加法满足分配律：\(\forall x,y,z\in\mathbb{F}\ (x\times(y+z)=x\times y+x\times z)\)。

则称\(\mathbb{F}\)是一个域（field）。

\(\mathbb{Q},\mathbb{R},\mathbb{C}\)都对普通加法和普通乘法构成域。\(\mathbb{Z}_m\)对模\(m\)加法和模\(m\)乘法构成域；特别地，\(\{0,1\}\)对异或\(\oplus\)和与\(\land\)构成域。

向量空间

显然对任何域\(\mathbb{F}\)，\(\mathbb{F}^n\ (n\in\mathbb{N}^*)\)一定是向量空间，只需将加法定义为逐元素的\(\mathbb{F}\)上的加法。更一般地，对任何域\(\mathbb{F}\)和集合\(S\)，\(\mathbb{F}^S\)构成向量空间。\(\mathbb{F}^n\)可以视作\(\mathbb{F}^{\mathbb{N}\cap[1,n]}\)。

向量空间的元素称为向量（vector）或点（point）。

可以证明一个向量空间\(V\)必然具有唯一的加法恒等元\(0\)，\(\forall v\in V\ (0v=0)\)（第一个\(0\)是数，第二个\(0\)是向量）。每个元素具有唯一的加法逆元，记作\(-v\)。易证\(-v=(-1)v\)。

如果向量空间\(V\)的子集\(U\)是与\(V\)具有相同的加法恒等元、加法和标量乘法运算的向量空间，则称\(U\)是\(V\)的子空间（subspace）。另一个等价的定义是，若向量空间\(V\)的子集\(U\)满足\(0\in U\)且对加法和标量乘法封闭，则称\(U\)是\(V\)的子空间。

线性组合和张成空间

设\(v_1,\dots,v_m\)是\(V\)中的一组向量，称形如\(\sum_{i=1}^ma_iv_i\ (a_1,\dots,a_m\in\mathbb{F})\)的向量是\(v_1,\dots,v_m\)的线性组合（linear combination）。

向量空间\(V\)中的一个向量组\(v_1,\dots,v_m\)的所有线性组合构成的集合称为这组向量的张成空间（span），记作\(\operatorname{span}(v_1,\dots,v_m)\)，即

\[\operatorname{span}(v_1,\dots,v_m)\triangleq\left\{\sum_{i=1}^ma_iv_i:\forall i\in\mathbb{N}\cap[1,m]\ (a_i\in\mathbb{F})\right\}. \]

可以证明这个张成空间一定是所有包含该向量组中所有向量的\(V\)的子空间中最小的。

如果一个向量空间可以由有限个向量张成，则称为有限维向量空间，否则称为无限维向量空间。\(\mathbb{F}^n\)是有限维向量空间，而系数都在\(\mathbb{F}\)中的所有多项式构成的集合是无限维向量空间。

对于向量空间\(V\)中的一个向量组\(v_1,\dots,v_m\)，若\(\sum_{i=1}^ma_iv_i=0\ (a_1,\dots,a_m\in\mathbb{F})\)当且仅当\(\forall i\in\mathbb{N}\cap[1,m]\ (a_i=0)\)，则称这组向量是线性无关（linearly independent）的，否则称为线性相关（linearly dependent）的。

线性相关性引理：如果一组向量线性相关，则其中存在某一个向量，使将其移除后，这组向量张成的空间不变。

线性空间\(V\)中一个线性无关且张成\(V\)的向量组称为\(V\)的一个基（basis）。另一个等价的定义是，如果向量空间\(V\)中的一组向量满足每个\(v\in V\)都可以由这组向量唯一地线性表示，则这组向量是\(V\)的一个基。

向量空间的每个张成组都能被削减成该向量空间的一个基；有限维向量空间都有基，且其每个线性无关向量组都可被扩充成该向量空间的一个基。

一个有限维向量空间\(V\)的基都由相同个数的向量组成，将这个数量称作\(V\)的维数（dimension），记作\(\dim V\)。

若\(U\)是\(V\)的子空间，则\(\dim U\leq\dim V\)，当且仅当\(U=V\)时取等。

\(\dim V\)个线性无关的向量构成\(V\)的一个基。

内积和范数

如果向量空间\(V\)上的运算\(\cdot:V^2\to\mathbb{F}\)满足

• 正性：\(\forall v\in V\ (v\cdot v\geq0)\)。
• 定性：\(\forall v\in V\ (v\cdot v=0\iff v=0)\)（前式中的\(0\)是数，而后式中的则是向量）。
• 第一个运算数的可加性：\(\forall u,v,w\in V\ ((u+v)\cdot w=u\cdot w+v\cdot w)\)。
• 第一个运算数的齐次性：\((\forall u,v\in V)\ (\forall\lambda\in\mathbb{F})\ ((\lambda u)\cdot v=\lambda(u\cdot v))\)。
• 共轭对称性：\(\forall u,v\in V\ (u\cdot v=\overline{v\cdot u})\)。

则称\(\cdot\)是\(V\)上的内积（inner product）。对\(v\in V\)，定义其范数（norm）为\(\Vert v\Vert=\sqrt{v\cdot v}\)。

范数的基本性质：

• \(\Vert v\Vert=0\iff v=0\)（第一个\(0\)是数，第二个\(0\)是向量）。
• \(\Vert\lambda v\Vert=|\lambda|\Vert v\Vert\)。

如果两个向量的内积是\(0\)，则称这两个向量是正交（orthogonal）的。

毕达哥拉斯定理（Pythagorean theorem）：若\(u\perp v\)，则\(\Vert u+v\Vert^2=\Vert u\Vert^2+\Vert v\Vert^2\)。

从线性方程组到矩阵

一次方程也称为线性方程（linear equation），由线性方程组成的方程组称为线性方程组。

线性方程组的行视角和列视角

以方程组

\[\begin{cases} x-2y=1\\ 3x+2y=11 \end{cases} \]

为例。

• 行视角：两条直线求交点。

• 列视角：怎么样用向量\((1,3)\)和\((-2,2)\)线性组合成\((1,11)\)？

  \[x \begin{bmatrix} 1\\ 3 \end{bmatrix} +y \begin{bmatrix} -2\\ 2 \end{bmatrix} = \begin{bmatrix} 1\\ 11 \end{bmatrix}. \]

矩阵-向量乘法的行视角和列视角

将上面的方程组写作矩阵乘向量的形式：

\[\begin{bmatrix} 1 & -2\\ 3 & 2 \end{bmatrix} \begin{bmatrix} x\\ y \end{bmatrix} = \begin{bmatrix} 1\\ 11 \end{bmatrix}. \]

• 行视角：

  \[\begin{bmatrix} (1,-2)\cdot(x,y)\\ (3,2)\cdot(x,y) \end{bmatrix} = \begin{bmatrix} 1\\ 11 \end{bmatrix}. \]

• 列视角：

  \[x \begin{bmatrix} 1\\ 3 \end{bmatrix} +y \begin{bmatrix} -2\\ 2 \end{bmatrix} = \begin{bmatrix} 1\\ 11 \end{bmatrix}. \]

由此我们可以得到矩阵-向量乘法的两种等价定义：

• 行视角：矩阵各行与被乘向量的内积。
• 列视角：矩阵各列按照被乘向量线性组合。

矩阵乘法

并非任意两个矩阵都能相乘。两个矩阵能相乘，当且仅当第一个矩阵的列数等于第二个矩阵的行数。

• 列视角：前矩阵依次乘后矩阵的各列。

• 行视角：前矩阵的各行依次乘后矩阵。

• 设\(A\in\mathbb{F}^{l\times m},\ B\in\mathbb{F}^{m\times n},\ C=AB\)，则\(C\in\mathbb{R}^{l\times n}\)，且

  \[C(i,j)=\sum_{k=1}^mA(i,k)B(k,j)\ (i,j\in\mathbb{N}^*,\ i\leq l,\ j\leq n). \]

容易注意到矩阵乘法对两个矩阵的顺序有严格要求。可以证明，矩阵乘法不满足交换律，但满足结合律。

计算矩阵乘法的朴素方法需要\(n^3\)次乘法。

The best at this moment is \(n^{2.376}\). But the algorithm is so awkward that scientific computing is done the regular way: \(n^2\) dot products in \(AB\), and \(n\) multiplications for each one.

矩阵的初等变换

单位矩阵

\[I= \begin{bmatrix} 1\\ & \ddots\\ & & 1 \end{bmatrix} \]

是矩阵乘法的单位元，它与任何矩阵相乘（无论左乘还是右乘）的结果都等于被乘矩阵。

如果我们将单位矩阵的第\(i\)行的\(1\)改为\(\lambda\)得到矩阵\(M\)，则\(MA\)的结果是\(A\)的第\(i\)行统统变为原来的\(\lambda\)倍，而其余行不变；\(AM\)的结果是\(A\)的第\(i\)列统统变为原来的\(\lambda\)倍，而其余行不变。

如果我们交换单位矩阵的两行，例如

\[P= \begin{bmatrix} 1\\ & & 1\\ & 1 \end{bmatrix}, \]

则\(PA\)的结果是交换\(A\)的第\(2,3\)行，其余行不变；\(AP\)的结果是交换\(A\)的第\(2,3\)列，其余列不变。进一步地，如果我们任意排列单位矩阵的行得到\(P\)，则\(PA\)是按照\(P\)的行来重新排列\(A\)的各行，\(AP\)是按照\(P\)的列重新排列\(A\)的各列。

如果将矩阵\(A\)左乘

\[E= \begin{bmatrix} 1\\ & 1\\ -1 & & 1 \end{bmatrix}, \]

所得结果是\(A\)的第三行减去第一行，其余行不变；而\(AE\)则是将\(A\)的第\(1\)列减去第\(3\)列，其余列不变。

交换两行、将某一行变为原来的非零倍、将某一行减去另一行统称为矩阵的初等行变换，类似地可以定义初等列变换，初等行变换和初等列变换统称为矩阵的初等变换。这些变换都可以用矩阵来描述，称为初等矩阵。由于矩阵乘法的可结合性（也可以说初等行变换的可结合性这正是矩阵乘法可结合的原因），\(EPA=E(PA)=(EP)A\)，也就是说，若干个初等行变换可以结合到一起，再一次性对被乘矩阵起作用。初等矩阵的乘积仍为初等矩阵。

如果一个矩阵\(A\)可以经过若干次初等行变换（即左乘若干个初等矩阵；当然，由于结合性，等价于左乘一个初等矩阵）得到矩阵\(B\)，则称\(A\)与\(B\)是行等价的。类似地可以定义列等价的概念。

矩阵的其他基本运算

相同形状的矩阵可以相加，定义为对应位置元素相加。矩阵乘法对矩阵加法满足分配律，即

\[(A+B)C=AC+BC,\ A(B+C)=AB+AC. \]

矩阵同样可以数乘，定义为各位置元素分别乘上乘数。数乘对矩阵加法、矩阵乘法对普通加法都满足分配律。

把整个矩阵沿主对角线对称一下，这个操作称为矩阵转置（transpose）。更严谨地，设\(A\in\mathbb{F}^{m\times n}\)，则它的转置\(A^\top\in\mathbb{F}^{n\times m}\)，且

\[A^\top(j,i)=A(i,j)\ (i,j\in\mathbb{N},\ i\leq m,\ j\leq n). \]

称\(A^\top B\)为内积，称\(AB^\top\)为外积。内积还是内积，但外积却跟叉乘完全不同。

旋转变换矩阵

矩阵-向量乘法的本质是线性变换。将向量左乘一个矩阵，本质上是按照这个矩阵所描述的规则对这个向量进行旋转和拉伸。

将一个二维向量逆时针旋转角\(\delta\)的矩阵为

\[R(\delta)= \begin{bmatrix} \cos\delta & -\sin\delta\\ \sin\delta & \cos\delta \end{bmatrix}. \]

解线性方程组的高斯消元法

“消元法”实际上有更一般的描述，高斯消元法（Gaussian elimination）、克劳特分解法（Crout decomposition method）只是其中特殊的两种，它们分别对普适的消元法中的参数选取方式做出了规定。从数值分析的角度来说，消元法没有方法误差，但舍入误差不可避免，所以需要根据数据的具体情况选择合适的参数。在算法竞赛中一般不需要考虑这个问题。

设我们要解的线性方程组是

\[\begin{cases} a_{1,1}x_1+a_{1,2}x_2+\dots+a_{1,n}x_n=b_1\\ a_{2,1}x_1+a_{2,2}x_2+\dots+a_{2,n}x_n=b_2\\ \vdots\\ a_{m,1}x_1+a_{n,2}x_2+\dots+a_{m,n}x_n=b_m \end{cases}, \]

写成矩阵的形式是

\[A\vec{x}=\vec{b}, \]

其中

\[A= \begin{bmatrix} a_{1,1} & a_{1,2} & \dots & a_{1,n}\\ a_{2,1} & a_{2,2} & \dots & a_{2,n}\\ \vdots & \vdots & \ddots & \vdots\\ a_{m,1} & a_{m,2} & \dots & a_{m,n} \end{bmatrix} \in\mathbb{F}^{m\times n},\ \vec{x}= \begin{bmatrix} x_1\\ x_2\\ \vdots\\ x_n \end{bmatrix},\ \vec{b}= \begin{bmatrix} b_1\\ b_2\\ \vdots\\ b_n \end{bmatrix}. \]

消元

先让我们考虑最理想的情况，\(m=n\)，即\(A\)是一个方阵，而且我们假设不会发生消不了的情况。

在这种理想情况下，用第\(1\)行消掉第\(2\sim n\)行\(x_1\)的系数，即从第\(i\in\mathbb{N}\cap[2,n]\)行中减去\(\dfrac{a_{i,1}}{a_{1,1}}\)倍第\(1\)行；用第\(2\)行消掉第\(3\sim m\)行\(x_2\)的系数，即从第\(i\in\mathbb{N}\cap[3,n]\)行中减去\(\dfrac{a_{i,2}}{a_{2,2}}\)倍第\(2\)行；……这样直到第\(n\)行，会得到一个上三角形的方程组：

\[\left\{ \begin{array}{rrrrr} u_{1,1}x_1 & +u_{1,2}x_2 & +\dots & +u_{1,n}x_n & =c_1\\ & u_{2,2}x_2 & +\dots & +u_{2,n}x_n & =c_2\\ & & \ddots & \vdots\\ & & & u_{n,n}x_n &=c_n \end{array} \right., \]

这一步称为消元（elimination）。用来消元的元素称为主元（pivot）。消元后，留在主对角线上的就是主元。

回代

由最后一行可以解出\(x_n=\dfrac{c_n}{u_{n,n}}\)，然后反代回倒数第二行得\(x_{n-1}=\dfrac{c_{n-1}-u_{n-1,n}x_n}{u_{n-1,n-1}}\)，……直到把所有未知数全部求出。这一步称为回代（back substitution）。

以上是最最理想的情况，但现实情况没有这么好。如果在消元到某一行时，我们发现主元等于\(0\)了怎么办？例如

\[\left\{ \begin{aligned} 1x+1y+1z & =6\\ {\color{red}0y}+2z & =4\\ 3y-2z & =5 \end{aligned} \right. \]

如果直接按照上面的步骤走，会导致除数为\(0\)，然后爆炸。但实际上，这组方程是有解的（\(x=1,y=3,z=2\)）。显然我们可以在此时将后面的行中这一列系数不为\(0\)的行换过来，反正交换两行不会影响方程组的解。进一步地，如果是以浮点数的形式求解，则选择后面的行中该列最大的系数作为主元是使误差最小的。这就是选主元的高斯消元法。

到目前为止，我们的消元法可以解决所有有唯一解的\(n\)个\(n\)元线性方程了。甚么时候会无解？当方程之间有矛盾时，方程组无解。但“矛盾”是不能直接知道的；如果我们解到最后一个方程发现，左边是\(0x\)，右边却不为\(0\)，方程组就会无解。那甚么时候不止有一组解呢？当我们发现，消元到某一行时，从当前行往后当前列都是\(0\)——我们选不出主元，则这个未知数可以任意取值，此时有无穷多组解。

在消元过程中，\(\vec{b}\)需要和\(A\)一起参与初等行变换。为了方便，通常将它们合并到一起，称为增广矩阵

\[\tilde{A}= \begin{bmatrix} A & \vec{b} \end{bmatrix} = \left[ \begin{array}{cccc:c} a_{1,1} & a_{1,2} & \dots & a_{1,n} & b_1\\ a_{2,1} & a_{2,2} & \dots & a_{2,n} & b_2\\ \vdots & \vdots & \ddots & \vdots & \vdots\\ a_{m,1} & a_{m,2} & \dots & a_{m,n} & b_m \end{array} \right]. \]

消元的过程就是分解

\[\tilde{A}=L \begin{bmatrix} U & \vec{c} \end{bmatrix} \]

的过程，其中\(L\)是所有初等行变换的结合。由于\(L\)只与\(A\)有关而与\(\vec{b}\)无关，所以我们甚至可以同时求解多个系数矩阵相同的方程组，只需将方程右边统统放进增广矩阵即可。根据我们的消元步骤，\(L\)必定是一个下三角矩阵。因此消元的过程本质上也是将\(A\)分解为一个下三角矩阵\(L\)和一个上三角矩阵\(U\)之积，称为矩阵的LU分解。

那如果未知数的数量和方程的数量不一样呢？这时系数矩阵不是方阵，那我们注定没办法将方程组化为完全的上三角形，我们只能尽可能化简。

满足下面两个条件的矩阵称为行阶梯形（row echelon form）矩阵：

• 非零行最左边的首个非零元素，严格地比上一行的首个非零元素更靠右。
• 全零行都在矩阵的底部。

各行的首个非零元素都是\(1\)而且都是所在列唯一非零元素的行阶梯形矩阵称作简化行阶梯形（reduced row echelon form）矩阵。

通过将增广矩阵化为简化行阶梯形，我们可以求出方程组的通解或证明其无解。消元的过程基本相同，唯在找不到主元时将我们考虑的列右移一列；回代时，由于我们的目的不只是求出各个未知数的值，而且要知道未知数之间的关系，所以需要用下面的行消掉主元所在列的其他元素。

模板题：P2455 [SDOI2006] 线性方程组

注意此题要求区分无解和无穷多组解。

inline void solve() {
	constexpr double EPS = 1e-9;
	size_t n;
	std::cin >> n;
	std::vector<std::vector<double>> mat(n, std::vector<double>(n + 1));
	for (auto &row : mat) {
		for (auto &e : row) {
			std::cin >> e;
		}
	}

	// Elimination
	for (size_t i = 0, col = 0; i < n; ++i, ++col) {
		while (col < n) {
			// Pivoting
			size_t pvt_row = i;
			for (size_t j = i + 1; j < n; ++j) {
				if (std::abs(mat[j][col]) > std::abs(mat[pvt_row][col])) {
					pvt_row = j;
				}
			}
			std::swap(mat[i], mat[pvt_row]);
			if (std::abs(mat[i][col]) >= EPS) {
				break;
			}
			++col;
		}

		if (col == n) {
			int ans = 0;
			for (size_t j = i; j < n; ++j) {
				if (std::abs(mat[j][n]) >= EPS) {
					ans = -1; // no solution
					break;
				}
			}
			std::cout << ans << '\n';
			return;
		}

		for (size_t j = i + 1; j < n; ++j) {
			double fac = mat[j][col] / mat[i][col];
			for (size_t k = col; k < n + 1; ++k) {
				mat[j][k] -= fac * mat[i][k];
			}
		}
	}

	// Back substitution
	std::vector<double> res(n);
	for (size_t i = n; i-- > 0;) {
		res[i] = mat[i][n];
		for (size_t j = i + 1; j < n; ++j) {
			res[i] -= mat[i][j] * res[j];
		}
		res[i] /= mat[i][i];
	}

	std::cout << std::setprecision(10) << std::fixed;
	for (size_t i = 0; i < n; ++i) {
		std::cout << 'x' << (i + 1) << '=' << res[i] << '\n';
	}
}

时间复杂度：\(\varTheta(n^3)\)。
空间复杂度：\(\varTheta(n)\)额外空间保存解。

对于大型/稀疏矩阵，可以采用迭代算法求近似解。算法竞赛中不会涉及。

模板题：P2447 [SDOI2010] 外星千足虫

这是一道GF(2)域上的消元。

inline void solve() {
	constexpr size_t N = 1000;
	constexpr char EVEN[] = "Earth", ODD[] = "?y7M#", INF_SLNS[] = "Cannot Determine";

	size_t m, n;
	std::cin >> n >> m;
	std::vector<std::bitset<N + 1>> mat(m);
	for (auto &row : mat) {
		char ch;
		for (size_t i = 0; i <= n; ++i) {
			std::cin >> ch;
			row[i] = (ch == '1');
		}
	}

	size_t r = n;
	for (size_t i = 0; i < n; ++i) {
		size_t pvt_row = m;
		for (size_t j = i; j < m; ++j) {
			if (mat[j][i]) {
				pvt_row = j;
				break;
			}
		}
		if (pvt_row == m) {
			std::cout << INF_SLNS << '\n';
			return;
		}
		if (pvt_row >= r) {
			r = pvt_row + 1;
		}
		std::swap(mat[i], mat[pvt_row]);
		for (size_t j = i + 1; j < m; ++j) {
			if (mat[j][i]) {
				mat[j] ^= mat[i];
			}
		}
	}

	std::bitset<N + 1> res{};
	for (size_t i = n; i--;) {
		res[i] = mat[i][n] ^ ((mat[i] & res).count() & 1);
	}

	std::cout << r << '\n';
	for (size_t i = 0; i < n; ++i) {
		std::cout << (res[i] ? ODD : EVEN) << '\n';
	}
}

时间复杂度：\(\varTheta\left(\dfrac{m\cdot n^2}{W}\right)\)，其中\(W\)是计算机的字长。

逆矩阵

设\(A\)为矩阵，若存在矩阵\(A^{-1}\)满足\(A^{-1}A=I\)，则称\(A^{-1}\)为\(A\)的左逆矩阵；若\(AA^{-1}=I\)，则称\(A^{-1}\)为\(A\)的右逆矩阵。

若\(A\)的左右逆矩阵都存在，则它们必相等，称为\(A\)的逆矩阵，此时称\(A\)是可逆的。一个矩阵是可逆的，当且仅当它是满秩的方阵。

容易证明若\(A,B\)都可逆，则\((AB)^{-1}=B^{-1}A^{-1}\)。

利用消元法可以求逆矩阵。首先取增广矩阵\(\tilde{A}=\begin{bmatrix}A & I\end{bmatrix}\)，利用消元法将\(\tilde{A}\)化为简化行阶梯形。\(A\)是可逆的当且仅当\(\operatorname{rref}\tilde{A}\)的左半部分是\(I\)，此时右半部分必为\(A^{-1}\)，因为

\[A^{-1}\tilde{A}= \begin{bmatrix} A^{-1}A & A^{-1}I \end{bmatrix} = \begin{bmatrix} I & A^{-1} \end{bmatrix}. \]

取\(\operatorname{rref}\tilde{A}\)的右半部分即可。

模板题：P4783 【模板】矩阵求逆

inline void solve() {
	constexpr uint32_t MOD = 1e9 + 7;

	size_t n;
	std::cin >> n;
	std::vector<std::vector<uint32_t>> mat(n, std::vector<uint32_t>(n << 1));
	for (size_t i = 0; i < n; ++i) {
		mat[i][n + i] = 1;
		for (size_t j = 0; j < n; ++j) {
			std::cin >> mat[i][j];
		}
	}

	for (size_t i = 0; i < n; ++i) {
		size_t pvt_row = -1;
		for (size_t j = i; j < n; ++j) {
			if (mat[j][i]) {
				pvt_row = j;
				break;
			}
		}
		if (pvt_row == size_t(-1)) {
			std::cout << "No Solution\n";
			return;
		}
		std::swap(mat[i], mat[pvt_row]);

		uint32_t inv = modMulInv(mat[i][i], MOD);
		for (size_t j = i; j < (n << 1); ++j) {
			mat[i][j] = SC<uint64_t>(mat[i][j]) * inv % MOD;
		}

		for (size_t j = i + 1; j < n; ++j) {
			uint64_t fac = mat[j][i];
			mat[j][i] = 0;
			for (size_t k = i + 1; k < (n << 1); ++k) {
				mat[j][k] = (MOD - fac * mat[i][k] % MOD + mat[j][k]) % MOD;
			}
		}
	}

	for (size_t i = n; i--;) {
		for (size_t j = i; j--;) {
			uint64_t fac = mat[j][i];
			mat[j][i] = 0;
			for (size_t k = i + 1; k < (n << 1); ++k) {
				mat[j][k] = (MOD - fac * mat[i][k] % MOD + mat[j][k]) % MOD;
			}
		}
	}

	for (size_t i = 0; i < n; ++i) {
		for (size_t j = n; j < (n << 1); ++j) {
			std::cout << mat[i][j] << ' ';
		}
		std::cout << '\n';
	}
}

时间复杂度：\(\varTheta(n^3)\)。
空间复杂度：\(\varTheta(n^3)\)。

行列式

行列式的递归定义：

• 一个\(n\)阶方阵\(A\)的行列式（determinant）定义为它的某一行每个元素依次乘以对应的代数余子式之和，记作\(\det A\)。
• 矩阵中一个元素\(a_{i,j}\)的余子阵（minor）\(M_{i,j}\)定义为原矩阵去掉第\(i\)行和第\(j\)列，而其余元素相对位置不变构成的矩阵。
• \(a_{i,j}\)的代数余子式（cofactor）定义为\((-1)^{i+j}\det M_{i,j}\)。

直接按照定义计算行列式，复杂度为\(\varTheta(n!)\)。

容易证明上三角矩阵、下三角矩阵的行列式都等于其主对角线上元素之积，而消元法恰好可以帮我们将矩阵化为上三角形。只不过，消元的过程会稍稍改变其行列式，需要我们人为地维护：

• 交换两行会导致行列式变为相反数。
• 将一行的\(\lambda\in\mathbb{F}\)倍加到另一行上不改变行列式的值。
• 将一行变为原来的\(\lambda\in\mathbb{F}\)倍会导致行列式也变为原来的\(\lambda\)倍。

没那么模板的模板题：P7112 【模板】行列式求值

Wrong Answer

这道题不但可能以合数为模数，而且可能出现逆元不存在的情况。如果像我一样直接用exGcd求逆元会导致将不是\(0\)的误判为\(0\)，然后获得WA掉7个点。

inline void solve() {
	size_t n;
	uint32_t mod;
	std::cin >> n >> mod;
	std::vector<std::vector<uint32_t>> mat(n, std::vector<uint32_t>(n));
	for (auto &row : mat) {
		for (auto &e : row) {
			std::cin >> e;
		}
	}

	uint32_t det = 1;
	for (size_t i = 0; i < n; ++i) {
		// Pivoting
		size_t pvt_row = n;
		for (size_t j = i; j < n; ++j) {
			if (mat[j][i]) {
				pvt_row = j;
				break;
			}
		}
		if (pvt_row == n) {
			std::cout << "0\n";
			return;
		}
		if (pvt_row != i) {
			std::swap(mat[i], mat[pvt_row]);
			det = (mod - det) % mod; // det *= -1
		}

		uint32_t inv = modMulInv(mat[i][i], mod);
		for (size_t j = i + 1; j < n; ++j) {
			uint64_t fac = static_cast<uint64_t>(mat[j][i]) * inv % mod;
			for (size_t k = i; k < n; ++k) {
				mat[j][k] = (mod - fac * mat[i][k] % mod + mat[j][k]) % mod;
			}
		}
	}

	for (size_t i = 0; i < n; ++i) {
		det = SC<uint64_t>(det) * mat[i][i] % mod;
	}

	std::cout << det << '\n';
}

正解

正确的做法是采用fraction-free (Euclidean) Gaussian elimination，其优点在于避免除法（只用了整除而且保证安全）。我们来看看它是怎么实现的。

假设目前我们在第\(i\)行（设为\(r_i\)），想要用第\(i\)行消掉第\(j\ (j>i)\)行（设为\(r_j\)）的第\(i\)列。本来应该是

\[r_j\gets r_j-\frac{a_{j,i}}{a_{i,i}}r_i, \]

但这里可能\(a_{j,i}^{-1}\)不存在。正确的做法是：

1. 计算\(q\gets\left\lfloor\dfrac{a_{i,i}}{a_{j,i}}\right\rfloor\)。
2. 令\(r_j\gets r_j-q\cdot r_i\)。
3. 交换\(r_i\)与\(r_j\)。
4. 重复以上步骤直到\(a_{j,i}=0\)。

这个算法的正确性来源于辗转相除法的正确性。在辗转相除法中，出于\(\gcd(a,b)=\gcd(b,a\bmod b)\)，我们每次令\((a,b)\gets(b,a\bmod b)\)直至\(b=0\)，而\(a\bmod b=a-\left\lfloor\dfrac{a}{b}\right\rfloor\cdot b\)。

如果\(a_{i,i}=0\)而存在\(a_{j,i}\neq0\)的话，一定会在遇到\(j\)时将两行直接交换，从而将主元换上来。如果遍历完所有的\(j\)仍然\(a_{i,i}=0\)，说明从第\(i\)行往下这一列都是\(0\)，就可以宣告行列式为\(0\)了。

进行“辗转相除”时注意维护行列式的值。

	uint32_t det = 1;
	for (size_t i = 0; i < n; ++i) {
		for (size_t j = i + 1; j < n; ++j) {
			while (mat[j][i]) {
				uint64_t q = mat[i][i] / mat[j][i];
				for (size_t k = i; k < n; ++k) {
					mat[i][k] = (mod - q * mat[j][k] % mod + mat[i][k]) % mod;
				}
				std::swap(mat[i], mat[j]);
				det = (mod - det) % mod; // det = -det
			}
		}
		if (!mat[i][i]) {
			det = 0;
			break;
		}
		det = SC<uint64_t>(det) * mat[i][i] % mod;
	}

	std::cout << det << '\n';

时间复杂度：\(O(n^3+n^2\log^2(mod))\)，因为在关于某一列执行辗转相除消元时复杂度均摊到\(O(n)\)行上，对于固定的i，while (mat[j][i])循环的总次数上界为\(O(n+\log^2(mod))\)。

注意此题的时间限制是\(5\ \mathrm{s}\)。

线性基

如果给你了一组向量，如何从中提取出它们的基/极大线性无关组？

普适方法

给定\(\vec{v}_1,\dots,\vec{v}_m\in\mathbb{F}^n\ (m,n\in\mathbb{N}^*)\)，求\(\operatorname{span}(\vec{v}_1,\dots,\vec{v}_m)\)的一个基。

1. 构造矩阵

   \[M= \begin{bmatrix} \vec{v}_1^\top\\ \vdots\\ \vec{v}_m^\top \end{bmatrix}. \]

2. 执行消元法，将\(M\)化为行阶梯形矩阵（不用到简化行阶梯形矩阵）。

3. 提取基底：矩阵中剩下的非全零行就是原向量组张成空间的基。

这样得到的基可能不全是原本的向量。如果想知道怎么完全从原来的向量中抽出基（而非它们的线性组合），只需在消元时额外维护每一行最初的位置即可。

要检验一个向量能否由一个基线性表示，只需用基中的向量依次对该向量消元即可。

GF(2)域

给定\(v_1,\dots,v_m\in \mathrm{GF}(2)^n\ (m,n\in\mathbb{N}^*)\)，求其基。

GF(2)域仅有两个元素，这使得对于一组基向量，每个基向量会贡献至少一个二进制位（别的基向量都没拥有这一位）。于是可以开辟一个长为\(n\)的数组\(b\)，用于保存基，第\(i\)项表示贡献了第\(i\)位的基。初始时都设为\(0\)，表示没有基贡献这一位。依次遍历各个向量，检查其非\(0\)位。如果\(v_{i,j}=1\)：

1. 若\(b_j=0\)：目前的基中还没有能贡献第\(j\)位的，可以将\(v_i\)加入基中，即令\(b_j\gets v_i\)。退出循环（即使它有其它贡献的位也不能重复记录）。
2. 若\(b_j\neq0\)：目前的基中已有向量\(b_j\)贡献了第\(j\)位，跳过。

给定\(v_1,\dots,v_m\in \mathrm{GF}(2)^n\ (m,n\in\mathbb{N}^*)\)，在其中选取任意个，求其最大异或和。

在求线性基的基础上，我们需要从高到低遍历非\(0\)位，且在发生第2.种情况时令\(v_i\gets v_i\oplus b_j\)。这样得到的\(b\)，满足\(b_i\)的最高位是第\(i\)位。要求所有子集的最大异或和，只需从高位到低位考虑，如果当前位的基向量对答案有贡献就合并到答案中，而低位的部分要么也由该向量贡献，要么将来可以随意线性组合，所以无需担心。

inline void solve() {
	constexpr size_t B = 50;
	std::array<uint64_t, B> bases{};
	size_t n;
	std::cin >> n;
	for (size_t i = 0; i < n; ++i) {
		uint64_t x;
		std::cin >> x;
		for (size_t j = B; j--;) {
			if ((x >> j) & 1) {
				if (!bases[j]) {
					bases[j] = x;
					break;
				}
				x ^= bases[j];
			}
		}
	}

	uint64_t res = 0;
	for (size_t i = B; i--;) {
		res = std::max(res, res ^ bases[i]);
	}
	std::cout << res << '\n';
}

时间复杂度：\(O\left(mn\max\left\{\dfrac{n}{W},1\right\}\right)\)，其中\(W\)为位长。

广义矩阵乘法

矩乘快速幂加速递推

快速求出Fibonacci数列的第\(n\)项。

将递推式写作矩阵乘法形式：

\[\begin{bmatrix} f_{n+1}\\ f_{n+2} \end{bmatrix} = \begin{bmatrix} 0 & 1\\ 1 & 1 \end{bmatrix} \begin{bmatrix} f_n\\ f_{n+1} \end{bmatrix}, \]

于是

\[\begin{bmatrix} f_n\\ f_{n+1} \end{bmatrix} = \begin{bmatrix} 0 & 1\\ 1 & 1 \end{bmatrix}^n \begin{bmatrix} f_0\\ f_1 \end{bmatrix}\ (n\in\mathbb{N}). \]

这个\(\begin{bmatrix}0 & 1 \\ 1 & 1\end{bmatrix}^n\)可以用快速幂求出。

只要递推式可以以线性组合的形式表示（不要求一定是普通加法乘法）且所构造出的广义矩阵乘法满足结合律，就可以利用矩乘快速幂计算某一项。有时为了降低常数，我们会手写矩乘而非用循环实现。

数据结构维护广义矩乘

CACC 2024 Regional C题

给定正整数\(n,q\in[1,5\times10^5]\)。一个长为\(n\)的序列\(\vec{a}=(a_0,...,a_{n-1})\)初始均为\(0\)。一共\(q\)次查询\(op\in\{0,1\}\)：

1. \(op=0\)：再给定整数\(b\in[0,n),\ e\in[b,n],odd\in\{0,1\},d\in \mathbb{Z}\)，若\(odd=0\)则将\(a_{[b,e)}\)中的所有偶数均加上\(d\)，否则将其中的所有奇数均加上\(d\)。
2. \(op=1\)：再给定整数\(b\in[0,n),\ e\in[b,n]\)，求出\(\sum\nolimits_{i=b}^{e-1}a_i\)。

对于一个区间，我们如果维护以下三个值，就可以实现题目的修改和查询要求：

• \(c_0\)：区间内偶数的个数。
• \(c_1\)：区间内奇数的个数。
• \(s\)：区间和。

用一个向量

\[\vec{v}= \begin{bmatrix} c_0\\ c_1\\ s \end{bmatrix} \]

保存这三个值。查询时，只需获取区间的\(s\)即为答案。

下面考虑如何将修改操作以矩阵乘法的形式描述：

• 给所有偶数加偶数：

  \[\vec{v}'= \begin{bmatrix} c_0\\ c_1\\ s+d\cdot c_0 \end{bmatrix} = \begin{bmatrix} 1\\ & 1\\ d & & 1 \end{bmatrix} \begin{bmatrix} c_0\\ c_1\\ s \end{bmatrix} \triangleq M_0\vec{v} . \]

• 给所有偶数加奇数：

  \[\vec{v}'= \begin{bmatrix} 0\\ c_0+c_1\\ s+d\cdot c_0 \end{bmatrix} = \begin{bmatrix} 0\\ 1 & 1\\ d & & 1 \end{bmatrix} \begin{bmatrix} c_0\\ c_1\\ s \end{bmatrix} \triangleq M_1\vec{v}. \]

• 给所有奇数加偶数：

  \[\vec{v}'= \begin{bmatrix} c_0\\ c_1\\ s+d\cdot c_1 \end{bmatrix} = \begin{bmatrix} 1\\ & 1\\ & d & 1 \end{bmatrix} \begin{bmatrix} c_0\\ c_1\\ s \end{bmatrix} \triangleq M_2\vec{v} . \]

• 给所有奇数加奇数：

  \[\vec{v}'= \begin{bmatrix} c_0+c_1\\ 0\\ s+d\cdot c_1 \end{bmatrix} = \begin{bmatrix} 1 & 1\\ & 0\\ & d & 1 \end{bmatrix} \begin{bmatrix} c_0\\ c_1\\ s \end{bmatrix} \triangleq M_3\vec{v}. \]

这里是定义在普通加法和普通乘法上的矩阵乘法，显然满足结合律，且有单位元\(I\)，构成幺半群，考虑用线段树维护。

直接这样实现的话，进行一次矩乘需要\(3^3=27\)次乘法，总复杂度\(27\cdot n\cdot q\approx2.56\times10^8\)，卡卡常应该能过（CACC的数据挺水的，而且还是类IOI赛制，狠狠混分）。

struct Matrix {
	std::array<std::array<int64_t, 3>, 3> data{};

	Matrix() {
		data[0][0] = data[1][1] = data[2][2] = 1;
	}

	friend Matrix operator*(const Matrix &lhs, const Matrix &rhs) {
		Matrix res;
		for (size_t i = 0; i < 3; ++i) {
			for (size_t j = 0; j < 3; ++j) {
				res.data[i][j] = 0;
				for (size_t k = 0; k < 3; ++k) {
					res.data[i][j] += lhs.data[i][k] * rhs.data[k][j];
				}
			}
		}
		return res;
	}

	bool isIdentity() const {
		for (size_t i = 0; i < 3; ++i) {
			for (size_t j = 0; j < 3; ++j) {
				if (data[i][j] != (i == j)) return false;
			}
		}
		return true;
	}
};

struct SegTree {
	static size_t leftChild(size_t rt) { return ((rt << 1) | 1); }
	static size_t rightChild(size_t rt) { return ((rt + 1) << 1); }

	struct Node {
		std::array<int64_t, 3> val{};
		Matrix lzy{};

		void applyLazy(const Matrix &lzy) {
			std::array<int64_t, 3> new_val{};
			for (size_t i = 0; i < 3; ++i) {
				for (size_t j = 0; j < 3; ++j) {
					new_val[i] += lzy.data[i][j] * val[j];
				}
			}
			val = new_val;
			this->lzy = lzy * this->lzy;
		}
	};

	std::vector<Node> nodes;

	SegTree(size_t n) : nodes(n << 2) { build(0, 0, n); }

	void pushUp(size_t rt) {
		size_t lch = leftChild(rt), rch = rightChild(rt);
		for (size_t i = 0; i < 3; ++i) {
			nodes[rt].val[i] = nodes[lch].val[i] + nodes[rch].val[i];
		}
	}

	void build(size_t rt, size_t beg, size_t end) {
		assert(beg < end);
		if (beg + 1 == end) {
			nodes[rt].val[0] = 1;
			return;
		}

		size_t mid = beg + ((end - beg) >> 1);
		build(leftChild(rt), beg, mid);
		build(rightChild(rt), mid, end);
		pushUp(rt);
	}

	void pushDown(size_t rt) {
		assert(rightChild(rt) < nodes.size());
		if (nodes[rt].lzy.isIdentity()) return;

		nodes[leftChild(rt)].applyLazy(nodes[rt].lzy);
		nodes[rightChild(rt)].applyLazy(nodes[rt].lzy);

		nodes[rt].lzy = Matrix();
	}

	void modify(size_t rt, size_t nd_beg, size_t nd_end,
				size_t beg, size_t end, bool is_odd, int64_t diff) {
		assert(nd_beg < nd_end && beg < end);

		if (beg <= nd_beg && nd_end <= end) {
			Matrix tag;
			tag.data[2][is_odd] = diff;
			if (diff & 1) {
				if (!is_odd) {
					tag.data[0][0] = 0;
					tag.data[0][1] = 1;
				} else {
					tag.data[1][1] = 0;
					tag.data[1][0] = 1;
				}
			}
			nodes[rt].applyLazy(tag);
			return;
		}

		pushDown(rt);
		size_t nd_mid = nd_beg + ((nd_end - nd_beg) >> 1);
		if (beg < nd_mid) {
			modify(leftChild(rt), nd_beg, nd_mid, beg, end, is_odd, diff);
		}
		if (nd_mid < end) {
			modify(rightChild(rt), nd_mid, nd_end, beg, end, is_odd, diff);
		}
		pushUp(rt);
	}

	int64_t query(size_t rt, size_t nd_beg, size_t nd_end,
				  size_t beg, size_t end) {
		assert(nd_beg < nd_end && beg < end);
		if (beg <= nd_beg && nd_end <= end) return nodes[rt].val[2];

		pushDown(rt);
		size_t nd_mid = nd_beg + ((nd_end - nd_beg) >> 1);
		int64_t res = 0;
		if (beg < nd_mid) res += query(leftChild(rt), nd_beg, nd_mid, beg, end);
		if (nd_mid < end) res += query(rightChild(rt), nd_mid, nd_end, beg, end);

		return res;
	}
};

参考资料

1. Gilbert Strang, Introduction to Linear Algebra, 4th Edition.
2. Sheldon Axler 著, 吴俊达、何阳 译, 线性代数应该这样学, 第四版.
3. 我的同学兼教练（？）Neutralized的课件.
4. OI Wiki.