【6】矩阵学习笔记

前言

矩阵在 OI 中运用广泛，是一个很重要的内容。正是因为我不会矩阵，所以 WC2024 时讲动态 DP 我没办法听，所以矩阵真的很重要。

由于我实力有限，这里只能介绍一些矩阵的基本应用。

矩阵定义

矩阵：将 $n\times m$ 个数排列成 $n$ 行 $m$ 列的形式称为一个 $n\times m$ 的矩阵。

矩阵一般用大写字母 $A,B,C$ 来表示，矩阵里第 $i$ 行第 $j$ 列的元素用对应的小写字母 $a_{i,j},b_{i,j},c_{i,j}$ 来表示。

例如：在矩阵 $A$ 中，$a_{1,1}=9,a_{3,2}=5,a_{4,3}=9$。

$$A=\begin{bmatrix}9&9&8\\2&4&4\\3&5&3\\9&9&8\end{bmatrix}$$

方阵：行数与列数相等的矩阵称为方阵。

例如：矩阵 $B$ 为方阵。

$$B=\begin{bmatrix}1&2\\3&4\end{bmatrix}$$

单位矩阵：对于方阵 $I$，如果主对角线每一个元素均为 $1$，其余元素为 $0$，这个矩阵被称为单位矩阵，记作 $I$。

$$I_2=\begin{bmatrix}1&0\\0&1\end{bmatrix},I_3=\begin{bmatrix}1&0&0\\0&1&0\\0&0&1\end{bmatrix},I_4=\begin{bmatrix}1&0&0&0\\0&1&0&0\\0&0&1&0\\0&0&0&1\end{bmatrix}\dots$$

在代码编写中，我们把矩阵定义为一个结构体，这样方便写在函数中作为参数或返回值。矩阵结构体包含矩阵的二维数组，有时包含矩阵的行数和列数。

struct matrix
{
   long long n,m;
   long long v[200][200];
};

矩阵运算

矩阵加减

对于两个 $n\times m$ 的矩阵，它们的和为一个 $n\times m$ 矩阵，矩阵中各位元素等于对应元素之和。

对于两个 $n\times m$ 的矩阵，它们的差为一个 $n\times m$ 矩阵，矩阵中各位元素等于对应元素之差。

例：

$$A=\begin{bmatrix}9&9&8\\2&4&4\\3&5&3\\9&9&8\end{bmatrix}, B=\begin{bmatrix}1&4&2\\2&0&4\\1&9&0\\8&1&8\end{bmatrix}$$

$$A+B=\begin{bmatrix}9+1&9+4&8+2\\2+2&4+0&4+4\\3+1&5+9&3+0\\9+8&9+1&8+8\end{bmatrix}=\begin{bmatrix}10&13&10\\4&4&8\\4&14&3\\17&10&16\end{bmatrix}$$

$$A-B=\begin{bmatrix}9-1&9-4&8-2\\2-2&4-0&4-4\\3-1&5-9&3-0\\9-8&9-1&8-8\end{bmatrix}=\begin{bmatrix}8&5&6\\0&4&0\\2&-4&3\\1&8&0\end{bmatrix}$$

矩阵加减法满足实数加减法的性质，比如加法交换律。

矩阵乘法

对于一个 $p\times q$ 的矩阵 $A$ 和一个 $q\times r$ 的矩阵 $B$，它们的积为一个 $p\times r$ 的矩阵 $C$。在这个矩阵中，$c_{i,j}$ 满足如下式子：

$$c_{i,j}=\sum_{k=1}^qa_{i,k}b_{k,j}$$

例：

$$A=\begin{bmatrix}1&2\\2&1\end{bmatrix}, B=\begin{bmatrix}1&3&0\\4&2&5\end{bmatrix}$$

$$AB=\begin{bmatrix}1 \times1+2\times4&1\times3+2\times2 & 1\times0+2\times5\\2\times1+1\times4&2\times3+1\times2&2\times0+1\times5\end{bmatrix}=\begin{bmatrix}9&7&10\\6&8&5\end{bmatrix}$$

矩阵乘法可以简单记为第 $i$ 行第 $j$ 列的元素等于第一个矩阵的第 $i$ 行与第二个矩阵的第 $j$ 列对应元素相乘的积的和。即先横着看，再竖着看。

矩阵乘法满足结合律，但是不满足交换律，这一点需要特别注意。

单位矩阵 $I$ 的性质：单位矩阵乘以任何一个矩阵(可以乘的)，积为这个矩阵。单位矩阵在矩阵乘法中，就相当于实数乘法中的 $1$。

$$IA=A$$

在代码编写中，我们按照矩阵乘法的定义求出每一个 $c_{i,j}$。这个过程时间复杂度较高，为 $O(n^3)$，其中 $n$ 为矩阵的边长(假设 $n,m$ 同阶)。

代码中 $n$ 为方阵边长。如果是矩阵，把循环中的 n 改为对应的 a.n,a.m,b.n,b.m 即可。

struct matrix mul(struct matrix a,struct matrix b)
{
	struct matrix c;
	for(int i=1;i<=n;i++)
	    for(int j=1;j<=n;j++)
	        c.v[i][j]=0;
	for(int i=1;i<=n;i++)
	    for(int j=1;j<=n;j++)
	        for(int k=1;k<=n;k++)
	            c.v[i][j]=c.v[i][j]+a.v[i][k]*b.v[k][j];
	return c;
}

矩阵的幂

对于一个方阵 $A$，记这个方阵的幂为 $A^k$，表示 $k$ 个 $A$ 连乘的积。

