| 4 | typedef vector<int> VI; |
| 5 | const DOUBLE EPS = 1e-9; |
| 6 | struct LPSolver { |
| 7 | int m, n; |
| 8 | VI B, N; |
| 9 | VVD D; |
| 10 | LPSolver(const VVD &A, const VD &b, const VD &c) : |
| 11 | m(b.size()), n(c.size()), |
| 12 | N(n + 1), B(m), D(m + 2, VD(n + 2)) { |
| 13 | for (int i = 0; i < m; i++) for (int j = 0; j < n; j++) |
| 14 | D[i][j] = A[i][j]; |
| 15 | for (int i = 0; i < m; i++) { B[i] = n + i; D[i][n] = -1; |
| 16 | D[i][n + 1] = b[i]; } |
| 17 | for (int j = 0; j < n; j++) { N[j] = j; D[m][j] = -c[j]; } |
| 18 | N[n] = -1; D[m + 1][n] = 1; } |
| 19 | void Pivot(int r, int s) { |
| 20 | double inv = 1.0 / D[r][s]; |
| 21 | for (int i = 0; i < m + 2; i++) if (i != r) |
| 22 | for (int j = 0; j < n + 2; j++) if (j != s) |
| 23 | D[i][j] -= D[r][j] * D[i][s] * inv; |
| 24 | for (int j = 0; j < n + 2; j++) if (j != s) D[r][j] *= inv; |
| 25 | for (int i = 0; i < m + 2; i++) if (i != r) D[i][s] *= -inv; |
| 26 | D[r][s] = inv; |
| 27 | swap(B[r], N[s]); } |
| 28 | bool Simplex(int phase) { |
| 29 | int x = phase == 1 ? m + 1 : m; |
| 30 | while (true) { |
| 31 | int s = -1; |
| 32 | for (int j = 0; j <= n; j++) { |
| 33 | if (phase == 2 && N[j] == -1) continue; |
| 34 | if (s == -1 || D[x][j] < D[x][s] || |
| 35 | D[x][j] == D[x][s] && N[j] < N[s]) s = j; } |
| 36 | if (D[x][s] > -EPS) return true; |
| 37 | int r = -1; |
| 38 | for (int i = 0; i < m; i++) { |
| 39 | if (D[i][s] < EPS) continue; |
| 40 | if (r == -1 || D[i][n + 1] / D[i][s] < D[r][n + 1] / |
| 41 | D[r][s] || (D[i][n + 1] / D[i][s]) == (D[r][n + 1] / |
| 42 | D[r][s]) && B[i] < B[r]) r = i; } |
| 43 | if (r == -1) return false; |
| 44 | Pivot(r, s); } } |
| 45 | DOUBLE Solve(VD &x) { |
| 46 | int r = 0; |
| 47 | for (int i = 1; i < m; i++) if (D[i][n + 1] < D[r][n + 1]) |
| 48 | r = i; |
| 49 | if (D[r][n + 1] < -EPS) { |
| 50 | Pivot(r, n); |
| 51 | if (!Simplex(1) || D[m + 1][n + 1] < -EPS) |
| 52 | return -numeric_limits<DOUBLE>::infinity(); |
| 53 | for (int i = 0; i < m; i++) if (B[i] == -1) { |
| 54 | int s = -1; |
| 55 | for (int j = 0; j <= n; j++) |
| 56 | if (s == -1 || D[i][j] < D[i][s] || |
| 57 | D[i][j] == D[i][s] && N[j] < N[s]) |
| 58 | s = j; |
| 59 | Pivot(i, s); } } |
| 60 | if (!Simplex(2)) return numeric_limits<DOUBLE>::infinity(); |
| 61 | x = VD(n); |
| 62 | for (int i = 0; i < m; i++) if (B[i] < n) |
| 63 | x[B[i]] = D[i][n + 1]; |
nothing calls this directly
no outgoing calls
no test coverage detected