【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\)。
方阵：行数与列数相等的矩阵称为方阵。
例如：矩阵 \(B\) 为方阵。
单位矩阵：对于方阵 \(I\)，如果主对角线每一个元素均为 \(1\)，其余元素为 \(0\)，这个矩阵被称为单位矩阵，记作 \(I\)。
在代码编写中，我们把矩阵定义为一个结构体，这样方便写在函数中作为参数或返回值。矩阵结构体包含矩阵的二维数组，有时包含矩阵的行数和列数。
struct matrix
{
long long n,m;
long long v[200][200];
};
矩阵运算
矩阵加减
对于两个 \(n\times m\) 的矩阵，它们的和为一个 \(n\times m\) 矩阵，矩阵中各位元素等于对应元素之和。
对于两个 \(n\times m\) 的矩阵，它们的差为一个 \(n\times m\) 矩阵，矩阵中各位元素等于对应元素之差。
例：
矩阵加减法满足实数加减法的性质，比如加法交换律。
矩阵乘法
对于一个 \(p\times q\) 的矩阵 \(A\) 和一个 \(q\times r\) 的矩阵 \(B\)，它们的积为一个 \(p\times r\) 的矩阵 \(C\)。在这个矩阵中，\(c_{i,j}\) 满足如下式子：
例：
矩阵乘法可以简单记为第 \(i\) 行第 \(j\) 列的元素等于第一个矩阵的第 \(i\) 行与第二个矩阵的第 \(j\) 列对应元素相乘的积的和。即先横着看，再竖着看。
矩阵乘法满足结合律，但是不满足交换律，这一点需要特别注意。
单位矩阵 \(I\) 的性质：单位矩阵乘以任何一个矩阵(可以乘的)，积为这个矩阵。单位矩阵在矩阵乘法中，就相当于实数乘法中的 \(1\)。
在代码编写中，我们按照矩阵乘法的定义求出每一个 \(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}\)。我们把这些写入矩阵，并在等式另一边写出每个元素递推一次后的矩阵。
接下来，我们考虑构造中间的转移矩阵。设状态矩阵的第一行为 \(x\;y\;z\)，根据矩阵乘法的运算法则，得到如下式子：
因为等号左边均为 \(f_i\)，我们把这个式子和题目所给的递推式进行对比：
观察系数不难发现 \(x=5,y=0,z=-3\)。所以我们可以填出转移矩阵的第一行：
接下来，设状态矩阵的第二行为 \(x\;y\;z\)，根据矩阵乘法的运算法则，得到如下式子：
很显然，\(x=1,y=0,z=0\) 时就有 \(f_{i-1}=f_{i-1}\)，所以我们可以填出转移矩阵的第二行：
设状态矩阵的第三行为 \(x\;y\;z\)，根据矩阵乘法的运算法则，得到如下式子：
同第二行，\(x=0,y=1,z=0\) 时就有 \(f_{i-2}=f_{i-2}\)，所以我们可以填出转移矩阵的第三行：
这样，我们就成功用矩阵表示了 \(f_i=5f_{i-1}-3f_{i-3}\) 这个递推关系。
这么做的好处是，我们可以用矩阵快速幂在 \(O(\log n)\) 时间内推出 \(f_n\) 的值。我们可以简单推导一下：
使用矩阵快速幂，求出 \(\begin{bmatrix}5&0&-3\\1&0&0\\0&1&0\\\end{bmatrix}^{n-3}\)，再乘以初始状态就能在状态矩阵中推出 \(f_n\)。由于转移矩阵很小，矩阵乘法视为常数，最后总的时间复杂度为 \(O(\log n)\)。
高斯消元法
我们可以用矩阵表示方程组，例如这个例子：
我们只把各项未知数的系数写进矩阵，就得到了这个方程的矩阵表示。例如上面这个方程组，就可以这样表示：
有的时候我们把 \(b_i\) 也写入这个矩阵，我们把这样的矩阵称为增广矩阵。
我们可以使用高斯消元法来求出这个方程的解。高斯消元法具体步骤如下：
假设现在处理到第 \(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\)。所以最后矩阵一定是这个形式，我们称之为上三角形式。
在这个矩阵中，直接从 \(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\)：
矩阵快速幂模板题，不多赘述。
#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\)：
递推需要用到了量有 \(F_{n-1},F_{n-2}\)，把它们写入状态矩阵。
对于每一个量，列出递推关系：
把这些递推关系逐行利用待定系数法写入矩阵，得到用矩阵表示的递推关系：
初始状态 \(\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\)：
我认为这题是变相大模拟。递推式比较复杂，梳理一下，用到的量有 \(a_{k},a_{k+1},b_{k},b_{k+1},c_{k},c_{k+1},k^2,k,1,w^k,z^k\) 共 \(11\) 个，把它们写入状态矩阵。
这里状态矩阵中没有使用 \(w^k,z^k\)，是因为 \(w^{k-1},z^{k-1}\) 通过改变系数可以做到一样的效果。两种写法都可以。
对于每一个量，列出递推关系：
把这些递推关系逐行利用待定系数法写入矩阵，得到用矩阵表示的递推关系：
根据题意，写出 \(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\)：
高斯消元模板题，不多赘述。
#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\)：
由于每次找到的最大值可能为 \(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\)：
设这个球的球心坐标为 \((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})\) 为例，有如下式子：
两边同时平方得：
把平方展开得：
移项合并得：
这样，就产生了一个关于 \((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\)：
我们发现，如果两项颜色相同，则把两项都删去，这很符合异或的性质。再结合后面一条两项颜色不同，将这两项替换为与这两种颜色不同的颜色，我们发现需要找到三个数 \(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;
}
后记
矩阵这个知识点，推导简洁，码量较小，是一个比较简单的知识点。当然，也有可能是我水平有限，没有学习矩阵较难的内容。
不过，确实只有推过莫比乌斯反演和写过高级数据结构，才会发现一个推导简洁，码量较小的知识点确实很难得。
浙公网安备 33010602011771号