[数论]高斯消元学习笔记
\(Q:\)高斯消元是什么？听起来好高级啊？？
\(A:\)二元一次方程组解过吗？那就是高斯消元。
- 高斯消元（\(Gauss\)），是一种用来求解线性方程组的方法，在\(OI\)竞赛中广泛使用。
首先对高斯消元做一些准备：
\(Q:\)什么是线性方程组？
\(A:\)鸡兔同笼方程组 由\(M\)个\(N\)元一次方程所构成的方程组。
为了简化表达，做出如下定义：
-
系数矩阵
把线性方程组未知数的系数写成一个\(M\times N\)的矩阵，即为“系数矩阵”。
-
增广矩阵
把线性方程组等号右边的常数写成一个\(M\times 1\)的矩阵，拼在系数矩阵右边，即为一个大小为\(M\times (N+1)\)的“增广矩阵”。
\(Example:\)
假如现有一方程组：
其系数矩阵即为：
(？？？)
增广矩阵为：
\(Q:\)这些很简单啊，接下来要干什么呢？
\(A:\)接着我们对方程组进行求解，包括如下\(3\)个步骤：
-
对其中一个方程乘以一个非\(0\)常数
-
用一个方程加上另一个方程
-
交换两个方程
称以上操作为矩阵的“初等行变换”。
\(Q:\)这些毫无意义的语言 都是什么意思啊？
以一个二元一次方程组为例：
增广矩阵：
- 对其中一个方程乘以一个非\(0\)常数:
\((1)*=-9:\)
增广矩阵：
- 用一个方程加上另一个方程
\((1)+=(2):\)
增广矩阵：
- 交换两个方程
\(Swap((1),(2)):\)
增广矩阵：
怎么样，是不是很简单啊？
初等行变换后的矩阵被称为“阶梯型矩阵”（是不是很像？并没有）。
系数部分称为“上三角矩阵”。
最后从下到上更新一遍即可求出方程组的解。
另外，若在求解过程中出现了\(0=a\)的情况，则方程组无解。
若过程中无法消元（找不到一个方程，\(x_1\sim x_{i-1}\)系数为\(0\)，\(x_i\)系数不为\(0\)），则解不唯一，称此类\(x_i\)为自由元。
时间复杂度 \(O(n^3)\)
空间复杂度 \(O(n^2)\)
接下来是实际应用。
这就是裸的模板题了。
时间复杂度 \(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;
}
稍微转换一下题目
设已知点为\(a_i\)，求一个点\((x_1,x_2,\dots,x_n)\)，使得存在常数\(Dis\)，满足
这样就有了\(n+1\)个方程，但发现无法求解。
考虑构造\(n\)个\(n\)元一次方程组进行高斯消元。
用相邻的两个方程做差，得：
高斯消元即可。
时间复杂度 \(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;
}
难题就不继续探讨了，毕竟我也不会
参考资料：
各省省选题目
《算法竞赛进阶指南》 -- 李煜东