特别的，有 $A^0=I$，其中 $I$ 为单位矩阵。

我们都知道，可以用快速幂来求一个实数的 $k$ 次方，那我们也可以用类似的思想来求矩阵的 $k$ 次方，这就是矩阵快速幂。

相较于快速幂模板，矩阵快速幂只需要把实数改为矩阵，把实数运算改为矩阵运算，并将 $ans$ 的初始值设为单位矩阵 $I$ 即可。

时间复杂度为 $O(n^3\log n)$。

代码中 $n$ 为方阵边长。

struct matrix power(struct matrix a,long long p)
{
	struct matrix ans,x=a;
	for(int i=1;i<=n;i++)
	    for(int j=1;j<=n;j++)
	        {
	        if(i!=j)ans.v[i][j]=0;
	        else ans.v[i][j]=1;
	        }
	while(p)
	   {
	   	if(p&1)ans=mul(ans,x);
	   	p>>=1;
	   	x=mul(x,x);
	   }
	return ans;
}

矩阵加速递推

矩阵可以用来表示递推关系。当用矩阵表示递推关系时，我们把状态写入一个 $n\times 1$ 的矩阵，这个矩阵称为状态矩阵；用一个 $n\times n$ 的方阵来表示递推关系，这个方阵称为转移矩阵。

由于矩阵乘法独特的计算方式，我们可以通过转移矩阵里各位置的值来确定状态转移时的系数。当我们转移不需要某一个值时，就把这一位设为 $0$。

例如：用矩阵表示 $f_i=5f_{i-1}-3f_{i-3}$。

首先，我们先把状态写入一个 $n\times 1$ 的矩阵，并在等式另一边写上递推一次后的状态矩阵。这个递推式直接出现的状态为 $f_{i-1},f_{i-3}$，额外会用到的状态为 $f_{i-2}$。我们把这些写入矩阵，并在等式另一边写出每个元素递推一次后的矩阵。

$$\begin{bmatrix}f_i\\f_{i-1}\\f_{i-2}\end{bmatrix}=\begin{bmatrix}\\\\\\\end{bmatrix}\begin{bmatrix}f_{i-1}\\f_{i-2}\\f_{i-3}\end{bmatrix}$$

接下来，我们考虑构造中间的转移矩阵。设状态矩阵的第一行为 $x\;y\;z$，根据矩阵乘法的运算法则，得到如下式子：

$$f_i=xf_{i-1}+yf_{i-2}+zf_{i-3}$$

因为等号左边均为 $f_i$，我们把这个式子和题目所给的递推式进行对比：

$$f_i=xf_{i-1}+yf_{i-2}+zf_{i-3}$$

$$f_i=5f_{i-1}-3f_{i-3}=5f_{i-1}+0f_{i-1}-3f_{i-3}$$

观察系数不难发现 $x=5,y=0,z=-3$。所以我们可以填出转移矩阵的第一行：

$$\begin{bmatrix}f_i\\f_{i-1}\\f_{i-2}\end{bmatrix}=\begin{bmatrix}5&0&-3\\\\\\\end{bmatrix}\begin{bmatrix}f_{i-1}\\f_{i-2}\\f_{i-3}\end{bmatrix}$$

接下来，设状态矩阵的第二行为 $x\;y\;z$，根据矩阵乘法的运算法则，得到如下式子：

$$f_{i-1}=xf_{i-1}+yf_{i-2}+zf_{i-3}$$

很显然，$x=1,y=0,z=0$ 时就有 $f_{i-1}=f_{i-1}$，所以我们可以填出转移矩阵的第二行：

$$\begin{bmatrix}f_i\\f_{i-1}\\f_{i-2}\end{bmatrix}=\begin{bmatrix}5&0&-3\\1&0&0\\\\\end{bmatrix}\begin{bmatrix}f_{i-1}\\f_{i-2}\\f_{i-3}\end{bmatrix}$$

设状态矩阵的第三行为 $x\;y\;z$，根据矩阵乘法的运算法则，得到如下式子：

$$f_{i-2}=xf_{i-1}+yf_{i-2}+zf_{i-3}$$

同第二行，$x=0,y=1,z=0$ 时就有 $f_{i-2}=f_{i-2}$，所以我们可以填出转移矩阵的第三行：

$$\begin{bmatrix}f_i\\f_{i-1}\\f_{i-2}\end{bmatrix}=\begin{bmatrix}5&0&-3\\1&0&0\\0&1&0\\\end{bmatrix}\begin{bmatrix}f_{i-1}\\f_{i-2}\\f_{i-3}\end{bmatrix}$$

这样，我们就成功用矩阵表示了 $f_i=5f_{i-1}-3f_{i-3}$ 这个递推关系。

这么做的好处是，我们可以用矩阵快速幂在 $O(\log n)$ 时间内推出 $f_n$ 的值。我们可以简单推导一下：

$$\begin{aligned} \begin{bmatrix}f_n\\f_{n-1}\\f_{n-2}\end{bmatrix}&=\begin{bmatrix}5&0&-3\\1&0&0\\0&1&0\\\end{bmatrix}\begin{bmatrix}f_{n-1}\\f_{n-2}\\f_{n-3}\end{bmatrix}\\ &=\begin{bmatrix}5&0&-3\\1&0&0\\0&1&0\\\end{bmatrix}^2\begin{bmatrix}f_{n-2}\\f_{n-3}\\f_{n-4}\end{bmatrix}\\ &=\dots\\ &=\begin{bmatrix}5&0&-3\\1&0&0\\0&1&0\\\end{bmatrix}^{n-3}\begin{bmatrix}f_{3}\\f_{2}\\f_{1}\end{bmatrix} \end{aligned}$$

