Submission #62818570

---

Source Code Expand

#include <bits/stdc++.h>

using namespace std;
using DOUBLE = __float128;
using VI = vector<int>;
using VD = vector<double>;
using VVD = vector<VD>;

const DOUBLE EPS = 1e-10;

struct LPSolver {
    int m, n;
    VI B, N;
    VVD D;

    LPSolver(const VVD &A, const VD &b, const VD &c) :
        m(b.size()), n(c.size()), N(n + 1), B(m), D(m + 2, VD(n + 2)) {
        for (int i = 0; i < m; i++) for (int j = 0; j < n; j++) D[i][j] = A[i][j];
        for (int i = 0; i < m; i++) { B[i] = n + i; D[i][n] = -1; D[i][n + 1] = b[i]; }
        for (int j = 0; j < n; j++) { N[j] = j; D[m][j] = -c[j]; }
        N[n] = -1; D[m + 1][n] = 1;
    }

    void Pivot(int r, int s) {
        double inv = 1.0 / D[r][s];
        for (int i = 0; i < m + 2; i++) if (i != r)
        for (int j = 0; j < n + 2; j++) if (j != s)
            D[i][j] -= D[r][j] * D[i][s] * inv;
        for (int j = 0; j < n + 2; j++) if (j != s) D[r][j] *= inv;
        for (int i = 0; i < m + 2; i++) if (i != r) D[i][s] *= -inv;
        D[r][s] = inv;
        swap(B[r], N[s]);
    }

    bool Simplex(int phase) {
        int x = phase == 1 ? m + 1 : m;
        while (true) {
        int s = -1;
        for (int j = 0; j <= n; j++) {
            if (phase == 2 && N[j] == -1) continue;
            if (s == -1 || D[x][j] < D[x][s] || D[x][j] == D[x][s] && N[j] < N[s]) s = j;
        }
        if (D[x][s] > -EPS) return true;
        int r = -1;
        for (int i = 0; i < m; i++) {
            if (D[i][s] < EPS) continue;
            if (r == -1 || D[i][n + 1] / D[i][s] < D[r][n + 1] / D[r][s] ||
            (D[i][n + 1] / D[i][s]) == (D[r][n + 1] / D[r][s]) && B[i] < B[r]) r = i;
        }
        if (r == -1) return false;
        Pivot(r, s);
        }
    }

    DOUBLE Solve(VD &x) {
        int r = 0;
        for (int i = 1; i < m; i++) if (D[i][n + 1] < D[r][n + 1]) r = i;
        if (D[r][n + 1] < -EPS) {
        Pivot(r, n);
        if (!Simplex(1) || D[m + 1][n + 1] < -EPS) return -numeric_limits<DOUBLE>::infinity();
        for (int i = 0; i < m; i++) if (B[i] == -1) {
            int s = -1;
            for (int j = 0; j <= n; j++)
            if (s == -1 || D[i][j] < D[i][s] || D[i][j] == D[i][s] && N[j] < N[s]) s = j;
            Pivot(i, s);
        }
        }
        if (!Simplex(2)) return numeric_limits<DOUBLE>::infinity();
        x = VD(n);
        for (int i = 0; i < m; i++) if (B[i] < n) x[B[i]] = D[i][n + 1];
        return D[m][n + 1];
    }
};

const int MAXN = 11;

int board[MAXN][MAXN];

/*
U_1, U_2, ..., U_{2 * N * (N-1)}
ezek c-je -1

-D_i_j + D_i_j+1 - U <= -a_i_j+1 + a_i_j
 D_i_j - D_i_j+1 - U <=  a_i_j+1 - a_i_j

 D_i_j - K_i_j <= 10^12
-D_i_j - K_i_j <= -10^12

sum K <= P/Q

K >= 0
U >= 0
D'_i_j >= 0

n = 2 * N * (N-1) + N * N + N * N
m = 4 * N * (N-1) + 2 * N * N + 1
*/

const DOUBLE INF = 1000;

