Template of Simplex Algorithm

任务

$$\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()

输入

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)
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)
}
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;
}


418 I'm a teapot