有 \(n\) 个形如 \(a_1x_1+a_2x_2+\cdots+a_nx_n=b\) 的方程。解方程组。
有唯一解则将其求出，否则输出No Solution。
把方程组用矩阵形式写出来
（这个东西叫增广矩阵）
考虑选 \(x_1\) 为主元。
若所有方程中 \(x_1\) 的系数均为 0，则 \(x_1\) 是个自由未知量，方程组无唯一解。直接结束即可。
否则，选出 \(|a_{i,1}|\) 最大的第 \(i\) 行（这样精度误差最小，虽然我也不知道为什么）
然后用第 \(i\) 个方程和其它所有方程做一遍加减消元，消去 \(x_1\)。
不妨设 \(i=1\)，则得到新的增广矩阵大概是下面这样
0 表示该项已被消去，1 表示其它。
对 \(x_2\) 执行相同操作。注意之前选过的行不能再选。
设这次选了第 \(j(j\ne i)\) 行，以 \(j=2\) 为例，新的增广矩阵是
重复上述过程，最终得到
到这里就算是解完了。
#include<stdio.h>
#include<string.h>
#define mcpy(a,b) memcpy(a,b,sizeof ta)
const int N=110; int n; double ta[N],a[N][N];
inline double abs(double x) { return x<0?-x:x; }
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]); // a[i][n+1]==b[i]
for (int i=1,t=1; i<=n; t=++i) {
for (int j=i+1; j<=n; ++j)
abs(a[j][i])>abs(a[t][i])&&(t=j);
mcpy(ta,a[i]),mcpy(a[i],a[t]),mcpy(a[t],ta);
// abs(a[t][i]) 最大，将第 t 行与第 i 行交换，方便操作
if (!a[i][i]) return puts("No Solution"),0; // 无唯一解
for (int j=1; j<=n; ++j) if (j!=i) { // 在其它行中消去 xi
double t2=a[j][i]/a[i][i];
for (int k=i+1; k<=n+1; ++k)
a[j][k]-=a[i][k]*t2;
}
}
for (int i=1; i<=n; ++i)
printf("%.2lf\n",a[i][n+1]/a[i][i]);
return 0;
}
时间复杂度 \(O(n^3)\)。
另一道模板题：洛谷P2455 [SDOI2006]线性方程组
这次需要区分无解和有无数组解。
其实也很简单
正常做消元，如果遇到某一列系数全为 0，就直接摆烂（反正也消不出啥东西），直接跳到下一列
这样得到的增广矩阵依然是
但这次写着 1 的位置可能会出现 0。
首先检查一下有没有 \(0x\ne 0\) 这样的方程出现。有则说明无解。
否则再检查有没有 \(0x=0\) 这样的方程出现。有则说明有无数组解。
否则说明有唯一解。
然后我们获得了 80pts 的好成绩。
看一眼报错信息，第三个点无穷解判成无解了，第八个点好像把 0.00 输出成了 -0.00。
后者是诡异的 feature，输出时加上一个很小的数即可解决。
前者更加离谱，我 debug 了一晚上也没找到错。
三天后我终于在题解区翻到一组 hack 数据：
2
0 2 3
0 0 0
答案为 0，而程序会输出 -1
因为第一次消元时选了第一行，然后第二行没法消了，但原本是可以消的。所以寄了。
原因是消元时只是贪心地选择了当前列系数最大的一行。没有考虑到对后面造成的影响。
补救措施：微调一下选择行的方案，如果两行的系数绝对值相等，选后面的数更小的。
#include<stdio.h>
#include<string.h>
#define mcpy(a,b) memcpy(a,b,sizeof ta)
#define rep(i,l,r) for (int i=l; i<=r; ++i)
const int N=52; const double eps=1e-6;
int n,f1,f2; double ta[N],a[N][N];
inline double abs(double x) { return x<0?-x:x; }
inline int F(double x) { return (int)(x*100); }
bool cmp(int i,int j,int k,bool f) { // 新的选择策略
double t1=abs(a[i][k]),t2=abs(a[j][k]);
if (t1!=t2) return f^(t1>t2);
// 首先保证当前列上的系数绝对值最大，其次保证后面的列上的系数绝对值最小
return k==n?0:cmp(i,j,k+1,1);
}
int main() {
scanf("%d",&n);
rep(i,1,n) rep(j,1,n+1) scanf("%lf",&a[i][j]);
for (int i=1,t=1; i<=n; t=++i) {
rep(j,i+1,n) cmp(j,t,i,0)&&(t=j);
mcpy(ta,a[i]),mcpy(a[i],a[t]),mcpy(a[t],ta);
if (!F(a[i][i])) continue;
rep(j,1,n) if (j!=i) {
double t=a[j][i]/a[i][i];
rep(k,i+1,n+1) a[j][k]-=t*a[i][k];
}
}
rep(i,1,n) F(a[i][i])||(f1=1,F(a[i][n+1])&&(f2=1));
if (f1) return puts(f2?"-1":"0"),0;
rep(i,1,n) printf("x%d=%.2lf\n",i,a[i][n+1]/a[i][i]+eps);
return 0;
}
做完之后还是感觉这代码丑的一批。
不过反正能过，复杂度也没假，就这样吧（
补个例题
给出一个 \(n\) 维球面上 \(n+1\) 个点的坐标，确定球心位置。
设球心位置为 \((x_{0,1},x_{0,2},\cdots,x_{0,n})\)，半径为 \(r\)。
对于第 \(i\) 个点，我们有方程：\(\sum\limits_{j=1}^n(x_{0,j}-x_{i,j})^2=r^2\)。
令 \(k=\left(\sum\limits_{j=1}^n{x_{0,j}}^2\right)-r^2,S_i=\sum\limits_{j=1}^nx_{i,j}^2\)，
则 \(2x_{i,1}\cdot x_{0,1}+2x_{i,2}\cdot x_{0,2}+\cdots+2x_{i,n}\cdot x_{0,n}-k=S_i\)
共 \(n+1\) 个方程，\(n+1\) 个未知数。高斯消元解方程组即可。