使用矩阵快速幂，求出 $\begin{bmatrix}5&0&-3\\1&0&0\\0&1&0\\\end{bmatrix}^{n-3}$，再乘以初始状态就能在状态矩阵中推出 $f_n$。由于转移矩阵很小，矩阵乘法视为常数，最后总的时间复杂度为 $O(\log n)$。

高斯消元法

我们可以用矩阵表示方程组，例如这个例子：

$$\begin{cases} a_{1, 1} x_1 + a_{1, 2} x_2 + \cdots + a_{1, n} x_n = b_1 \\ a_{2, 1} x_1 + a_{2, 2} x_2 + \cdots + a_{2, n} x_n = b_2 \\ \cdots \\ a_{n,1} x_1 + a_{n, 2} x_2 + \cdots + a_{n, n} x_n = b_n \end{cases}$$

我们只把各项未知数的系数写进矩阵，就得到了这个方程的矩阵表示。例如上面这个方程组，就可以这样表示：

$$\begin{bmatrix}a_{1,1}&a_{1,2}&\dots&a_{1,n}\\a_{2,1}&a_{2,2}&\dots&a_{2,n}\\\dots&\dots&\dots&\dots\\a_{n,1}&a_{n,2}&\dots&a_{n,n}\\\end{bmatrix}$$

有的时候我们把 $b_i$ 也写入这个矩阵，我们把这样的矩阵称为增广矩阵。

$$\begin{bmatrix}a_{1,1}&a_{1,2}&\dots&a_{1,n}&\mid b_1\\a_{2,1}&a_{2,2}&\dots&a_{2,n}&\mid b_2\\\dots&\dots&\dots&\dots&\dots\\a_{n,1}&a_{n,2}&\dots&a_{n,n}&\mid b_n\\\end{bmatrix}$$

我们可以使用高斯消元法来求出这个方程的解。高斯消元法具体步骤如下：

假设现在处理到第 $i$ 行。

$1$：在第 $i\sim n$ 行的第 $i$ 列寻找绝对值最大值，并将绝对值最大值所在的行与第 $i$ 行交换。这么做是为了避免在第 $i$ 行 $x_i$ 的系数为 $0$ 导致的一些奇怪情况。由于前 $i-1$ 行已经处理了，化成了最终需要的形式，所以不能交换。

$2$：化 $x_i$ 的系数为 $1$，把第 $i$ 行的每一个元素(包括 $b_i$)除以 $x_{i,i}$。

$3$：消元。我们希望第 $i+1\sim n$ 行 $x_i$ 的系数为 $0$，所以我们用加减消元的方式进行消元，具体操作见代码。

根据第 $3$ 步，从第 $i+1$ 行开始，$x_i$ 的系数均为 $0$。所以最后矩阵一定是这个形式，我们称之为上三角形式。

$$\begin{bmatrix}1&a_{1,2}&\dots&a_{1,n}\\0&1&\dots&a_{2,n}\\\dots&\dots&\dots&\dots\\0&0&\dots&1&\\\end{bmatrix}$$

在这个矩阵中，直接从 $x_n\to x_1$ 不断代入就可以求出每一个未知数的值。

时间复杂度为 $O(n^3)$。

代码中有一些循环的范围容易记错，需要特别注意。另外，如果某次找到的最大值为 $0$，那么证明这个方程组并不是唯一解，可能是无解或者无穷解。因为如果找到的最大值为 $0$，第 $i$ 行的方程就不会被用上(不然你怎么化系数为 $1$)。根据初中数学，一个 $n$ 元方程组必须至少有 $n$ 个方程才有唯一解，所以这个方程组必定不是唯一解。

bool gauss()
{
	double div=0;
	for(int i=1;i<=n;i++)
	    {
	    long long mx=i;
	    for(int j=i+1;j<=n;j++)
	        if(fabs(a[j][i])>fabs(a[mx][i]))mx=j;
	    if(fabs(a[mx][i])<eps)return 0; 
	    if(i!=mx)
	       {
	       for(int j=1;j<=n;j++)swap(a[i][j],a[mx][j]);
	       swap(b[i],b[mx]);
	       }
	    div=a[i][i],b[i]/=div;
	    for(int j=1;j<=n;j++)a[i][j]/=div;
	    for(int j=i+1;j<=n;j++)
	        {
	        	div=a[j][i],b[j]-=b[i]*div;
	        	for(int k=i;k<=n;k++)a[j][k]-=a[i][k]*div;
			}
	    }
	ans[n]=b[n];
	for(int i=n-1;i>=1;i--)
	    {
	    ans[i]=b[i];
	    for(int j=n;j>i;j--)ans[i]-=a[i][j]*ans[j];
	    }
	return 1;
}

对于高斯消元判定方程组解的情况，具体请看例题 $5$。

这是更具有普适性的高斯消元算法模板。

