[数论]高斯消元学习笔记

\(Q:\)高斯消元是什么？听起来好高级啊？？

\(A:\)二元一次方程组解过吗？那就是高斯消元。

• 高斯消元（\(Gauss\)），是一种用来求解线性方程组的方法，在\(OI\)竞赛中广泛使用。

首先对高斯消元做一些准备：

\(Q:\)什么是线性方程组？

\(A:\)鸡兔同笼方程组 由\(M\)个\(N\)元一次方程所构成的方程组。

为了简化表达，做出如下定义：

• 系数矩阵

  把线性方程组未知数的系数写成一个\(M\times N\)的矩阵，即为“系数矩阵”。

• 增广矩阵

  把线性方程组等号右边的常数写成一个\(M\times 1\)的矩阵，拼在系数矩阵右边，即为一个大小为\(M\times (N+1)\)的“增广矩阵”。

\(Example:\)

假如现有一方程组：

\[\begin{cases}x_1+x_2+4x_3=17\\5x_1+x_2+4x_3=25\\8x_1+x_2+0x_3=19\end{cases} \]

其系数矩阵即为：

\[\begin{bmatrix}1&1&4\\5&1&4\\8&1&0\end{bmatrix} \]

(？？？)

增广矩阵为：

\[\begin{bmatrix}1&1&4&17\\5&1&4&25\\8&1&0&19\end{bmatrix} \]

\(Q:\)这些很简单啊，接下来要干什么呢？

\(A:\)接着我们对方程组进行求解，包括如下\(3\)个步骤：

• 对其中一个方程乘以一个非\(0\)常数

• 用一个方程加上另一个方程

• 交换两个方程

称以上操作为矩阵的“初等行变换”。

\(Q:\)这些毫无意义的语言 都是什么意思啊？

---

以一个二元一次方程组为例：

\[\begin{cases}x_1+9x_2=21\quad(1)\\9x_1+x_2=29\quad(2)\end{cases} \]

增广矩阵：

\[\begin{bmatrix}1&9&21\\9&1&29\end{bmatrix} \]

---

• 对其中一个方程乘以一个非\(0\)常数:

\((1)*=-9:\)

\[\begin{cases}-9x_1+-81x_2&=-189\quad(1)\\9x_1+x_2&=29\quad(2)\end{cases} \]

增广矩阵：

\[\begin{bmatrix}-9&-81&-189\\9&1&29\end{bmatrix} \]

---

• 用一个方程加上另一个方程

\((1)+=(2):\)

\[\begin{cases}-80x_2&=-160\\9x_1+x_2&=29\end{cases} \]

\[\begin{cases}x_2&=2\\9x_1+x_2&=29\end{cases} \]

增广矩阵：

\[\begin{bmatrix}0&1&2\\9&1&29\end{bmatrix} \]

---

• 交换两个方程

\(Swap((1),(2)):\)

\[\begin{cases}9x_1+x_2&=29\\x_2&=2\end{cases} \]

增广矩阵：

\[\begin{bmatrix}9&1&29\\0&1&2\end{bmatrix} \]

---

怎么样，是不是很简单啊？

初等行变换后的矩阵被称为“阶梯型矩阵”（是不是很像？并没有）。

系数部分称为“上三角矩阵”。

最后从下到上更新一遍即可求出方程组的解。

另外，若在求解过程中出现了\(0=a\)的情况，则方程组无解。

若过程中无法消元（找不到一个方程，\(x_1\sim x_{i-1}\)系数为\(0\)，\(x_i\)系数不为\(0\)），则解不唯一，称此类\(x_i\)为自由元。

时间复杂度 \(O(n^3)\)

空间复杂度 \(O(n^2)\)

接下来是实际应用。

---

P3389 【模板】高斯消元法

这就是裸的模板题了。

时间复杂度 \(O(n^3)\)

空间复杂度 \(O(n^2)\)