void solve() {
    int n; long double p, q; cin >> n >> p >> q;
    DOUBLE lim = DOUBLE(p)/q;

    lim = min(DOUBLE(1000), lim);

    for (int i = 0; i < n; i++) {
        for (int j = 0; j < n; j++) {
            cin >> board[i][j];
        }
    }

    const int NUM_VAR = 2 * n * (n-1) + n * n + n * n;
    const int NUM_INEQ = 4 * n * (n-1) + 2 * n * n + 1;

    VVD A(NUM_INEQ, VD(NUM_VAR));
    VD B(NUM_INEQ), C(NUM_VAR);

    for (int i = 0; i < 2 * n * (n-1); i++) {
        C[i] = -1;
    }

    auto get_index_d = [&](int i, int j) {
        return 2 * n * (n-1) + i*n + j;
    };

    auto get_index_k = [&](int i, int j) {
        return 2 * n * (n-1) + n*n + i*n + j;
    };

    auto get_index_u_2 = [&](int i, int j) {
        return (n-1)*i+j;
    };

    auto get_index_u_1 = [&](int i, int j) {
        return n * (n-1) + i+(n-1)*j;
    };

    for (int i = 0; i < n; i++) {
        for (int j = 0; j < n-1; j++) {
            int ptr = get_index_u_2(i, j)*2;

            A[ptr][get_index_d(i, j)] = -1;
            A[ptr][get_index_d(i, j+1)] = 1;
            A[ptr][get_index_u_2(i, j)] = -1;
            B[ptr] = -board[i][j+1] + board[i][j];

            ptr++;

            A[ptr][get_index_d(i, j)] = 1;
            A[ptr][get_index_d(i, j+1)] = -1;
            A[ptr][get_index_u_2(i, j)] = -1;
            B[ptr] = board[i][j+1] - board[i][j];
        }
    }
    for (int i = 0; i < n-1; i++) {
        for (int j = 0; j < n; j++) {
            int ptr = get_index_u_1(i, j)*2;
            
            A[ptr][get_index_d(i, j)] = -1;
            A[ptr][get_index_d(i+1, j)] = 1;
            A[ptr][get_index_u_1(i, j)] = -1;
            B[ptr] = -board[i+1][j] + board[i][j];

            ptr++;

            A[ptr][get_index_d(i, j)] = 1;
            A[ptr][get_index_d(i+1, j)] = -1;
            A[ptr][get_index_u_1(i, j)] = -1;
            B[ptr] = board[i+1][j] - board[i][j];
        }
    }

    for (int i = 0; i < n; i++) {
        for (int j = 0; j < n; j++) {
            int ptr = 4 * n * (n-1) + (i*n + j) * 2;

            A[ptr][get_index_d(i, j)] = 1;
            A[ptr][get_index_k(i, j)] = -1;
            B[ptr] = INF;

            ptr++;
            A[ptr][get_index_d(i, j)] = -1;
            A[ptr][get_index_k(i, j)] = -1;
            B[ptr] = -INF;
        }
    }

    int ptr = 4 * n * (n-1) + 2 * n * n;
    for (int i = 0; i < n; i++) {
        for (int j = 0; j < n; j++) {
            A[ptr][get_index_k(i, j)] = 1;
        }
    }
    B[ptr] = p/q;

    LPSolver lp(A, B, C);
    VD res;
    cout << fixed << setprecision(20);
    cout << (long double)-lp.Solve(res) << "\n";
    for (int i = 0; i < n; i++) {
        for (int j = 0; j < n; j++) {
            cout << (long double)((DOUBLE)board[i][j] + res[get_index_d(i, j)]-INF) << " ";
        }
        cout << "\n";
    }

    // int n, m; cin >> n >> m;
    // VVD A(m, VD(n));
    // VD B(m), C(n);
    // for (int i = 0; i < m; i++) {
    //     for (int j = 0; j < n; j++) {
    //         cin >> A[i][j];
    //     }
    //     cin >> B[i];
    // }
    // for (int i = 0; i < n; i++) cin >> C[i];
    // LPSolver lp(A, B, C);
    
}

int main() {
    ios::sync_with_stdio(0); cin.tie(0);
    int t = 1; 
    // cin >> t;
    while (t--) solve();
    return 0;
}

Submission Info

Compile Error

Main.cpp: In constructor ‘LPSolver::LPSolver(const VVD&, const VD&, const VD&)’:
Main.cpp:13:11: warning: ‘LPSolver::N’ will be initialized after [-Wreorder]
   13 |     VI B, N;
      |           ^
Main.cpp:13:8: warning:   ‘VI LPSolver::B’ [-Wreorder]
   13 |     VI B, N;
      |        ^
Main.cpp:16:5: warning:   when initialized here [-Wreorder]
   16 |     LPSolver(const VVD &A, const VD &b, const VD &c) :
      |     ^~~~~~~~
Main.cpp: In member function ‘bool LPSolver::Simplex(int)’:
Main.cpp:41:68: warning: suggest parentheses around ‘&&’ within ‘||’ [-Wparentheses]
   41 |             if (s == -1 || D[x][j] < D[x][s] || D[x][j] == D[x][s] && N[j] < N[s]) s = j;
Main.cpp:48:64: warning: suggest parentheses around ‘&&’ within ‘||’ [-Wparentheses]
   48 |             (D[i][n + 1] / D[i][s]) == (D[r][n + 1] / D[r][s]) && B[i] < B[r]) r = i;
Main.cpp: In member function ‘DOUBLE LPSolver::Solve(VD&)’:
Main.cpp:64:68: warning: suggest parentheses around ‘&&’ within ‘||’ [-Wparentheses]
   64 |             if (s == -1 || D[i][j] < D[i][s] || D[i][j] == D[i][s] && N[j] < N[s]) s = j;

Judge Result