long long gauss()
{
	double div=0;
	for(int i=1;i<=n;i++)
	    {
	    long long mx=now;
	    for(int j=now+1;j<=n;j++)
	        if(fabs(a[j][i])>fabs(a[mx][i]))mx=j;
	    if(fabs(a[mx][i])<eps)continue; 
	    if(now!=mx)
	       {
	       for(int j=1;j<=n;j++)swap(a[now][j],a[mx][j]);
	       swap(b[now],b[mx]);
	       }
	    div=a[now][i],b[now]/=div;
	    for(int j=1;j<=n;j++)a[now][j]/=div;
	    for(int j=now+1;j<=n;j++)
	        {
	        	div=a[j][i],b[j]-=b[now]*div;
	        	for(int k=i;k<=n;k++)a[j][k]-=a[now][k]*div;
			}
		now++;
	    }
	if(now-1==n)
		{
		ans[n]=b[n];
		for(int i=n-1;i>=1;i--)
		    {
		    ans[i]=b[i];
		    for(int j=n;j>i;j--)ans[i]-=a[i][j]*ans[j];
		    }
		}
	else
	   {
	   	for(int i=now;i<=n;i++)
	   	    if(fabs(b[i])>=eps)return -1;
		return 1;
	   }
	return 0;
}

例题

例题 $1$：

P3390 【模板】矩阵快速幂

矩阵快速幂模板题，不多赘述。

#include <bits/stdc++.h>
using namespace std;
struct matrix
{
	long long v[200][200];
}a;
long long n,k,mod=1000000007;
struct matrix mul(struct matrix a,struct matrix b)
{
	struct matrix c;
	for(int i=1;i<=n;i++)
	    for(int j=1;j<=n;j++)
	        c.v[i][j]=0;
	for(int i=1;i<=n;i++)
	    for(int j=1;j<=n;j++)
	        for(int k=1;k<=n;k++)
	            c.v[i][j]=(c.v[i][j]+a.v[i][k]*b.v[k][j])%mod;
	return c;
}

struct matrix power(struct matrix a,long long p)
{
	struct matrix ans,x=a;
	for(int i=1;i<=n;i++)
	    for(int j=1;j<=n;j++)
	        {
	        if(i!=j)ans.v[i][j]=0;
	        else ans.v[i][j]=1;
	        }
	while(p)
	   {
	   	if(p&1)ans=mul(ans,x);
	   	p>>=1;
	   	x=mul(x,x);
	   }
	return ans;
}

int main()
{
	scanf("%lld%lld",&n,&k);
	for(int i=1;i<=n;i++)
	    for(int j=1;j<=n;j++)
	        scanf("%lld",&a.v[i][j]);
    a=power(a,k);
	for(int i=1;i<=n;i++)
	    {
	    for(int j=1;j<=n;j++)
	        printf("%lld ",a.v[i][j]);
	    printf("\n");
	    }
	return 0;
}

例题 $2$：

P1962 斐波那契数列

递推需要用到了量有 $F_{n-1},F_{n-2}$，把它们写入状态矩阵。

$$\begin{bmatrix}F_i\\F_{i-1}\end{bmatrix}=\begin{bmatrix}\\\\\end{bmatrix}\begin{bmatrix}F_{i-1}\\F_{i-2}\end{bmatrix}$$

对于每一个量，列出递推关系：

$$F_i=F_{i-1}+F_{i-2}$$

$$F_{i-1}=F_{i-1}$$

把这些递推关系逐行利用待定系数法写入矩阵，得到用矩阵表示的递推关系：

$$\begin{bmatrix}F_i\\F_{i-1}\end{bmatrix}=\begin{bmatrix}1&1\\1&0\\\end{bmatrix}\begin{bmatrix}F_{i-1}\\F_{i-2}\end{bmatrix}$$

初始状态 $\begin{bmatrix}F_2\\F_1\end{bmatrix}=\begin{bmatrix}1\\1\end{bmatrix}$，使用矩阵快速幂递推 $n-2$ 次，就可以推出 $F_{n}$。

注意 mul(power(a,n-2),ans) 不能写成 mul(ans,power(a,n-2))，因为矩阵乘法不满足交换律。也要注意 $n=1$ 时的特判，否则快速幂会进入死循环。

#include <bits/stdc++.h>
using namespace std;
struct matrix
{
	long long v[200][200];
}a,ans;
long long n,k,mod=1000000007;
struct matrix mul(struct matrix a,struct matrix b)
{
	struct matrix c;
	for(int i=1;i<=2;i++)
	    for(int j=1;j<=2;j++)
	        c.v[i][j]=0;
	for(int i=1;i<=2;i++)
	    for(int j=1;j<=2;j++)
	        for(int k=1;k<=2;k++)
	            c.v[i][j]=(c.v[i][j]+a.v[i][k]*b.v[k][j])%mod;
	return c;
}

struct matrix power(struct matrix a,long long p)
{
	struct matrix ans,x=a;
	for(int i=1;i<=2;i++)
	    for(int j=1;j<=2;j++)
	        {
	        if(i!=j)ans.v[i][j]=0;
	        else ans.v[i][j]=1;
	        }
	while(p)
	   {
	   	if(p&1)ans=mul(ans,x);
	   	p>>=1;
	   	x=mul(x,x);
	   }
	return ans;
}

int main()
{
	scanf("%lld",&n);
    ans.v[1][1]=1,ans.v[2][1]=1;
	a.v[1][1]=1,a.v[1][2]=1,a.v[2][1]=1,a.v[2][2]=0;
    if(n>1)
        {
        ans=mul(power(a,n-2),ans);
    	printf("%lld",ans.v[1][1]%mod);
        }
    else printf("1\n");
	return 0;
}

例题 $3$：

P1707 刷题比赛

我认为这题是变相大模拟。递推式比较复杂，梳理一下，用到的量有 $a_{k},a_{k+1},b_{k},b_{k+1},c_{k},c_{k+1},k^2,k,1,w^k,z^k$ 共 $11$ 个，把它们写入状态矩阵。

$$\begin{bmatrix}a_{k+1}\\a_{k+2}\\{(k+1)}^2\\k+1\\1\\b_{k+1}\\b_{k+2}\\w^k\\c_{k+1}\\c_{k+2}\\z^k\end{bmatrix}=\begin{bmatrix}\\\\\\\\\\\\\\\\\\\\\\\end{bmatrix}\begin{bmatrix}a_{k}\\a_{k+1}\\k^2\\k\\1\\b_{k}\\b_{k+1}\\w^{k-1}\\c_{k}\\c_{k+1}\\z^{k-1}\end{bmatrix}$$

这里状态矩阵中没有使用 $w^k,z^k$，是因为 $w^{k-1},z^{k-1}$ 通过改变系数可以做到一样的效果。两种写法都可以。

对于每一个量，列出递推关系：

$$a_k=a_k$$

$$a_{k+2}=pa_{k+1}+qa_k+b_{k+1}+c_{k+1}+rk^2+tk+1$$

$$b_k=b_k$$

$$b_{k+2}=ub_{k+1}+vb_k+a_{k+1}+c_{k+1}+w^{k-1}\times w$$

$$c_k=c_k$$

$$c_{k+2} = xc_{k+1}+yc_k + a_{k+1} + b_{k+1} + z^{k-1}\times z+k+2$$

$$(k+1)^2=k^2+2k+1$$

$$k+1=k+1$$

$$1=1$$

$$w^k=w^{k-1}\times w$$

$$z^k=z^{k-1}\times z$$

把这些递推关系逐行利用待定系数法写入矩阵，得到用矩阵表示的递推关系：

$$\begin{bmatrix}a_{k+1}\\a_{k+2}\\{(k+1)}^2\\k+1\\1\\b_{k+1}\\b_{k+2}\\w^k\\c_{k+1}\\c_{k+2}\\z^k\end{bmatrix} = \begin{bmatrix}0&1&0&0&0&0&0&0&0&0&0\\q&p&r&t&1&0&1&0&0&1&0\\0&0&1&2&1&0&0&0&0&0&0\\0&0&0&1&1&0&0&0&0&0&0\\0&0&0&0&1&0&0&0&0&0&0\\0&0&0&0&0&0&1&0&0&0&0\\0&1&0&0&0&v&u&w&0&1&0\\0&0&0&0&0&0&0&w&0&0&0\\0&0&0&0&0&0&0&0&0&1&0\\0&1&0&1&2&0&1&0&y&x&z\\0&0&0&0&0&0&0&0&0&0&z\end{bmatrix} \begin{bmatrix}a_{k}\\a_{k+1}\\k^2\\k\\1\\b_{k}\\b_{k+1}\\w^{k-1}\\c_{k}\\c_{k+1}\\z^{k-1}\end{bmatrix}$$

根据题意，写出 $k=1$ 时等式右边的初始状态。使用矩阵快速幂递推 $n-1$ 次，就可以推出 $a_k,b_k,c_k$。

注意，在递推中有常数时，要把这个常数也写入矩阵状态。因为递推时所有信息都得包含在转移矩阵中，常数也不例外。

由于模数很大，需要使用龟速乘。总体时间复杂度为 $O(\log^2n)$。

#include <bits/stdc++.h>
using namespace std;
struct matrix
{
	long long v[200][200];
}a,ans;
long long n,k,p,q,r,t,u,v,w,x,y,z,mod;
void init()
{
	ans.v[1][1]=1,ans.v[2][1]=3,ans.v[3][1]=1,ans.v[4][1]=1,ans.v[5][1]=1,ans.v[6][1]=1;
	ans.v[7][1]=3,ans.v[8][1]=1,ans.v[9][1]=1,ans.v[10][1]=3,ans.v[11][1]=1;
    a.v[1][2]=1,a.v[2][1]=q,a.v[2][2]=p,a.v[2][3]=r,a.v[2][4]=t,a.v[2][5]=1,a.v[2][7]=1;
    a.v[2][10]=1,a.v[3][3]=1,a.v[3][4]=2,a.v[3][5]=1,a.v[4][4]=1,a.v[4][5]=1,a.v[5][5]=1;
    a.v[6][7]=1,a.v[7][2]=1,a.v[7][6]=v,a.v[7][7]=u,a.v[7][8]=w,a.v[7][10]=1,a.v[8][8]=w;
    a.v[9][10]=1,a.v[10][2]=1,a.v[10][4]=1,a.v[10][5]=2,a.v[10][7]=1,a.v[10][9]=y,a.v[10][10]=x;
    a.v[10][11]=z,a.v[11][11]=z;
}

long long slow(long long a,long long p)
{
	long long ans=0,x=a;
	while(p)
	   {
	   	if(p&1)ans=(ans+x)%mod;
	   	p>>=1;
	   	x=(x+x)%mod;
	   }
	return ans;
}

struct matrix mul(struct matrix a,struct matrix b)
{
	struct matrix c;
	for(int i=1;i<=11;i++)
	    for(int j=1;j<=11;j++)
	        c.v[i][j]=0;
	for(int i=1;i<=11;i++)
	    for(int j=1;j<=11;j++)
	        for(int k=1;k<=11;k++)
	            c.v[i][j]=(c.v[i][j]+slow(a.v[i][k],b.v[k][j]))%mod;
	return c;
}

struct matrix power(struct matrix a,long long p)
{
	struct matrix ans,x=a;
	for(int i=1;i<=11;i++)
	    for(int j=1;j<=11;j++)
	        {
	        if(i!=j)ans.v[i][j]=0;
	        else ans.v[i][j]=1;
	        }
	while(p)
	   {
	   	if(p&1)ans=mul(ans,x);
	   	p>>=1;
	   	x=mul(x,x);
	   }
	return ans;
}

