# cat /blog/2019-02-28-simplex-cpp-template/

Simplex（单纯形）算法模板

2019-03-01 13:00 +0800

浪费两个小时一晚上+一上午找一个能正常跑的Simplex模板之后我只想祝提供假模板的人NM$L。

以下 米青采彡 模板可以求是否有可行解、是否无界，并求出最优解。

使用方法：将$A[0][j]$设置为目标函数的参数，整个函数为 $$ f = \sum\limits_{j = 1}^{n} A[0][j] $$

第$i$个约束条件为 $$ \sum\limits_{j = 1}^n A[i][j] \leqslant A[i][0] $$

运行simplex()的结果若为-1则无解，1则无界，0则代表有解。最优解为$A[0][0]$，$v[j]$为各个参数的值。

感谢唾面自干同学指出这个模板的烂bug。

const double EPS = 1e-8;

int id[MAXN+MAXM] = {};
double v[MAXN] = {};
double a[MAXM][MAXN] = {};
 
int sgn(double x) {
  if (x < -EPS) return -1;
  return x > EPS ? 1 : 0;
}
 
void pivot(int r, int c) {
  swap(id[n + r], id[c]);
  double x = -a[r][c];
  a[r][c] = -1;
  for (int i = 0; i <= n; ++i) a[r][i] /= x;
  for (int i = 0; i <= m; ++i) {
    if (sgn(a[i][c]) && i != r) {
      x = a[i][c];
      a[i][c] = 0;
      for (int j = 0; j <= n; ++j) a[i][j] += x * a[r][j];
    }
  }
}
 
int simplex() {
  /* important: revert symbols of conditions */
  /* bug fixed thanks to TuoMianZiGan */
  for (int i = 1; i <= m; ++i) {
    for (int j = 1; j <= n; ++j) {
      a[i][j] *= -1;
    }
  }
  for (int i = 1; i <= n; ++i) id[i] = i;
  /* initial-simplex */
  while (true) {
    int x = 0, y = 0;
    for (int i = 1; i <= m; ++i) {
      if (sgn(a[i][0]) < 0) { x = i; break; }
    }
    if (!x) break;
    for (int i = 1; i <= n; ++i) {
      if (sgn(a[x][i]) > 0) { y = i; break; }
    }
    if (!y) return -1; // infeasible
    pivot(x, y);
  }
  /* solve-simplex */
  while (true) {
    int x = 0, y = 0;
    for (int i = 1; i <= n; ++i) {
      if (sgn(a[0][i]) > 0) { x = i; break; }
    }
    if (!x) break; // finished
    double w = 0, t = 0; bool f = true;
    for (int i = 1; i <= m; ++i) {
      if (sgn(a[i][x]) < 0) {
        t = -a[i][0] / a[i][x];
        if (f || t < w) {
          w = t, y = i, f = false;
        }
      }
    }
    if (!y) { return 1; } // unbounded
    pivot(y, x);
  }
  for (int i = 1; i <= n; ++i) v[i] = 0;
  for (int i = n + 1; i <= n + m; ++i) v[id[i]] = a[i - n][0];
  return 0;
}

References: