for (int j = 0; j < n-1; j++) {
        int k;
        double max_val = -1e100;
        for (int i = j; i < n; i++) {
            if (max_val < a[i][j]) {
                max_val = a[i][j];
                k = i;
            }
        }

        swap(a[j], a[k]);
        swap(b[j], b[k]);

        for (k = j+1; k < n; k++) {
            double m = a[k][j] / a[j][j];
            for (int i = j; i < n; i++)
                a[k][i] -= m * a[j][i];
            b[k] -= m * b[j];
        }
    }