int main()
{
	scanf("%lld%lld",&n,&mod);
	scanf("%lld%lld%lld%lld",&p,&q,&r,&t);
	scanf("%lld%lld%lld",&u,&v,&w);
	scanf("%lld%lld%lld",&x,&y,&z);
    init();
    ans=mul(power(a,n-1),ans);
	printf("nodgd %lld\nCiocio %lld\nNicole %lld\n",ans.v[1][1],ans.v[6][1],ans.v[9][1]);
	return 0;
}

例题 $4$：

P3389 【模板】高斯消元法

高斯消元模板题，不多赘述。

#include <bits/stdc++.h>
using namespace std;
long long n;
double a[200][200],b[200],ans[200],eps=1e-10; 
bool gauss()
{
	double div=0;
	for(int i=1;i<=n;i++)
	    {
	    long long mx=i;
	    for(int j=i+1;j<=n;j++)
	        if(fabs(a[j][i])>fabs(a[mx][i]))mx=j;
	    if(fabs(a[mx][i])<eps)return 0; 
	    if(i!=mx)
	       {
	       for(int j=1;j<=n;j++)swap(a[i][j],a[mx][j]);
	       swap(b[i],b[mx]);
	       }
	    div=a[i][i],b[i]/=div;
	    for(int j=1;j<=n;j++)a[i][j]/=div;
	    for(int j=i+1;j<=n;j++)
	        {
	        	div=a[j][i],b[j]-=b[i]*div;
	        	for(int k=i;k<=n;k++)a[j][k]-=a[i][k]*div;
			}
	    }
	ans[n]=b[n];
	for(int i=n-1;i>=1;i--)
	    {
	    ans[i]=b[i];
	    for(int j=n;j>i;j--)ans[i]-=a[i][j]*ans[j];
	    }
	return 1;
}

int main()
{
	scanf("%lld",&n);
	for(int i=1;i<=n;i++)
	    {
	    for(int j=1;j<=n;j++)
	        scanf("%lf",&a[i][j]);
	    scanf("%lf",&b[i]);
	    }
	if(gauss())
	   {
	   	for(int i=1;i<=n;i++)
	   	    printf("%.2lf\n",ans[i]);
	   }
	else printf("No Solution\n");
	return 0;
}

例题 $5$：

P2455 [SDOI2006] 线性方程组

由于每次找到的最大值可能为 $0$，所以我们分别记录目前处理的未知数 $x_i$ 和目前用到第 $now$ 行的方程。如果找到的绝对值最大值为 $0$，那么我们跳过这个元素，而不是跳过这一行。具体而言，就是把代码中用到 $i$ 作为行数时直接换成 $now$，并在每一次成功操作完一行之后将 $now$ 的值增加。

如果最后我们用到了所有方程，也就是 $now=n+1$，那么这个方程组有唯一解。注意这里是 $n+1$，因为 $now$ 表示的是处理到，但是还没有处理。

否则，这个方程组必定不是唯一解。对于没有用到的方程 $i$ 中未知数 $x_j$ 的系数 $x_{i,j}$，要么整个第 $j$ 列绝对值最大值为 $0$，要么被消元消成 $0$，要么被换走再消元消成 $0$，所以这个方程的未知数系数一定全为 $0$。

我们遍历未被处理的第 $i$ 个方程，根据上文，这个方程的未知数系数一定全为 $0$，所以等号左边一定为 $0$。如果 $b_i\ne0$，就是要求 $0\ne0$，这个方程组显然无解。否则，方程组有解，且由于没有用足 $n$ 个方程，一定无法确定某些未知数，所以方程组有无穷解。

#include <bits/stdc++.h>
using namespace std;
long long n,now=1;
double a[200][200],b[200],ans[200],eps=1e-10; 
long long gauss()
{
	double div=0;
	for(int i=1;i<=n;i++)
	    {
	    long long mx=now;
	    for(int j=now+1;j<=n;j++)
	        if(fabs(a[j][i])>fabs(a[mx][i]))mx=j;
	    if(fabs(a[mx][i])<eps)continue; 
	    if(now!=mx)
	       {
	       for(int j=1;j<=n;j++)swap(a[now][j],a[mx][j]);
	       swap(b[now],b[mx]);
	       }
	    div=a[now][i],b[now]/=div;
	    for(int j=1;j<=n;j++)a[now][j]/=div;
	    for(int j=now+1;j<=n;j++)
	        {
	        	div=a[j][i],b[j]-=b[now]*div;
	        	for(int k=i;k<=n;k++)a[j][k]-=a[now][k]*div;
			}
		now++;
	    }
	if(now-1==n)
		{
		ans[n]=b[n];
		for(int i=n-1;i>=1;i--)
		    {
		    ans[i]=b[i];
		    for(int j=n;j>i;j--)ans[i]-=a[i][j]*ans[j];
		    }
		}
	else
	   {
	   	for(int i=now;i<=n;i++)
	   	    if(fabs(b[i])>=eps)return -1;
		return 1;
	   }
	return 0;
}

int main()
{
	scanf("%lld",&n);
	for(int i=1;i<=n;i++)
	    {
	    for(int j=1;j<=n;j++)
	        scanf("%lf",&a[i][j]);
	    scanf("%lf",&b[i]);
	    }
	long long flag=gauss();
	if(flag==0)
	   {
	   	for(int i=1;i<=n;i++)
	   	    printf("x%d=%.2lf\n",i,ans[i]);
	   }
	else if(flag==1)printf("0\n");
	else printf("-1\n");
	return 0;
}

