线性代数是 OI 数学知识中的一个重要板块,主要分为矩阵,高斯消元,线性基三个部分,接下来主要介绍这三个板块的知识。

矩阵

定义

对于矩阵(Matrix\mathrm{Matrix}),我们可以形象地把它理解成一个二维数组。由 m×nm\times n 个数 ai,ja _{i,j} 排成的 mm 行 nn 列的数表称为 mm 行 nn 列的矩阵,我们把这个矩阵设为 AA,则其可以如下这样表示。

A=[a1,1a1,2⋯a1,na2,1a2,2⋯a2,n⋮⋮⋱⋮am,1am,2⋯am,n] A=\begin{bmatrix} a _{1,1} & a _{1,2} & \cdots & a _{1,n}\\ a _{2,1} & a _{2,2} & \cdots & a _{2,n}\\ \vdots & \vdots & \ddots & \vdots\\ a _{m,1} & a _{m,2} & \cdots & a _{m,n} \end{bmatrix}

对于一个 n×nn\times n 的方阵 AA,主对角线是指所有 Ai,i(1≤i≤n)A _{i,i}(1\le i\le n) 构成的对角线。对于主对角线上全为 11,其余位置为 00 的矩阵,我们称作 nn 阶单位矩阵,用 II 表示。

对于一个方阵 AA,如果存在矩阵 BB 使得 A×B=IA\times B=I,则 BB 是 AA 的逆矩阵,可写作 B=A−1B=A ^{-1}。矩阵的逆是唯一的。我们假设 AB=I,CA=IAB=I,CA=I,可以推出:

B=IB=(CA)B=C(AB)=CI=C \begin{aligned} B & =IB=(CA)B\\ & =C(AB)=CI=C \end{aligned}

矩阵运算

矩阵加减

矩阵的加减法只对两个相同规模的矩阵成立,即把其对应的元素相加减。设 A,BA,B 为两个 m×nm\times n 的矩阵,C=A±BC=A\pm B,那么 CC 也是一个 m×nm\times n 的矩阵,且满足:

Ci,j=Ai,j±Bi,j C _{i,j}=A _{i,j}\pm B _{i,j}

矩阵加法满足交换律和结合律,即:
A+B=B+AA+(B+C)=(A+B)+C \begin{aligned} & A+B=B+A\\ & A+(B+C)=(A+B)+C \end{aligned}

矩阵乘法

对于一个数 λ\lambda 乘以一个矩阵,直接把这个数和矩阵的每一个元素相乘即可。对于两个矩阵相乘,只有在第一个矩阵的列数和第二个矩阵的行数相同时才有意义。设 AA 为一个 m×nm\times n 的矩阵,BB 为一个 n×pn\times p 的矩阵,C=A×BC=A\times B,那么 CC 为一个 m×pm\times p 的矩阵,且满足:

Ci,j=∑k=1nAi,k×Bk,j C _{i,j}=\sum _{k=1} ^{n} A _{i,k}\times B _{k,j}

矩阵乘法满足结合律和分配律(对数和矩阵都满足),但不满足交换律(只对矩阵),即
(pq)A=p(qA)(BC)A=B(CA)(p+q)A=pA+qAp(A+B)=pA+pBAB≠BA \begin{aligned} & (pq)A=p(qA) & (BC)A=B(CA)\\ & (p+q)A=pA+qA & p(A+B)=pA+pB\\ & AB\neq BA \end{aligned}

同理地,我们可以用快速幂地思想优化矩阵乘法。由于线性递推式可以表示成矩阵乘法地形式,以通常用矩阵快速幂加速递推来求得递推数列的某项。

两个矩阵相乘,我们还有另外一种理解方式。我们设 A×B=CA\times B=C,如下所示:

[010101100]×[123456789]=[45681012123] \begin{bmatrix} 0 & 1 & 0\\ 1 & 0 & 1\\ 1 & 0 & 0\\ \end{bmatrix}\times \begin{bmatrix} 1 & 2 & 3\\ 4 & 5 & 6\\ 7 & 8 & 9\\ \end{bmatrix}= \begin{bmatrix} 4 & 5 & 6\\ 8 & 10 & 12\\ 1 & 2 & 3\\ \end{bmatrix}

我们把 BB 的每一行看作一个向量,那么 CC 矩阵的每一行就是 BB 的 33 个行向量的一个线性组合。


线性组合

即对一个向量集合 (x,y,z)(x,y,z),以这个向量集合作为基张成的向量空间中的任何一个元素都是这个向量集合的一个线性组合,其实就相当于对每个向量配一个实数系数相加,ax+by+czax+by+cz 为线性组合。


例如对于 CC 的第一行 (4,5,6)(4,5,6),它其实是以 AA 矩阵的第一行作为系数对 BB 矩阵的 33 个行向量的一个线性组合,即 (4,5,6)=0×(1,2,3)+1×(4,5,6)+0×(7,8,9)(4,5,6)=0\times (1,2,3)+1\times (4,5,6)+0\times (7,8,9),同理可得:

(8,10,12)=1×(1,2,3)+0×(4,5,6)+1×(7,8,9)(1,2,3)=1×(1,2,3)+0×(4,5,6)+0×(7,8,9) \begin{aligned} & (8,10,12)=1\times (1,2,3)+0\times (4,5,6)+1\times (7,8,9)\\ & (1,2,3)=1\times (1,2,3)+0\times (4,5,6)+0\times (7,8,9) \end{aligned}

高斯消元

高斯消元主要用于解决解线性方程组的方程。对于一个 nn 个方程 nn 个未知数的方程组,高斯消元可以在 O(n3)O(n ^{3}) 的时间内解出每一个未知数。

首先先要了解线性方程组,它是由 mm 个 nn 元一次方程组构成的。线性方程组的所有系数可以写成一个 mm 行 nn 列的系数矩阵,再加上每个方程右侧的常数,可以写成一个 mm 行 n+1n+1 列的增广矩阵,如下所示:

