高斯消元
Gaussian Elimination
给定 \(n\) 个形如 \(\displaystyle b_i=\sum_{i=1}^na_ix_i\) 的关于 \(x_i\) 方程，求解。
首先可以构建一个 \(n\times(n+1)\) 的矩阵
那么每一行就可以表示一个方程，同时每一列也表示同一个元。
那么我们可以考虑用第一行来消掉第一个元，第二行来消掉第二个元……到最后剩下一行一个元，就是解一元一次方程了。
然后再把消了的元拼回去。
消元的过程用前两行举个例子：用第一行消去第一个元
然后第一个和式同时乘以 \(rate=\dfrac{a_{2,i}}{a_{1,i}}\) ，然后第二个再减他，那么出来的结果 \(x_1\) 前面的系数就会变成 \(0\) 。用同样的方法处理第三行、第四行……这样剩下的 \(n-1\) 行中 \(x_1\) 就都被消去了。
大概就是这样
for(int i=0;i<n;++i) { // 用第 i 行
for(int j=i+1;j<n;++j) { // 消去第 j 行中第 i 个元
const double rate=matrix[j][i]/matrix[i][i];
for(int k=i;k<=n;++k) matrix[j][k]-=rate*matrix[i][k]; // 注意最后有一个b, 所以循环到n
}
}
然后我们期望出来的矩阵应该是这样的：每行都比上一行少一个元，最后一行一个元一个和就是一元一次。
（第一行其实没有动过所以不打撇）第一行没有零，然后后面每一行都比前面多一个零，这就是我们所期待的结果。
然而现实很残酷。
先看代码里面的 bug ，上面代码第三行 rate=matrix[j][i]/matrix[i][i]; 写了这么久代码，应该下意识反应过来：除数是 \(0\) 怎么办？？？
因为给定方程组不一定没有零，而且也不能保证消元的过程中不会出现零，所以这个担心是必要的。
那么其实我们往后面的行找一下，说不定就能找到这个元系数不是零的式子呢？
所以我们把代码修改一下
for(int i=0;i<n;++i) {
int flag=0;
if(matrix[i][i]==0) for(int j=i;j<n;++j) if(matrix[j][i]!=0) { swap(matrix[j],matrix[i]); break; }
if(matrix[i][i]==0) continue;
for(int j=i+1;j<n;++j) {
const double rate=matrix[j][i]/matrix[i][i];
for(int k=i;k<=n;++k) matrix[j][k]-=rate*matrix[i][k];
}
}
注意到 官网 在对 swap 的复杂度(Complexity) 的分析时称，如果对象是数组，那么时间复杂度在长度上线性，那么这样一来时间就多了很多了。
官网上给出了 C++11 时候对于数组对象的 swap 重载
The behavior of this function template is equivalent to:
template <class T> void swap (T& a, T& b) { T c(std::move(a)); a=std::move(b); b=std::move(c); } template <class T, size_t N> void swap (T (&a)[N], T (&b)[N]) { for (size_t i = 0; i<N; ++i) swap (a[i],b[i]); }
可以很明显看见时间复杂度和数组元素个数成线性关系。虽然在 C++98 中没有这样写重载，但是根据实例测试
int a[2][1000000];
for(int i=1;i<1000;++i) swap(a[0],a[1]);
这段代码在我的机子上跑了 \(3s\) ，这足以说明所谓的“常数”已经非常大。实际上他不是常数，而是和数组长度成线性关系。
能不能用结构体 ? 也不行，根据官网中底层代码，不是数组对象的话复杂度是“常数”。但是看清楚这句话：
Non-array: Constant: Performs exactly one construction and two assignments (although, notice that each of these operations works on its own complexity)
什么意思？大概就是说虽然非数组对象的复杂度是“常数”，但是底层进行三个操作：一个拷贝构造（没有就普通构造），两个赋值。虽然如此，注意到这三个操作各自的复杂度不能忽略！如果赋值操作是 \(\mathcal O(n)\) 的，那么这个 swap 就得老老实实是 \(\mathcal O(n)\) 了。
但是注意到上面有这样两句话：
Notice how this function involves a copy construction and two assignment operations, which may not be the most efficient way of swapping the contents of classes that store large quantities of data, since each of these operations generally operate in linear time on their size.
Large data types can provide an overloaded version of this function optimizing its performance. Notably, all standard containers specialize it in such a way that only a few internal pointers are swapped instead of their entire contents, making them operate in constant time.
需要注意的是这个功能包含了一个拷贝构造和两个赋值操作。这些操作可能并不是交换储存了大量数据的文本或者对象最有效的方式。因为这些操作一般来说都需要在长度上线性的时间。
（但是）大数据类型可以提供这个功能的重载来优化他的执行效率。值得注意的是，所有标准容器（stl）提供了一种专门的方式，只交换内部的一些指针而不是整个文本对象，这使得操作的时间复杂度是常数。
也就是说，我们可以用 stl 代替数组！最简单就是用一个 vector。
但是有没有更好的方法？
有的，我们先定义一个 整形指针 数组，然后对于每个指针再指向一个数组，那么每次交换就只是指针的交换，是彻彻底底的 \(\mathcal O(3)\)
好了扯远了
再放一下上面讲到的代码
for(int i=0;i<n;++i) {
if(matrix[i][i]==0) for(int j=i;j<n;++j) if(matrix[j][i]!=0) { swap(matrix[j],matrix[i]); break; }
if(matrix[i][i]==0) continue;
for(int j=i+1;j<n;++j) {
const double rate=matrix[j][i]/matrix[i][i];
for(int k=i;k<=n;++k) matrix[j][k]-=rate*matrix[i][k];
}
}
我们暂且不看效率问题，现在把 matrix[j][i]!=0 的那一行调了上来，那么就可以用这一行消元了。
但是，如果搜索完了后面 \(n-i\) 行都没有一个不是零呢？
那么就说明这个元是一个自由 元！他的取值和后面的元的解无关，而前面的元就因为这个元没有唯一解，那么整个方程就没有唯一解。
同时我们也可以发现，如果有一个自由元，那么整个方程组就可以变成 \(n\) 条方程， \(n-1\) 个未知数，有可能无解。
而这个无解的表现就是后面某一步发现 \(0 x=b, b\neq 0\)，这样就是无解的。
根据我们消元的步骤，这样消元下来下一行至少比上一行少一个元。如果少的多于一个元那么多出来的就是自由元。
也就是说最后出来的矩阵可以分为两部分，其中左下角部分为零。这两个部分有一个分界线。
当且仅当这条分界线是直线，也就是每一行都严格比上一行少一个元，或者书该矩阵为上三角矩阵的时候，该方程组有唯一解。
那么我们就可以从最后一行做起，先解一元一次方程。得到一个解之后倒数第二行的方程就从二元变成了一元，也可以求解了……
如此类推，一直推到第一行，解 \(n\) 个一元一次方程，就能解出所有未知数。
如果某一行比下一行多了不止一个元，那么有些元就是自由元，这些自由元直接全部给 \(0\) 就可以了。
如果过程中发现某个一元一次方程无解，那么整个方程组就无解。
int Elimination() {
/// @return : the count of the free element(s), or -1 when no solution.
int free=0;
for(int i=0;i<n;++i) {
if(fabs(matrix[i][i])<eps) for(int j=i+1;j<n;++j) if(matrix[j][i]!=0) Swap(matrix[j],matrix[i]);
if(fabs(matrix[i][i])<eps) { ++free; continue; }
for(int j=i+1;j<n;++j) {
if(j==i) continue;
const long double rate=matrix[j][i]/matrix[i][i];
for(int k=i;k<=n;++k) matrix[j][k]-=matrix[i][k]*rate;
}
}
int now=n-1;
for(int i=n-1;i>=0;--i) {
while(now and fabs(matrix[i][now-1])>eps) result[now--]=0;
double sum=0;
for(int j=now+1;j<n;++j) sum+=matrix[i][j]*result[j];
if(fabs(matrix[i][now])<eps) { if(fabs(matrix[i][n]-sum)>eps) return -1; else continue; }
result[now]=(matrix[i][n]-sum)/matrix[i][now];
now--;
}
while(now>=0) result[now--]=0;
return free;
}
模
取模意义下的高斯消元怎么做？
其实就是运算替换成取模意义下的等价运算即可。
加减乘不说，除法。如果模数是质数直接乘法逆元上去。
如果模数不是质数，考虑使用更相减损法。
比如说对于任意两行 \(a,\dots\\b,\dots\) ，现在我们要把 \(b\) 这两个系数消去。
不取模的话，首先第一行全部乘以 \(\dfrac ba\) ，然后去减第二行，就可以消去 \(b\) 了。
而在模意义下 \(a\) 不一定有逆元，这个时候就使用更相减损法。
设 \(a<b\) ，那么就第二行不停减去第一行，一直减去 \(\left\lfloor\dfrac ba\right\rfloor\) 个，这个时候可以保证 \(a>b\) ，或者说这一步操作就是 b%=a;，但是由于要整行操作所以先乘再减而不直接用取模。注意这里的加减乘操作需要在模意义下进行。
此时 \(a>b\)，然后交换两行，此时又有 \(a<b\) ，就可以进行上述操作。
直到减完某次之后 \(b=0\) 就可以结束了，这个时候消元结束。
根据辗转相除法，每次取模然后交换一定会出现 \(b=0\) 的结果，因此一定可以正常消元。
与此同时，\(a\) 也会比原来小，再去消下一行的时候运行的次数也会变小。
也可以看做每次都是在做 \(\gcd\) ，数据足够大时 \(\gcd\) 大概率会非常小，然后再求当前的数和前面的数的 \(\gcd\) 迭代次数也接近一两次就可以了。
所以均摊下来时间复杂度上不会有太大的增加。
代码略（还没打完）
