高斯约旦WA#8

#include <cstdio>
#include <iostream>
#include <cmath>
#include <vector>
using namespace std;
#define esp 1e-16
#define ERR_NO_SOLVTION 1
#define ERR_MULTI_SOLVTION 2
#define SUCC 0
typedef  vector<vector<double> > Matrix;

// void print_matrix(const Matrix &a) {
//     int n = a.size();
//     printf("DEBUG(print_matrix): \n");
//     for (int i = 0; i < n; i++) {
//         for (int j = 0; j < n; j++)
//             printf("%.2f ", a[i][j]);
//         printf("| %.2f\n", a[i][n]);
//     }
// }

#define abs_bigger(a, b) (fabs(a) > fabs(b))
#define feq(a, b) (fabs(a - b) < esp)
#define abs_eq(a, b) (fabs(fabs(a) - fabs(b)) < esp)

// 看la和lb哪个更适合拿来做xi的行
int better(const Matrix &A, int la, int lb, int xi) {
    int n = A.size();
    // 如果他们不相等, 那么要未知数系数大的
    if (!abs_eq(A[la][xi], A[lb][xi]))
        return abs_bigger(A[la][xi], A[lb][xi]) ? la : lb;
    // 如果相等, 找后面未知数的系数小的
    for (int i = xi + 1; i < n; i++)
        if (!abs_eq(A[la][i], A[lb][i]))
            return abs_bigger(A[la][i], A[lb][i]) ? lb : la;
    return la;
}

// 第一个返回值标识成功, 无解, 多解, 第二个返回值为解
pair<int, vector<double> > GSE(Matrix A) {
    int n = A.size();
    bool multi_solution = false, no_solvtion = false;
    // printf("DEBUG(n): %d\n", n);
    // 对于每一个xj, 我们想要它放在第j行第j列
    for (int j = 0; j < n; j++) {
        // 寻找最大的系数
        int pivot_index = j;
        for (int i = j+1; i < n; i++)
            pivot_index = better(A, pivot_index, i, j);
        // printf("DEBUG(j, pi): %d %d\n", j, pivot_index);
        // 往第j行放
        for (int i = 0; i <= n; i++)
            swap(A[j][i], A[pivot_index][i]);

        // 主元系数0, bi不为0, 无解
        if (feq(A[j][j], 0) && !feq(A[j][n], 0)) {
            no_solvtion = true;
        }

        // 主元系数0, bi为0, 多解
        if (feq(A[j][j], 0) && feq(A[j][n], 0)) {
            multi_solution = true;
        }

        // printf("DEBUG1\n");
        // 这时我们把xj放到了第j行去了, 可以把j当行号用了
        // 然后去化简其他行
        for (int i = 0; i < n; i++) {
            if (i == j || A[j][j] == 0) continue;
            double div = A[i][j] / A[j][j];
            for (int k = j+1; k <= n; k++)
                A[i][k] -= A[j][k] * div;
        }
        // print_matrix(A);
    }

    if (no_solvtion)
        return make_pair(ERR_NO_SOLVTION, vector<double>());

    if (multi_solution)
        return make_pair(ERR_MULTI_SOLVTION, vector<double>());

    // printf("DEBUG\n");
    vector<double> X;
    for (int i = 0; i < n; i++)
        X.push_back(A[i][n] / A[i][i]);
    return make_pair(SUCC, X);
}

int main() {
    int n;
    cin >> n;
    Matrix A;
    for (int i = 0; i < n; i++)
        A.push_back(vector<double>(n+1));
    for (int i = 0; i < n; i++) 
        for (int j = 0; j <= n; j++)
            cin >> A[i][j];
    // print_matrix(A);
    pair<int, vector<double> > p = GSE(A);
    if (p.first == SUCC) {
        // 防止输出-0.00
        for (int i = 0; i < n; i++)
            printf("x%d=%.2f\n", i+1, feq(p.second[i], 0) ? 0 : p.second[i]);
    } else if (p.first == ERR_NO_SOLVTION) {
        printf("-1");
    } else if (p.first == ERR_MULTI_SOLVTION) {
        printf("0");
    }

    return 0;
}
/* 已过:
Input:
3
0 0 1 2
0 0 1 1
0 0 0 0
Output: -1

Input:
3
5 1 5 3 
5 4 1 3 
5 5 2 3 
Output:
x1=0.60
x2=0.00
x3=0.00

Input:
2
0 2 3
0 0 0
Output: 0
*/

#8WA: read -1, expected 0 link

好绝望啊...