學習了《挑程》中的高斯消元法，它是求解最基礎的線性方程組（未知數個數和方程個數相等，並且有唯一解）的演算法。

首先舉一個例子：求解如下方程組：

我們手算一下這個方程組，過程如下：
1. ② - 4*①；③ - 7*①，得到如下式子：

我們這是平時手算的過程，但是一旦學習過線性代數這門課程，我們就不用再寫方程式了，直接寫矩陣的形式，還是這個例子，我們便可以寫成如下矩陣：

此矩陣為增廣矩陣，豎線之前的矩陣為方程組的係數矩陣 A，而豎線之後的列矩陣為常數矩陣 b。
所以方程組的運算就是對增廣矩陣進行相關操作，而在我們的高階語言中，可以通過二維陣列來表示此矩陣，當然也可以使用vector來表示。

具體程式碼如下：

#include<cstdio>
#include<vector>
#include<cmath>
#include<algorithm>
using namespace std;
const double EPS = 1E-8;
typedef vector<double> vec;
typedef vector<vec> mat;

vec gauss_jordan(const mat &A, const vec &b)
{
    int n = A.size();           //n個未知數 
 

    mat B(n, vec(n+1));         //初始化增廣矩陣 
    for(int i=0; i<n; i++)
    { 
        for(int j=0; j<n; j++)
        { 
            B[i][j] = A[i][j];  //係數矩陣 
        }
    }
    for(int i=0; i<n; i++)
    {
        B[i][n] = b[i];         //常數矩陣 
    }   

    //把正在處理的未知數係數的絕對值最大的式子換到第i行 
    for(int i=0; i<n; i++)
    {
        int pivot = i;
        for(int j=i; j<n; j++)
        {
            if(abs(B[j][i]) > abs(B[pivot][i])) 
                pivot = j;  
        }   
        swap(B[i], B[pivot]);
        //無解或者有無窮解 
        if(abs(B[i][i]) < EPS) return vec();

        //把正在處理的未知數的係數變為1
        double k = B[i][i];
        for(int j=i; j<=n; j++) B[i][j] /= k;
        for(int j=0; j<n; j++)
        {
            if(i != j)
            {
                //從第j個式子中消去第i個未知數 
                for(int k=i+1; k<=n; k++)
                {
                    B[j][k] -= B[j][i] * B[i][k];
                } 
            }   
        }
//      for(int j=0; j<n; j++)
//      {
//          for(int k=0; k<n+1; k++)
//          {
//              printf("%f ",B[j][k]);
//          }   
//          printf("\n");
//      } 
//      printf("--------------------\n");
    }
    vec x(n);
    for(int i=0; i<n; i++)
    {
        x[i] = B[i][n];
    }
    return x;
}


int main()
{
    freopen("input.txt", "r", stdin);
    int n;
    scanf("%d", &n);
    mat A(n, vec(n));
    vec b(n);
    double u;
    for(int i=0; i<n; i++)
    {
        for(int j=0; j<n; j++)
        {
            scanf("%lf", &A[i][j]);
        }
    }
    for(int i=0; i<n; i++)
    {
        scanf("%lf", &b[i]);
    }
    vec c = gauss_jordan(A, b);
    for(int i=0; i<c.size(); i++)
    {
        printf("%f ",c[i]);
    }
    return 0;
}