任务
给定一个$n$维向量$a$、一个$m \times n$的矩阵$b$和一个$m$维向量$c$。要求构造一个$n$维向量$x$,使得在满足
$$
\forall i \in [1, m], \; \sum_{j=1}^n b_{i, j}x_j \le c_i
$$
的前提下
$$
\sum_{i=1}^n a_ix_i
$$
最大。
接口
int solve()
返回$-1$表示无解,$1$表示无界,$0$表示存在最大值。
输入
int n, m
如题。
double a[][]
$(m+1) \times (n+1)$的矩阵。第$0$行表示$a$向量,第$0$列表示$c$向量,第$1 \ldots m$行、$1 \ldots n$列表示$b$矩阵。
输出
double ans[]
当
solve()
返回$0$时,表示$n$维向量$x$。此时的最大值为a[0][0]
的相反数。
代码
/* * Life is like a diaper -- short and loaded. */ #include <algorithm> #include <cstdio> #include <cstring> static int const N = 45; static long double const EPS = 1e-8; int t; int pm(long double const &n) { return (n < 0 ? -n : n) < EPS ? 0 : n < 0 ? -1 : 1; } class Array { private: long double arr[N]; public: int size; long double &operator[](int const &n) { return arr[n]; } long double max() { return *std::max_element(arr + 1, arr + size + 1); } void operator/=(long double n) { for (int i = 0; i <= size; ++i) arr[i] /= n; } void add(Array a, long double n) { for (int i = 0; i <= size; ++i) arr[i] += a[i] * n; } }; class Simplex { private: int id[N]; Array tmp; void pivot(int e, int v) { a[e] /= a[e][v]; for (int i = 0; i <= m; ++i) if (i != e) a[i].add(a[e], -a[i][v]); id[e] = v; } int simplex(int n, int m) { for (; pm(a[0].max()) > 0;) { int col = -1, row = -1; for (int i = 1; i <= n; ++i) if (pm(a[0][i]) > 0) { col = i; break; } long double mn = 1e100, tmp; for (int i = 1; i <= m; ++i) if (pm(a[i][col]) > 0 && ((tmp = pm(a[i][0] / a[i][col] - mn)) < 0 || tmp == 0 && id[i] < id[row])) mn = a[i][0] / a[i][col], row = i; if (row == -1) return 1; pivot(row, col); } return 0; } public: int n, m; Array ans, a[N]; int solve() { for (int i = 1; i <= m; ++i) a[i][n + i] = 1, id[i] = n + i; int row = -1; long double mn = 1e100; for (int i = 1; i <= m; ++i) if (a[i][0] < mn) mn = a[i][0], row = i; if (pm(mn) < 0) { for (int i = 1; i <= n; ++i) tmp[i] = a[0][i], a[0][i] = 0; for (int i = 0; i <= m; ++i) a[i][n + m + 1] = -1; for (int i = 0; i <= m; ++i) a[i].size = n + m + 1; pivot(row, n + m + 1), simplex(n + m + 1, m); if (pm(a[0][0]) != 0) return -1; for (int i = 1; i <= m; ++i) if (id[i] == n + m + 1) { for (int j = 1; j <= n + m; ++j) if (pm(a[i][j]) != 0) { pivot(i, j); break; } break; } for (int i = 1; i <= n; ++i) a[0][i] = tmp[i]; for (int i = 1; i <= m; ++i) a[0].add(a[i], -a[0][id[i]]); } for (int i = 0; i <= m; ++i) a[i].size = n + m; int ret = simplex(n + m, m); if (!ret) for (int i = 1; i <= m; ++i) if (id[i] <= n) ans[id[i]] = a[i][0]; return ret; } } solver; int main() { scanf("%d%d%d", &solver.n, &solver.m, &t); for (int i = 1; i <= solver.n; ++i) scanf("%Lf", &solver.a[0][i]); for (int i = 1; i <= solver.m; scanf("%Lf", &solver.a[i++][0])) for (int j = 1; j <= solver.n; ++j) scanf("%Lf", &solver.a[i][j]); int res = solver.solve(); if (res == -1) puts("Infeasible"); else if (res == 1) puts("Unbounded"); else { printf("%.10Lf\n", -solver.a[0][0]); if (t == 1) { for (int i = 1; i <= solver.n; ++i) printf("%.10Lf ", solver.ans[i]); puts(""); } } return 0; }