例题 $6$：

P4035 [JSOI2008] 球形空间产生器

设这个球的球心坐标为 $(x_1,x_2\dots x_n)$，根据球心到球面上任意一点距离都相等的，我们可以列出一些等量关系式，组成一个方程组。

以第一个坐标 $(a_{1,1},a_{1,2}\dots x_{1,n})$ 和第二个坐标为 $(a_{2,1},a_{2,2}\dots x_{2,n})$ 为例，有如下式子：

$$\sqrt{(a_{1,1}-x_1)^2 + (a_{1,2}-x_2)^2 + \dots + (a_{1,n}-x_n)^2 }=\sqrt{(a_{2,1}-x_1)^2 + (a_{2,2}-x_2)^2 + \dots + (a_{2,n}-x_n)^2 }$$

两边同时平方得：

$$(a_{1,1}-x_1)^2 + (a_{1,2}-x_2)^2 + \dots + (a_{1,n}-x_n)^2 =(a_{2,1}-x_1)^2 + (a_{2,2}-x_2)^2 + \dots + (a_{2,n}-x_n)^2$$

把平方展开得：

$$a_{1,1}^2-2a_{1,1}x_1+x_1^2 + a_{1,2}^2-2a_{1,2}x_2+x_2^2 + \dots + a_{1,n}^2-2a_{1,n}x_n+x_n^2 =a_{2,1}^2-2a_{2,1}x_1+x_1^2 + a_{2,2}^2-2a_{2,2}x_2+x_2^2 + \dots + a_{2,n}^2-2a_{2,n}x_n+x_n^2$$

移项合并得：

$$2(a_{1,1}-a_{2,1})x_1+2(a_{1,2}-a_{2,2})x_2\dots2(a_{1,n}-a_{2,n})x_n=a_{1,1}^2+a_{1,2}^2+\dots a_{1,n}^2-a_{2,1}^2+a_{2,2}^2+\dots a_{2,n}^2$$

这样，就产生了一个关于 $(x_1,x_2\dots x_n)$ 的方程，且未知数的系数和等号右边的数都已知。对于任意两个的点的坐标，我们都可以这样写出一个方程。只要我们写出至少 $n$ 个这样的方程联立，就可以用高斯消元法求出每一个未知数，进而确定球心的坐标。

为了方便，对于第 $i(i\ne1)$ 个点，我们都利用这个点和第 $i-1$ 个点到球心的距离建立方程。这样我们就可以在输入时记录上一次输入的值与平方和，简便地建立方程。一共输入 $n+1$ 个点，这样刚好建立了 $n$ 个方程。

#include <bits/stdc++.h>
using namespace std;
long long n;
double a[200][200],b[200],ans[200],l[200],now[200],ls,ns,eps=1e-10; 
bool gauss()
{
	double div=0;
	for(int i=1;i<=n;i++)
	    {
	    long long mx=i;
	    for(int j=i+1;j<=n;j++)
	        if(fabs(a[j][i])>fabs(a[mx][i]))mx=j;
	    if(fabs(a[mx][i])<eps)return 0; 
	    if(i!=mx)
	       {
	       for(int j=1;j<=n;j++)swap(a[i][j],a[mx][j]);
	       swap(b[i],b[mx]);
	       }
	    div=a[i][i],b[i]/=div;
	    for(int j=1;j<=n;j++)a[i][j]/=div;
	    for(int j=i+1;j<=n;j++)
	        {
	        	div=a[j][i],b[j]-=b[i]*div;
	        	for(int k=i;k<=n;k++)a[j][k]-=a[i][k]*div;
			}
	    }
	ans[n]=b[n];
	for(int i=n-1;i>=1;i--)
	    {
	    ans[i]=b[i];
	    for(int j=n;j>i;j--)ans[i]-=a[i][j]*ans[j];
	    }
	return 1;
}

int main()
{
	scanf("%lld",&n);
	for(int i=1;i<=n+1;i++)
	    {
	    ns=0;
	    for(int j=1;j<=n;j++)
	        {
	        scanf("%lf",&now[j]);
	        ns+=now[j]*now[j];
	        if(i!=1)a[i-1][j]=2*(l[j]-now[j]);
	        }
	    if(i!=1)b[i-1]=ls-ns;
	    for(int j=1;j<=n;j++)l[j]=now[j];
	    ls=ns;
	    }
	gauss();
	for(int i=1;i<=n;i++)
	   	printf("%.3lf ",ans[i]);
	return 0;
}

例题 $7$：

CF1344F Piet's Palette

我们发现，如果两项颜色相同，则把两项都删去，这很符合异或的性质。再结合后面一条两项颜色不同，将这两项替换为与这两种颜色不同的颜色，我们发现需要找到三个数 $a,b,c$ 来代表三种颜色，并且满足 $a\oplus b=c,a\oplus c=b,b\oplus c=a$。当 $a=1,b=2,c=3$ 时，刚好满足条件。自然，$0$ 就代表了白色与空位。

接下来，我们考虑如何维护另外三种操作。每一个位置的颜色需要两个二进制位来存储，我们把每一个位置的两个二进制位拆成 $a_i$ 和 $b_i$，其中 $a_i$ 表示低位，$b_i$ 表示高位。我们令 $1$ 对应红色，$2$ 对应黄色，$3$ 对应蓝色。

对于 RY 操作，我们只需要把子序列内 $a_i$ 和 $b_i$ 交换，就把所有 R 变为 Y，所有 Y 变为 R，B 和空位置不变。

