《挑战程序设计竞赛》 4.1.1 P287

/** \brief 高斯消元法
 *
 * \param 复杂度 n的3次方（n是方程数）
 * \param 
 * \return 返回答案数组
 *
 */     

const double EPS= 1e-8 ;

typedef vector<double> vec;
typedef vector<vec> mat;

//求解Ax=b，其中A是方阵
//当方程组无解或有无穷多解时，返回一个长度为0的数组
vec gauss_jordan(const mat& A,const vec& b){
    int n=A.size();
    mat B(n,vec(n+1));
    //把A复制给B
    for (int i=0;i<n;i++)
        for (int j=0;j<n;j++) B[i][j]=A[i][j];
    //把正在处理的未知数系数的绝对值最大的式子换到第i行
    for (int i=0;i<n;i++) B[i][n]=b[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
        for (int j=i+1;j<=n;j++) B[i][j]/=B[i][i];
        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];
            }
        }
    }
    vec x(n);
    //存放在右边的b就是答案
    for (int i=0;i<n;i++) x[i]=B[i][n];
    return x;
}