學了將近兩個晚上。突然看懂時感覺自己是個傻逼。
寫個部落格記一下吧。

模板：洛谷P3389 【模板】高斯消元法

有 \(n\) 個形如 \(a_1x_1+a_2x_2+\cdots+a_nx_n=b\) 的方程。解方程組。
有唯一解則將其求出，否則輸出 No Solution。

把方程組用矩陣形式寫出來

\[\left( \begin{matrix} a_{1,1} & a_{1,2} & \cdots & a_{1,n} \\ a_{2,1} & a_{2,2} & \cdots & a_{2,n} \\ \vdots & \vdots & \ddots & \vdots \\ a_{n,1} & a_{n,2} & \cdots & a_{n,n} \end{matrix} \middle| \begin{matrix} b_1 \\ b_2 \\ \vdots \\ b_n \end{matrix} \right)\]

（這個東西叫增廣矩陣）

考慮選 \(x_1\) 為主元。
若所有方程中 \(x_1\) 的係數均為 0，則 \(x_1\) 是個自由未知量，方程組無唯一解。直接結束即可。

否則，選出 \(|a_{i,1}|\) 最大的第 \(i\) 行（這樣精度誤差最小，雖然我也不知道為什麼）
然後用第 \(i\) 個方程和其它所有方程做一遍加減消元，消去 \(x_1\)。

不妨設 \(i=1\)，則得到新的增廣矩陣大概是下面這樣

\[\left( \begin{matrix} 1 & 1 & 1 & \cdots & 1 \\ 0 & 1 & 1 & \cdots & 1 \\ 0 & 1 & 1 & \cdots & 1 \\ \vdots & \vdots & \vdots & \ddots & \vdots \\ 0 & 1 & 1 & \cdots & 1 \end{matrix} \middle| \begin{matrix} 1 \\ 1 \\ 1 \\ \vdots \\ 1 \end{matrix} \right)\]

0 表示該項已被消去，1 表示其它。

對 \(x_2\) 執行相同操作。注意之前選過的行不能再選。
設這次選了第 \(j(j\ne i)\) 行，以 \(j=2\) 為例，新的增廣矩陣是

\[\left( \begin{matrix} 1 & 0 & 1 & \cdots & 1 \\ 0 & 1 & 1 & \cdots & 1 \\ 0 & 0 & 1 & \cdots & 1 \\ \vdots & \vdots & \vdots & \ddots & \vdots \\ 0 & 0 & 1 & \cdots & 1 \end{matrix} \middle| \begin{matrix} 1 \\ 1 \\ 1 \\ \vdots \\ 1 \end{matrix} \right)\]

重複上述過程，最終得到

\[\left( \begin{matrix} 1 & 0 & 0 & \cdots & 0 \\ 0 & 1 & 0 & \cdots & 0 \\ 0 & 0 & 1 & \cdots & 0 \\ \vdots & \vdots & \vdots & \ddots & \vdots \\ 0 & 0 & 0 & \cdots & 1 \end{matrix} \middle| \begin{matrix} 1 \\ 1 \\ 1 \\ \vdots \\ 1 \end{matrix} \right)\]

到這裡就算是解完了。

#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，就直接擺爛（反正也消不出啥東西），直接跳到下一列

這樣得到的增廣矩陣依然是

\[\left( \begin{matrix} 1 & 0 & 0 & \cdots & 0 \\ 0 & 1 & 0 & \cdots & 0 \\ 0 & 0 & 1 & \cdots & 0 \\ \vdots & \vdots & \vdots & \ddots & \vdots \\ 0 & 0 & 0 & \cdots & 1 \end{matrix} \middle| \begin{matrix} 1 \\ 1 \\ 1 \\ \vdots \\ 1 \end{matrix} \right)\]

但這次寫著 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;
}

做完之後還是感覺這程式碼醜的一批。
不過反正能過，複雜度也沒假，就這樣吧（