Template of Simplex Algorithm

任务

给定一个$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;
}