{x+2y−z=−62x+y−3z=−9−x−y+2z=7⇒[12−121−3−1−12−6−97] \begin{cases} x+2y-z=-6\\ 2x+y-3z=-9\\ -x-y+2z=7 \end{cases}\Rightarrow \left[\begin{array}{c|c} \begin{matrix} 1 & 2 & -1\\ 2 & 1 & -3\\ -1 & -1 & 2\\ \end{matrix} & \begin{matrix} -6\\ -9\\ 7 \end{matrix} \end{array}\right]

要求解这个方程组需要对增广矩阵进行一个叫做初等行变换的操作。


初等行变换

矩阵的初等行变换分为 33 类:

  1. 交换两行;
  2. 用一个非零的数乘以某一行;
  3. 把其中一行的若干倍加到另一行上。

对一个矩阵进行初等行变换相当于左乘一个系数矩阵,可通过向量线性组合的方式来理解。


我们通过若干次初等行变换的操作后,可以将增广矩阵化成一个阶梯形矩阵,系数矩阵部分为一个上三角矩阵,如下图所示:

[12−1011001−613] \left[\begin{array}{c|c} \begin{matrix} 1 & 2 & -1\\ 0 & 1 & 1\\ 0 & 0 & 1\\ \end{matrix} & \begin{matrix} -6\\ 1\\ 3 \end{matrix} \end{array}\right]

此时我们已经知道了最后一个方程未知数的解,依次向上回代即可。

最后会得到如下图所示的一个矩阵:

[1000100011−23] \left[\begin{array}{c|c} \begin{matrix} 1 & 0 & 0\\ 0 & 1 & 0\\ 0 & 0 & 1\\ \end{matrix} & \begin{matrix} 1\\ -2\\ 3 \end{matrix} \end{array}\right]

此时系数矩阵部分是一个对角矩阵,说明我们已经求出了方程组的解。对于第 ii 个方程,我们把最后消元结束后,这一行得到的未知数称为这个方程的主元,一般第 ii 个未知数就是第 ii 个方程的主元。通过初等行变换先将系数矩阵消成上三角矩阵,再通过回代将上三角矩阵消成对角矩阵的算法就是高斯消元算法。

总结一下,高斯消元算法分为如下步骤:

  1. (找到在这一行后面这一列系数最大且系数大于 00 的方程,与这一行交换)(可以省略,目的是从系数最大的开始消元,会使精度得到提升);
  2. 从这一行来消后面每一行的方程,保证能把这一列的系数消成 00;
  3. 最后从后往前,每次把这一行除主元其他的系数消完,得到这一个方程的解。

上面给出的例子为比较理想的情况,但事实上在高斯消元中可能遇到不同的情况,比如如下所示:

[12−1011000−613] \left[\begin{array}{c|c} \begin{matrix} 1 & 2 & -1\\ 0 & 1 & 1\\ 0 & 0 & 0\\ \end{matrix} & \begin{matrix} -6\\ 1\\ 3 \end{matrix} \end{array}\right]

发现在消元过程中遇到了一个方程系数全为 00,但常数项不为 00,说明某些方程存在矛盾,此方程组无解。

另一种情况如下图所示:

[12−1011000−610] \left[\begin{array}{c|c} \begin{matrix} 1 & 2 & -1\\ 0 & 1 & 1\\ 0 & 0 & 0\\ \end{matrix} & \begin{matrix} -6\\ 1\\ 0 \end{matrix} \end{array}\right]

此时最后一个方程是一个不定方程,它没法解出任何一个未知数的解,这会使得方程组出现了方程个数小于未知数个数的情况,即无穷多解的情况。我们把最后一行这种不定方程的主元称为自由元。

对于这两种情况,我们只用在高斯消元消上三角矩阵结束后判断即可。由于我们是每次找系数大于 00 的方程交换,如果后面出现这种系数全为 00 的方程,那么我们高斯消元结束后枚举的行会不到 mm 就结束。我们把它后面到 mm 的方程依次判断一下常数项是否存在不为 00 的情况,如果不为 00 说明无解,如果不存在说明无穷解。

例 1

模板题。

Luogu3389 高斯消元法

Luogu2455 线性方程组

#include<cstdio>
#include<cstring>
#include<cstdlib>
#include<cmath>
#include<algorithm>
#define INF 1e9
using namespace std;
const int maxn=110;
int n;
double map[maxn][maxn],ans[maxn];
inline bool gauss()
{
	int i=1,j=1;
	while(i<=n&&j<=n)
	{
		int maxx=i;
		for(int k=i+1;k<=n;k++)
			if(fabs(map[k][j])>fabs(map[maxx][j]))
				maxx=k;
		if(fabs(map[maxx][j])>0)
		{
			swap(map[i],map[maxx]);
			for(int k=i+1;k<=n;k++)
			{
				double q=map[k][j]/map[i][j];
				for(int l=j;l<=n+1;l++)
					map[k][l]-=q*map[i][l];
			}
			i++;
		}
		j++;
	}
	if(i<=n)
		return false;
	for(i=n;i>=1;i--)
	{
		for(j=i+1;j<=n;j++)
			map[i][n+1]-=ans[j]*map[i][j];
		ans[i]=map[i][n+1]/map[i][i];
	}
	return true;
}
int main()
{
	scanf("%d",&n);
	for(int i=1;i<=n;i++)
		for(int j=1;j<=n+1;j++)
			scanf("%lf",&map[i][j]);
	bool flag=gauss();
	if(!flag)
		printf("No Solution\n");
	else
		for(int i=1;i<=n;i++)
			printf("%.2lf\n",ans[i]);
}
#include<cstdio>
#include<cstring>
#include<cstdlib>
#include<cmath>
#include<algorithm>
#define INF 1e9
using namespace std;
const int maxn=55;
int n;
double map[maxn][maxn],ans[maxn],eps=1e-5;
inline int gauss()
{
	int i=1,j=1;
	while(i<=n&&j<=n)
	{
		int maxx=i;
		for(int k=i+1;k<=n;k++)
			if(fabs(map[k][j])>fabs(map[maxx][j]))
				maxx=k;//找到系数最大的行
		if(fabs(map[maxx][j])>eps)//是否大于0
		{
			swap(map[i],map[maxx]);//交换
			for(int k=i+1;k<=n;k++)//后面的行
			{
				double q=map[k][j]/map[i][j];
				for(int l=j;l<=n+1;l++)
					map[k][l]-=q*map[i][l];//消去后面每一行
			}
			i++;//消下一行,只有这一行这一列大于0
		}
		j++;//消下一列
	}
	for(int k=i;k<=n;k++)
		if(fabs(map[k][n+1])>eps)//后面常数存在大于0的情况
			return -1;//无解
	if(i<=n)//不存在大于0的情况
		return 0;//多解
	for(i=n;i>=1;i--)//倒着回代
	{
		for(j=i+1;j<=n;j++)
			map[i][n+1]-=ans[j]*map[i][j];//把这个方程其他系数消掉
		ans[i]=map[i][n+1]/map[i][i];//得到这个方程主元的解
	}
	return 1;
}
int main()
{
	scanf("%d",&n);
	for(int i=1;i<=n;i++)
		for(int j=1;j<=n+1;j++)
			scanf("%lf",&map[i][j]);
	int flag=gauss();
	if(flag!=1)
		printf("%d\n",flag);
	else
		for(int i=1;i<=n;i++)
			printf("x%d=%.2lf\n",i,ans[i]+eps);
}

例 2

矩阵求逆。

Luogu3389 矩阵求逆

解析:对于一个矩阵 AA,我们要求它的逆矩阵 A−1A ^{-1}。由于我们知道初等行变换相当于左乘一个矩阵,我们发现将矩阵 AA 高斯消元后就变成了 II,相当于乘了一个 A−1A ^{-1}。那么我们维护两个矩阵,第一个矩阵为 AA,第二个矩阵为 II,对 AA 高斯消元,同时对 II 做与 AA 同样的初等行变换操作,因此 II 也会乘以 A−1A ^{-1} 变成 A−1A ^{-1},最后将第二个矩阵输出即可。注意到这是在模意义下进行高斯消元,那么就将高斯消元里面的所有运算都取模,把除法改成乘以逆元即可。同时还要注意无解的情况,一个矩阵不存在逆矩阵,说明高斯消元不能成功进行,即出现了无解或多解的情况,其实就相当于系数矩阵最后存在全为 00 的行。

/*
 * @Author: clorf 
 * @Date: 2020-09-09 15:32:27 
 * @Last Modified by: clorf
 * @Last Modified time: 2020-09-14 09:12:28
 */
#include<cstdio>
#include<cstring>
#include<cstdlib>
#include<cmath>
#include<algorithm>
#include<ctime>
#define INF 1e9
using namespace std;
const int maxn=410;
const int mod=1e9+7;
const double Pi=acos(-1.0);
// 数组一定要算好开!!!!边和点的关系
// 记得看开不开long long!!!
template<class T>void read(T &x)
{
    x=0;int f=0;char ch=getchar();
    while(ch<'0'||ch>'9') {f|=(ch=='-');ch=getchar();}
    while(ch>='0'&&ch<='9'){x=(x<<1)+(x<<3)+(ch^48);ch=getchar();}
    x=f?-x:x;
    return;
}
// n和m不要搞混!!!任何循环时候都要注意!!
int n;
long long a[maxn][maxn],b[maxn][maxn];
long long power(long long a,long long b,int p)
{
    long long ans=1%p;
    for(;b;b>>=1)
    {
        if(b&1)
            ans=ans*a%p;
        a=a*a%p;
    }
    return ans%p;
}
long long inv(long long x)
{
    return power(x,mod-2,mod);
}
bool gauss()
{
    int i=1,j=1;
    while(i<=n&&j<=n)
    {
        int maxx=i;
        for(int k=i+1;k<=n;k++)
            if(abs(a[k][j])>abs(a[maxx][j]))
                maxx=k;
        if(abs(a[maxx][j]))
        {
            swap(a[maxx],a[i]);
            swap(b[maxx],b[i]);
            for(int k=i+1;k<=n;k++)
            {
                long long num=(a[k][j]*inv(a[i][j]))%mod;
                for(int l=1;l<=n;l++)
                {
                    a[k][l]=(a[k][l]+(mod-num)*a[i][l])%mod;
                    b[k][l]=(b[k][l]+(mod-num)*b[i][l])%mod;
                }
            }
            i++;
        }
        j++;
    }
    if(i<=n)
        return 0;
    for(int k=n;k>=1;k--)
    {
        for(int p=1;p<=n;p++)
            b[k][p]=(b[k][p]*inv(a[k][k]))%mod;
        a[k][k]=1;
        for(int l=1;l<k;l++)
        {
            for(int p=1;p<=n;p++)
                b[l][p]=(b[l][p]+(mod-a[l][k])*b[k][p])%mod;
            a[l][k]=0;
        }
    }
    return 1;
}
int main()
{
    scanf("%d",&n);
    for(int i=1;i<=n;i++)
    {
        b[i][i]=1;
        for(int j=1;j<=n;j++)
            scanf("%lld",&a[i][j]);
    }
    bool flag=gauss();
    if(!flag)
        printf("No Solution");
    else
        for(int i=1;i<=n;i++)
        {
            for(int j=1;j<=n;j++)
                printf("%lld ",b[i][j]);
            printf("\n");
        }
    return 0;
}

例 3

给出一个 nn 维球体上的 n+1n + 1 个点,求球心坐标。

Luogu4035 球形空间产生器

解析:球体上所有点到球心的距离相等,设球心坐标为 (x1,x2,⋯ ,xn)(x _{1},x _{2},\cdots,x _{n}),满足对任意一个 i∈[1,n+1]i\in [1,n+1]:

∑j=1n(ai,j−xj)2=k \sum _{j=1} ^{n}(a _{i,j}-x _{j}) ^{2}=k

其中 kk 是常数,第 ii 个点的坐标为 (ai,1,ai,2,⋯ ,ai,n)(a _{i,1},a _{i,2},\cdots,a _{i,n})。我们为了把 kk 消掉,将 n+1n+1 个方程两两作差,得到 nn 个 nn 元一次方程:
∑j=1n(ai,j2−ai+1,j2−2xj(ai,j−ai+1,j))=0∑j=1n2(ai,j−ai+1,j)xj=∑j=1n(ai,j2−ai+1,j2) \begin{aligned} & \sum _{j=1} ^{n}(a _{i,j} ^{2}-a _{i+1,j} ^{2}-2x _{j}(a _{i,j}-a _{i+1,j}))=0\\ & \sum _{j=1} ^{n}2(a _{i,j}-a _{i+1,j})x _{j}=\sum _{j=1} ^{n}(a _{i,j} ^{2}-a _{i+1,j} ^{2}) \end{aligned}

接着高斯消元即可。

#include<cstdio>
#include<cmath>
#include<cstring>
#include<cstdlib>
#include<algorithm>
using namespace std;
int n;
double a[15][15],map[21][21];
double ans[21];
void gauss()
{
	for(int i=1;i<=n;i++)
	{
		int maxx=i;
		for(int j=i+1;j<=n;j++)
			if(fabs(map[j][i])>fabs(map[maxx][i]))
				maxx=j;
		if(maxx!=i)
			swap(map[maxx],map[i]);
		double num=map[i][i];
		for(int j=i+1;j<=n+1;j++)
			map[i][j]/=num;
		for(int j=i+1;j<=n;j++)
		{
			double num=map[j][i];
			for(int k=i;k<=n+1;k++)
				map[j][k]-=map[i][k]*num;
		}
	}
	ans[n]=map[n][n+1];
	for(int i=n-1;i>=1;i--)
	{
		ans[i]=map[i][n+1];
		for(int j=i+1;j<=n;j++)
			ans[i]-=map[i][j]*ans[j];
	}
	for(int i=1;i<=n;i++)
		printf("%.3lf ",ans[i]);
}
int main() 
{
	scanf("%d",&n);
	for(int i=1;i<=n+1;i++)
		for(int j=1;j<=n;j++)
			scanf("%lf",&a[i][j]);
	for(int i=1;i<=n;i++)
	{
		for(int j=1;j<=n;j++)
		{
			map[i][j]=2.0*(a[i][j]-a[i+1][j]);
			map[i][n+1]+=(a[i][j]*a[i][j]-a[i+1][j]*a[i+1][j]);
		}
	}
	gauss();
    return 0;
}

例 4

一个无向图,每次操作一个节点可以把它和相邻节点都反色,最少多少次能把所有节点都反色。

Luogu2962 Lights G

解析:

这道题有两种解法。

经典的异或方程组,比这道题还要简单一些。

POJ1830 开关问题

我们设 ai,ja _{i,j} 表示按下第 jj 个开关是否会影响到 ii 的状态,xix _{i} 代表是否按下第 ii 个开关,特别地,令 ai,i=1a _{i,i}=1。如果一个开关的起始状态 stist _{i},结束状态为 edied _{i}。那么可以列出如下方程组:

{a1,1x1 xor a1,2x2 xor ⋯ xor a1,nxn=st1 xor ed1a2,1x1 xor a2,2x2 xor ⋯ xor a2,nxn=st2 xor ed2⋮an,1x1 xor an,2x2 xor ⋯ xor an,nxn=stn xor edn \begin{cases} a _{1,1}x _{1}~\mathrm{xor}~a _{1,2}x _{2}~\mathrm{xor}~\cdots~\mathrm{xor}~a _{1,n}x _{n}=st _{1}~\mathrm{xor}~ed _{1}\\ a _{2,1}x _{1}~\mathrm{xor}~a _{2,2}x _{2}~\mathrm{xor}~\cdots~\mathrm{xor}~a _{2,n}x _{n}=st _{2}~\mathrm{xor}~ed _{2}\\ \vdots\\ a _{n,1}x _{1}~\mathrm{xor}~a _{n,2}x _{2}~\mathrm{xor}~\cdots~\mathrm{xor}~a _{n,n}x _{n}=st _{n}~\mathrm{xor}~ed _{n} \end{cases}

最后对这个异或方程组进行高斯消元即可,只需要把高斯消元中的运算改成异或即可。开关问题求得是不同得方案数,只需要求出自由元个数 xx,每个自由元都可以为 00 或 11,所以答案就是 2x2 ^{x}。

对于这道题的处理方法类似,但与开关问题不同的是这题求得是最少的操作次数,我们可以直接爆搜枚举每个自由元的情况,对所有不同的情况取最小值即可。

#include<cstdio>
#include<cstring>
#include<cstdlib>
#include<iostream>
#include<algorithm>
#include<cmath>
#define INF 2e9
using namespace std;
const long long maxn=41;
long long n,m,a[maxn][maxn];
long long b[maxn],ans=INF,now;
void gauss()
{
	for(long long i=1;i<=n;i++)
	{
		long long maxx=0;
		for(long long j=i+1;j<=n;j++)
			if(a[j][i]==1)
			{
				maxx=j;
				break;
			}
		if(a[maxx][i])
			swap(a[maxx],a[i]);
		for(long long k=1;k<=n;k++)
			if(k!=i&&a[k][i]==1)
				for(long long j=1;j<=n+1;j++)
					a[k][j]^=a[i][j];
	}
}
void dfs(long long x)
{
	if(now>ans)
		return ;
	if(!x)
	{
		ans=min(ans,now);
		return ;
	}
	if(a[x][x])
	{
		long long num=a[x][n+1];
		for(long long i=x+1;i<=n;i++)
			if(a[x][i])
				num^=b[i];
		b[x]=num;
		if(num==1)
			now++;
		dfs(x-1);
		if(num==1)
			now--;
	}
	else
	{
		dfs(x-1);
		b[x]=1;
		now++;
		dfs(x-1);
		b[x]=0;
		now--;
	}
}
int main()
{
	scanf("%lld%lld",&n,&m);
	for(long long i=1;i<=n;i++)
		a[i][i]=a[i][n+1]=a[n+1][i]=1;
	for(long long i=1;i<=m;i++)
	{
		long long x,y;
		scanf("%lld%lld",&x,&y);
		a[x][y]=a[y][x]=1;
	}
	gauss();
	dfs(n);
	printf("%lld",ans);
	return 0;
}

顾名思义,就是我们先把所有的灯的状态先状压一下,然后先找出前半的状态存储在 map 中,然后再搜索后半的状态,每搜出一种方案就把与它互补的前半段的方案来更新答案。这个算法也被称作 meet in middle。

/*
 * @Author: clorf 
 * @Date: 2020-09-05 22:08:59 
 * @Last Modified by: clorf
 * @Last Modified time: 2020-09-05 22:38:16
 */
// 数组一定要算好开!!!!边和点的关系
#include<cstdio>
#include<cstring>
#include<cstdlib>
#include<cmath>
#include<algorithm>
#include<ctime>
#include<map>
#define INF 1e9
using namespace std;
const int maxn=41;
const double Pi=acos(-1.0);
template<class T>void read(T &x)
{
    x=0;int f=0;char ch=getchar();
    while(ch<'0'||ch>'9') {f|=(ch=='-');ch=getchar();}
    while(ch>='0'&&ch<='9'){x=(x<<1)+(x<<3)+(ch^48);ch=getchar();}
    x=f?-x:x;
    return;
}
int n,m;
long long b[42],lim,last;
long long p[maxn],minn;
long long stepa,stepb,flag;
map<long long,int> h;
void dfs(int x,long long now,long long step)//第几个灯,目前状态,操作次数
{
    if(x==last+1)
    {
        if(now==lim)
            minn=min(minn,step);
        if(!flag)
        {
            if(!h[now]||h[now]>step)
                h[now]=step;
        }
        else
        {
            int t=h[lim-now];//互补状态的步数
            if(!t)
                return ;
            minn=min(minn,t+step);
        }
        return ;
    }
    dfs(x+1,now,step);
    dfs(x+1,now^p[x],step+1);
    return ;
}
int main()
{
    minn=INF;
    scanf("%d%d",&n,&m);
    b[1]=1;
    for(int i=2;i<=41;i++)
        b[i]=b[i-1]<<1;
    lim=b[n+1]-1;//目标状态
    for(int i=1;i<=m;i++)
    {
        int x,y;
        scanf("%d%d",&x,&y);
        p[x]+=b[y];//状压每个灯所影响的状态
        p[y]+=b[x];
    }
    for(int i=1;i<=n;i++)
        p[i]+=b[i];
    last=n/2;
    dfs(1,0,0);
    flag=1;
    last=n;
    dfs(n/2+1,0,0);//折半搜索
    printf("%lld",minn);
    return 0;
}

例 5

给定 nn 个整数,从中选出一些,它们的乘积为完全平方数,求方案数。(素因子都小于 500500)

UVA11542 Square

解析:首先要了解完全平方数的性质。将完全平方数质因数分解后,每个质因子的指数都是偶数。由于我们要从中选出一些数,那么我们就可以把未知数 xix _{i} 设成第 ii 个数选不选,系数 ai,ja _{i,j} 就直接设成第 ii 个数第 jj 个质因子的个数,发现 500500 以内的素因子只有 9595 个。因此我们需要求使得 ∀ j∈[1,95],∑i=1nai,j≡0(mod2)\forall~j\in[1,95],\sum _{i=1} ^{n} a _{i,j}\equiv 0\pmod{2} 的方案数。发现未知数只有 00 和 11 两种可能,因此我们需要对 ai,ja _{i,j} 做一些修改让其变成一个异或方程组。我们设 ai,ja _{i,j} 代表第 ii 个数第 jj 个质因子的个数是否为奇数,是就为 11,否则就为 00,最后由于要求是偶数,所以每行方程的常数项是 00。最后解出异或方程组的自由元个数 xx,由于要去掉空集的方案,答案即为 2x−12 ^{x}-1。

/*
 * @Author: clorf 
 * @Date: 2020-09-07 19:28:36 
 * @Last Modified by: clorf
 * @Last Modified time: 2020-09-07 21:24:48
 */
#include<cstdio>
#include<cstring>
#include<cstdlib>
#include<cmath>
#include<algorithm>
#include<ctime>
#define INF 1e9
using namespace std;
const int maxn=510;
const double Pi=acos(-1.0);
// 数组一定要算好开!!!!边和点的关系
// 记得看开不开long long!!!
template<class T>void read(T &x)
{
    x=0;int f=0;char ch=getchar();
    while(ch<'0'||ch>'9') {f|=(ch=='-');ch=getchar();}
    while(ch>='0'&&ch<='9'){x=(x<<1)+(x<<3)+(ch^48);ch=getchar();}
    x=f?-x:x;
    return;
}
// n和m不要搞混!!!任何循环时候都要注意!!
int t,n,m=95;
long long s[maxn];
long long a[maxn][maxn];
long long p[]={0,2,3,5,7,11,13,17,19,23,29,31,37,41,43,47,53,59,\
         61,67,71,73,79,83,89,97,101,103,107,109,113,127,\
         131,137,139,149,151,157,163,167,173,179,181,191,\
         193,197,199,211,223,227,229,233,239,241,251,257,\
         263,269,271,277,281,283,293,307,311,313,317,331,\
         337,347,349,353,359,367,373,379,383,389,397,401,\
         409,419,421,431,433,439,443,449,457,461,463,467,\
         479,487,491,499};
inline long long gauss()
{
    int i=1,j=1;
    while(i<=m&&j<=n)
    {
        int maxx=i;
        for(int k=i;k<=m;k++)
            if(a[k][j])
            {
                maxx=k;
                break;
            }
        if(a[maxx][j])
        {
            swap(a[i],a[maxx]);
            for(int k=i+1;k<=m;k++)
                if(a[k][j])
                    for(int l=j;l<=n+1;l++)
                        a[k][l]^=a[i][l];
            i++;
        }
        j++;
    }
    for(int k=i;k<=m;k++)
		if(a[k][n+1])
			return -1;
    return n-i+1;
}
int main()
{
    scanf("%d",&t);
    while(t--)
    {
        int maxx=0;
        memset(s,0,sizeof(s));
        memset(a,0,sizeof(a));
        scanf("%d",&n);
        for(int i=1;i<=n;i++)
        {
            scanf("%lld",&s[i]);
            for(int j=1;j<=95;j++)
                while(s[i]%p[j]==0)
                {
                    maxx=max(maxx,j);
                    s[i]/=p[j];
                    a[j][i]^=1;
                }
        }
        m=maxx;
        for(int i=1;i<=m;i++)
            a[i][n+1]=0;
        long long num=gauss();
        printf("%lld\n",(1ll<<(num))-1);
    }
    return 0;
}

线性基

线性基(Linear Basis\text{Linear Basis})是针对于异或空间的一个以数构成的集合,每个由数构成的序列或集合都会有至少一个线性基。要了解线性基,首先要了解基的定义等向量知识。在异或空间中,两个数(向量也行) x,yx,y 的线性组合为 ax xor byax~\mathrm{xor}~by。


线性表出

若 aa 是 AA 中元素的一个线性组合,那么也可以称作 aa 能被 AA 中的元素线性表出。

线性相关

集合内的元素线性相关,即存在一个元素是其他元素的线性组合,也等价于存在一个元素能被其他元素线性表出。

线性无关

集合内的元素线性无关,即每个元素都不是其他元素的线性组合,也等价于所有元素都不能被其他元素线性表出。

线性空间

由至少两个元素或多个元素所构成的所有线性组合的集合空间。 x,yx,y 所构成的线性空间即为所有 ax xor byax~\mathrm{xor}~by 构成的集合空间。

基(基底)

一个线性空间的基即为可以张成这个线性空间的线性无关的元素集合。

张成

由基中的元素的所有线性组合构成一个线性空间,这样的动作过程称为张成,即基张成线性空间。


一个序列的线性基即为这个序列的最小的线性无关集合,使它能张成的线性空间与所给序列的线性空间相同。这是线性基最重要的性质,通俗地讲,即线性基的元素能相互异或得到原集合的元素的所有相互异或得到的值,且线性基是满足这个性质最小的集合。

总结以下,线性基有以下性质:

  1. 线性基的元素能相互异或得到原集合的元素的所有相互异或得到的值;
  2. 线性基是满足 11 性质最小的集合,且不同的线性基的元素个数都相同;
  3. 线性基没有异或和为 00 的子集;
  4. 线性基中不同的异或组合异或出来的数都不一样;
  5. 线性基中每个元素的二进制最高位互不相同。

构造

其中 pp 数组是序列的线性基。由线性基的第三个性质可以得到,我们用 pip _{i} 来记录最高位为第 i+1i+1 位的线性基元素。 我们把原序列的数依次做 insert 操作即可。

inline void insert(long long x)
{
    for(int i=55;i>=0;i--)
    {
        if(!(x&(1ll<<i)))
            continue;
        if(!p[i])
        {
            p[i]=x;
            return ;
        }
        x^=p[i];
    }
}

那为什么这样是满足线性基性质的呢?我们尝试证明一下。

对于性质 11,我们尝试证明线性基的元素可以相互异或得到原序列的每一个元素,这样是显然满足性质 11 的,我们从插入的情况来考虑。如果某个原序列的元素 xx 不能插入线性基,那么说明插入它时异或若干个数后变成了 00,即 x xor pi xor ⋯ xor pk=0x~\mathrm{xor}~p _{i}~\mathrm{xor} ~\cdots~\mathrm{xor}~p _{k}=0,根据异或运算的性质,可以得到 pi xor ⋯ xor pk=xp _{i}~\mathrm{xor} ~\cdots~\mathrm{xor}~p _{k}=x,说明 xx 可以由线性基中的元素异或得到;如果某个原序列的元素 xx 成功插入了线性基,那么他在插入之前可能异或了若干个数,即 x xor pi xor ⋯ xor pk=plx~\mathrm{xor}~p _{i}~\mathrm{xor} ~\cdots~\mathrm{xor}~p _{k}=p _{l},同理可得 pl xor pi xor ⋯ xor pk=xp _{l}~\mathrm{xor}~p _{i}~\mathrm{xor} ~\cdots~\mathrm{xor}~p _{k}=x,此时 xx 也可以由线性基中的元素异或得到。综上,性质 11 得证。

对于性质 22,我们还是分类讨论一下。假设原序列里面的所有元素都能插入到线性基里面,那么不管什么顺序插入,最终的元素都是 nn 个;如果原序列的一些元素不能插入进去,那么它一定由线性基里面的某些元素异或得到,例如 pi xor pj xor pk=xp _{i}~\mathrm{xor}~p _{j}~\mathrm{xor}~p _{k}=x。我们可以适当改变插入顺序来使 xx 与 pi,pj,pkp _{i},p _{j},p _{k} 交换,因为它们中的其中一个都可以由其它三个异或得到。因此改变插入顺序后,虽然线性基中的元素变了,但元素个数不会变。所以元素的数量是一定的。在此基础上我们尝试删掉线性基里面任意一个元素,这样都会使得原序列里的一些元素无法异或得到,所以线性基里的每个元素都是必要的。所以线性基的元素个数一定是在性质 11 的前提下最少的。

证明了线性基这两个最最基本的性质后,其他的性质都可以相应的推出来,此处不再做过多证明。

求最大值

求一个序列中,取若干个数,使它们的异或和最大,输出最大异或和。

Luogu3812 线性基

求出线性基后,我们从线性基的最高位开始往下搜,对于每一个线性基元素,贪心地异或它,如果异或它能使答案变得更大,那么异或它即可。

inline long long find()
{
    long long x=0;
    for(int i=55;i>=0;i--)
        if((ans^p[i])>ans)
            ans^=p[i];
    return ans;
}

为什么这么贪心是对的呢?这要用到线性基的性质 55,线性基的最高位各不相同。我们从高位到低位枚举,由于 pip _{i} 的第 i+1i+1 位为 11,我们就需要判断是否要异或这个 11。如果异或后能让答案变大,说明 ansans 第 i+1i+1 位为 00,这一位由 00 变成了 11,贡献比 1∼i1\sim i 位的贡献都要大,所以不需要管 pip _{i} 后面位对答案的影响。否则说明 ansans 这一位为 11,异或后会使答案变小,我们不异或它即可。

求次大值

求一个序列的严格次大异或和。

BZOJ4269 再见Xor

求出最大值后,再找到最大值最低为 11 的那一位 ii,且那一位存在线性基元素,用最大值异或 pip _{i} 即可。

inline void find()
{
    ans1=ans2=0;
    for(int i=62;i>=0;i--)
        if((ans1^p[i])>ans1)
            ans1^=p[i];
    for(int i=0;i<=62;i++)
        if((ans1&(1<<i))&&p[i])
        {
            ans2=ans1^p[i];
            break;
        }
}

由于要严格次大值,我们肯定想把最大值的末尾为 11 的那一位异或掉就是次大值了。但是如果那一位不存在线性基元素,即 pi=0p _{i}=0,那么此时次大值等于最大值,并不满足严格次大值的定义。因此还要保证那一位存在线性基元素。

求最小值

分为两个子问题。

求一个序列的线性基内的元素能异或出的最小值和求一个序列内元素能异或出来的最小值。

子问题 1

求一个序列的线性基内的元素能异或出的最小值。

显然就是线性基中最小的元素,因为它无论异或谁都会变大。

子问题 2

求一个序列内元素能异或出来的最小值。

这个时候就要分类讨论一下了。如果序列有元素未插入线性基,说明线性基中的元素能异或得到它,再异或它本身就为 00,此时最小值为 00;如果没有,那么答案依然是线性基中最小的元素。

求 k 小值

求一个序列内元素能异或出来的第 kk 小值。

我们要先把求出的线性基进行一番处理,保证在所有的线性基元素中,只有 pip _{i} 的第 i+1i+1 位为 11,其他的第 i+1i+1 位都为 00。即对于每一个 pip _{i},对于 ∀ j∈[0,i]\forall~ j\in[0,i],如果 pip _{i} 的第 j+1j+1 位为 11,那么就让 pip _{i} 异或 pjp _{j}。这样做的目的是让第 i+1i+1 位的贡献只由 pip _{i} 决定。最后把 kk 转成二进制后,如果 kk 的第 ii 位为 11,那么答案就异或线性基的从小到大的第 ii 个元素(不是直接异或 pi−1p _{i-1},因为可能第 ii 位不存在线性基的元素)。

inline void change()
{
    for(int i=1;i<=55;i++)
        for(int j=0;j<i;j++)
            if(p[i]&(1ll<<j))
                p[i]^=p[j];
}
inline long long kth(int k)
{
    if(k==1&&cnt<n)//cnt指原序列插入线性基中的个数,如果有没插入进去的,那么最小值为0
        return 0;
    if(cnt<n)//把0的情况排除
        k--;
    change();
    long long ans=0;
    for(int i=0;i<=55;i++)
        if(p[i])
        {
            if(k&1)
                ans^=p[i];
            k>>=1;
        }
}

为什么是正确的呢?可以发现我们处理线性基过后,第 ii 位的贡献只能由 pi−1p _{i-1} 决定(如果存在 pi−1p _{i-1} 的话)。我们设线性基元素的选择方案的状态为 qq,即如果 qq 的第 ii 位为 11 那么就异或从小到大的第 ii 个线性基元素,发现 qq 增大了,那么异或出来的结果也会增大,我们可以分类讨论来证明。qq 增大 11 只有两种方式,一是 qq 进位了,某一位 ii 由 00 变成了 11,比这位低的所有位都变成了 00,例如 qq 从 1111 变成了 100100,从只异或线性基的第 11 个和第 22 个元素变成只异或线性基的第 33 个元素,此时答案变大了,第 ii 位由原来的 00 变成了 11,它的贡献已经比第 1∼i−11\sim i-1 位的贡献大了,所以即使第 1∼i−11\sim i-1 变成 00 答案依旧变大;第二种方式是 qq 的某位非最高位由 00 变成了 11,其他不变,即 qq 没进位,例如 1010 变成 1111,此时虽然答案的最高位没变,但同时它的相应那位的低位也从 00 变成了 11,答案依然增大了。因此 qq 和答案成一个一一对应的关系,qq 增大 11,答案也从第 kk 小变成第 k+1k+1 小。因此若选择方案为 qq,那么答案就是第 qq 小值。感性一点理解,其实就相当于处理后的线性基的作用就是只提供自己最高位的那个 11,提供出的 11 和 kk 的每一位上的 11 对应。

判断一个数是否可构造出来

把它尝试插入进线性基里面去,如果成功说明不能,否则可以得到。

删除

BZOJ4184 shallot

线性基的删除操作。我们采取离线做法,找到线性基中第 ii 个数,可以代替它的那些数的集合 did _{i}。然后使线性基优先存删除时间晚的即可,这样是为了避免删掉后其他序列中的数可以填进来代替它的问题,最后在删除时直接在线性基里面清零即可。其余的插入和求最大值操作一样。具体实现就是在插入一个新数时,比较一下新数的删除时间和当前位线性基元素的删除时间,如果比它晚就换掉,否则异或他继续枚举。代码看这里大佬的代码。

另一种方法是用线段树时间分治维护线性基,可以自己上网查到。

合并

把一个线性基的所有元素插入另一个线性基中即可。

区间询问

多次询问,每次询问一个区间 [l,r][l,r] 内的数的最大异或和。

CF1100F Ivan and Burgers

维护一个前缀线性基(因为异或运算满足可加性)。多开一个数组 posipos _{i} 代表线性基第 i+1i+1 个元素在原序列中对应的数在原序列中的编号。 如果我们求出了 [1,r][1,r] 的线性基,那么其中所有 pos≥lpos\ge l 的数构成的线性基就是 [l,r][l,r] 的线性基。发现在 rr 一定的情况下,pospos 值越大贡献越大(这个数所属的区间越多)。因此我们要保证 posipos _{i} 最大,插入时类比删除时比较删除时间一样依次比较即可。

区间修改

区间修改:对一个区间内的元素全异或同一个数。

CF587E Duff as a Queen

对于异或我们可以类比简单的加法运算(因为都具有可加性),加法运算的区间修改可以通过差分数组解决。这道题我们也用差分数组,设 bi=ai xor ai−1b _{i}=a _{i}~\mathrm{xor}~a _{i-1},b1=a1b _{1}=a _{1},这样每次只用修改 bib _{i} 的区间端点的两个值,同时维护它的前缀异或即可,对于询问操作 [l,r][l,r],答案即为 al,bl+1∼bra _{l},b _{l+1}\sim b _{r} 所组成的线性基。因此,bb 需要单点修改,区间查询,用线段树即可,注意线段树的每个元素是一个线性基,且支持合并操作。aa 需要区间修改,单点查询,又由于 bb 已经是 aa 的差分数组,直接用树状数组即可。

代码来自大佬的代码。

#include<bits/stdc++.h>
#define maxn 200005
using namespace std;
struct Base{
	int a[31],cnt;
	void clear(){memset(a,0,sizeof a),cnt=0;}
	void ins(int x){
		if(cnt==30) return;
		for(int i=29;i>=0&&x;i--) if(x>>i&1){
			if(a[i]) x^=a[i];
			else {a[i]=x,cnt++;return;}
		}
	}
	void merge(const Base &B){
		for(int i=29;i>=0;i--) if(B.a[i]) ins(B.a[i]);//线性基合并
	}
}t[maxn<<2];//线段树,维护b的线性基
int n,m,a[maxn],b[maxn];
int arr[maxn];//树状数组,维护a
void upd(int i,int v){for(;i<=n;i+=i&-i) arr[i]^=v;}
int qxor(int i){int s=0;for(;i;i-=i&-i) s^=arr[i];return s;}
void upd(int i){t[i]=t[i<<1],t[i].merge(t[i<<1|1]);}//线段树父亲的线性基是两个孩子线性基合并起来
void build(int i,int l,int r){
	if(l==r) return t[i].ins(b[l]);
	int mid=(l+r)>>1;
	build(i<<1,l,mid),build(i<<1|1,mid+1,r);
	upd(i);
}
void mdf(int i,int l,int r,int x){
	if(l==r) {t[i].clear(),t[i].ins(b[x]);return;}
	int mid=(l+r)>>1;
	x<=mid?mdf(i<<1,l,mid,x):mdf(i<<1|1,mid+1,r,x);
	upd(i);
}
void qry(int i,int l,int r,int x,int y){
	if(x<=l&&r<=y) return t[0].merge(t[i]);//将t[i]的线性基合并到用来储存答案的t[0]中
	int mid=(l+r)>>1;
	if(x<=mid) qry(i<<1,l,mid,x,y);
	if(y>mid) qry(i<<1|1,mid+1,r,x,y);
}
int main()
{
	scanf("%d%d",&n,&m);
	for(int i=1;i<=n;i++) scanf("%d",&a[i]),b[i]=a[i]^a[i-1],upd(i,b[i]);
	build(1,1,n);
	for(int op,l,r,x;m--;){
		scanf("%d%d%d",&op,&l,&r);
		if(op==1){
			scanf("%d",&x);
			upd(l,x),upd(r+1,x);//树状数组维护差分数组,区间修改
            b[l]^=x,b[r+1]^=x;//更新数组的值
			mdf(1,1,n,l); if(r<n) mdf(1,1,n,r+1);//将更新后的值重新更新到线段树中,单点修改
		}
		else{
			t[0].clear(),t[0].ins(qxor(l));//qxor,树状数组维护差分数组,单点查询a _{l},t[0]不属于线段树,这里只用来储存答案
            if(l<r) qry(1,1,n,l+1,r);//线段树区间查询
			printf("%d\n",1<<t[0].cnt);
		}
	}
}