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

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

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

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

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

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

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

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

• 系数矩阵

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

• 增广矩阵

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

(Example:)

假如现有一方程组：

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

其系数矩阵即为：

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

(？？？)

增广矩阵为：

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

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

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

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

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

• 交换两个方程

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

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

---

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

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

增广矩阵：

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

---

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

((1)*=-9:)

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

增广矩阵：

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

---

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

((1)+=(2):)

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

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

增广矩阵：

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

---

• 交换两个方程

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

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

增广矩阵：

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

---

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

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

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

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

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

若过程中无法消元（找不到一个方程，(x_1sim 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
",a[i][n+1]/a[i][i]);
	return 0;
}

---

P4035 [JSOI2008]球形空间产生器

稍微转换一下题目

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

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

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

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

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

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

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

[sum_{1le jle n}2(a_{i,j}-a_{i+1,j})x_j=sum_{1le jle 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?'
':' ');
    return 0;
}

---

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

参考资料：

各省省选题目

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