代码：

#include <cstdio>
#include <algorithm>

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

int n;
double a[105][105];
const double Eps=1e-7;//精度视情况而定，不要太小或太大

int main()
{
	scanf("%d",&n);
	for(int i=1;i<=n;++i)
		for(int j=1;j<=n+1;++j)
			scanf("%lf",&a[i][j]);//共同构成增广矩阵
	for(int i=1;i<=n;++i)//进行n次消元，消去第i项未知数
	{
		for(int j=1;j<=n;++j)//找一个x[i]系数不为0的方程进行消元
			if(Abs(a[j][i])>Eps)
				for(int k=1;k<=n+1;++k)
					std::swap(a[i][k],a[j][k]);
		if(Abs(a[i][i])<=Eps)return puts("No Solution"),0;//出现自由元，无解
		//消去其他方程x[i]系数
		for(int j=1;j<=n;++j)
			if(i!=j)
			{
				double Rate=a[j][i]/a[i][i];//相对比例
				for(int k=1;k<=n+1;++k)a[j][k]-=a[i][k]*Rate;//消元
			}
	}
	for(int i=1;i<=n;++i)printf("%.2f\n",a[i][n+1]/a[i][i]);
	return 0;
}

---

P4035 [JSOI2008]球形空间产生器

稍微转换一下题目

设已知点为\(a_i\)，求一个点\((x_1,x_2,\dots,x_n)\)，使得存在常数\(Dis\)，满足

\[\forall i\in[1,n+1],\sum_{1\le j\le n}(a_{i,j}-x_j)^2=Dis \]

这样就有了\(n+1\)个方程，但发现无法求解。

考虑构造\(n\)个\(n\)元一次方程组进行高斯消元。

用相邻的两个方程做差，得：

\[\sum_{1\le j\le n}(a_{i,j}-x_j)^2-\sum_{1\le j\le n}(a_{i+1,j}-x_j)^2=Dis-Dis \]

\[\sum_{1\le j\le n}(a_{i,j}-x_j)^2-(a_{i+1,j}-x_j)^2=0 \]

\[\sum_{1\le j\le n}2(a_{i,j}-a_{i+1,j})x_j=\sum_{1\le j\le n}(a_{i,j}^2-a_{i+1,j}^2) \]

高斯消元即可。

时间复杂度 \(O(n^3)\)

空间复杂度 \(O(n^2)\)

代码：

#include <cmath>
#include <cstdio>
#include <algorithm>

const double Eps=1e-8;

int n;
double p[15][15],b[15],c[15][15];

int main()
{
    scanf("%d",&n);
    for(int i=1;i<=n+1;++i)
        for(int j=1;j<=n;++j)
            scanf("%lf",&p[i][j]);
    for(int i=1;i<=n;++i)
        for(int j=1;j<=n;++j)
        {
            c[i][j]=(p[i][j]-p[i+1][j])*2;
            b[i]+=p[i][j]*p[i][j]-p[i+1][j]*p[i+1][j];//这里把系数和常数分开放
        }
    for(int i=1;i<=n;++i)
    {
        for(int j=i;j<=n;++j)
            if(fabs(c[j][i])>Eps)
            {
                for(int k=1;k<=n;++k)std::swap(c[i][k],c[j][k]);
                std::swap(b[i],b[j]);
            }
        for(int j=1;j<=n;++j)
        {
            if(i==j)continue;
            double Rate=c[j][i]/c[i][i];
            for(int k=i;k<=n;++k)c[j][k]-=c[i][k]*Rate;
            b[j]-=b[i]*Rate;
        }
    }
    for(int i=1;i<=n;++i)
        printf("%.3f%c",b[i]/c[i][i],i==n?'\n':' ');
    return 0;
}

---

难题就不继续探讨了，毕竟我也不会

参考资料：

各省省选题目

《算法竞赛进阶指南》 -- 李煜东