修改前颜色	修改前编号	修改后编号	修改后颜色
R	$01$	$10$	Y
Y	$10$	$01$	R
B	$11$	$11$	B
.	$00$	$00$	.

对于 RB 操作，我们只需要把子序列内 $b_i$ 改为 $a_i\oplus b_i$，就把所有 R 变为 B，所有 B 变为 R，Y 和空位置不变。

修改前颜色	修改前编号	修改后编号	修改后颜色
R	$01$	$11$	B
Y	$10$	$10$	Y
B	$11$	$01$	R
.	$00$	$00$	.

对于 YB 操作，我们只需要把子序列内 $a_i$ 改为 $a_i\oplus b_i$，就把所有 Y 变为 B，所有 B 变为 Y，R 和空位置不变。

修改前颜色	修改前编号	修改后编号	修改后颜色
R	$01$	$01$	R
Y	$10$	$11$	B
B	$11$	$10$	Y
.	$00$	$00$	.

我们需要对于每个位置维护两个变量，$sa_i$ 和 $sb_i$，分别表示低位是 $a_i,b_i,a_i\oplus b_i$ 中的哪一个和高位是 $a_i,b_i,a_i\oplus b_i$ 中的哪一个。我们利用位运算，来记录构成这一位的答案中是否存在 $a_i$ 或 $b_i$ 。

对于每一个 mix 操作，把子序列中的低位和高位依次异或，结果要等于给出的混合结果对应的编号的低位和高位，我们就能得到两个等式。于是，我们就得到了一个异或方程组。直接使用高斯消元算法求解即可。

#include <bits/stdc++.h>
using namespace std;
long long n,k,now[4100],y[4100],sa[4100],sb[4100],cnt=0;
bool a[4100][4100],b[4100],ans[4100];
string op;
void ry()
{
	long long m=0,now=0;
	cin>>m;
	for(int i=1;i<=m;i++)
	    {
	    cin>>now;
	    swap(sa[now],sb[now]);
	    }
}

void rb()
{
	long long m=0,now=0;
	cin>>m;
	for(int i=1;i<=m;i++)
		{
		cin>>now;
		sb[now]=sa[now]^sb[now];
	    }
}

void yb()
{
	long long m=0,now=0;
	cin>>m;
	for(int i=1;i<=m;i++)
	    {
	    cin>>now;
	    sa[now]=sa[now]^sb[now];
	    }
}

void mix()
{
	long long m=0,id=0;
	char col=0;
	cin>>m;
	for(int i=1;i<=m;i++)cin>>now[i];
	cin>>col;
	if(col=='W')id=0;
	if(col=='R')id=1;
	if(col=='Y')id=2;
	if(col=='B')id=3;
	cnt++,b[cnt]=id&1;
	for(int i=1;i<=m;i++)
	    {
	    if(sa[now[i]]&1)a[cnt][now[i]]=1;
	    if(sa[now[i]]&2)a[cnt][now[i]+n]=1;
	    }
	cnt++,b[cnt]=(id>>1)&1;
	for(int i=1;i<=m;i++)
	    {
	    if(sb[now[i]]&1)a[cnt][now[i]]=1;
	    if(sb[now[i]]&2)a[cnt][now[i]+n]=1;
	    }
}

bool gauss()
{
	long long now=1,res=0;
    for(int i=1;i<=n*2;i++)
        {
        long long mx=now;
        for(int j=now;j<=cnt;j++)
            if(a[j][i]!=0)
               {
			   mx=j,res=max(res,mx);
			   break;
		       }
        if(a[mx][i]==0)continue;
        if(now!=mx)
           {
           for(int j=1;j<=n*2;j++)swap(a[now][j],a[mx][j]);
           swap(b[now],b[mx]);
           }
        for(int j=now+1;j<=cnt;j++)
            if(a[j][i]==1)
	            {
	                b[j]^=b[now];
	                for(int k=i;k<=n*2;k++)a[j][k]^=a[now][k];
	            }
	    y[now]=i;
	    now++;
		}
	for(int i=now;i<=cnt;i++)
	    if(b[i])return 0;
    return 1;
}

void print(long long p)
{
	if(ans[p]==0&&ans[p+n]==0)printf(".");
	if(ans[p]==1&&ans[p+n]==0)printf("R");
	if(ans[p]==0&&ans[p+n]==1)printf("Y");
	if(ans[p]==1&&ans[p+n]==1)printf("B");
}

int main()
{
    cin>>n>>k;
    for(int i=1;i<=n;i++)sa[i]=1,sb[i]=2;
    for(int i=1;i<=k;i++)
        {
        	cin>>op;
        	if(op=="RY")ry();
        	else if(op=="RB")rb();
        	else if(op=="YB")yb();
        	else if(op=="mix")mix();
		}
	if(!gauss())printf("NO\n");
	else
	   {
	   	printf("YES\n");
	   	ans[y[n*2]]=b[n*2];
	   	for(int i=n*2-1;i>=1;i--)
	        {
	        ans[y[i]]=b[i];
	        for(int j=n*2;j>y[i];j--)
			    if(a[i][j])ans[y[i]]^=ans[j];
	        }
	    for(int i=1;i<=n;i++)print(i);
	    printf("\n");
	   }
	return 0;
}

后记

矩阵这个知识点，推导简洁，码量较小，是一个比较简单的知识点。当然，也有可能是我水平有限，没有学习矩阵较难的内容。

不过，确实只有推过莫比乌斯反演和写过高级数据结构，才会发现一个推导简洁，码量较小的知识点确实很难得。