Simplex
コードについての説明
シンプレックス法(単体法)は線形計画問題を効率的に解くアルゴリズムである. 多項式時間で終了する実装方法は知られていないが実用上は高速に動作することが知られている.
実際に verify 用問題のように係数行列が疎なら $2$ 万(制約式の数) $\times 200$ (変数の数) とかでも計算できる. 行列が疎な場合にある程度特化して実装しているのでそうでない場合は通常の実装より遅いと思う.
多面体の面上で目的関数を最大にする方向の頂点へと移動していく手法(雑)で, 単体法と双璧をなす内点法は多面体内を移動する.
愚直に更新を行うと有限回のステップで終了しないようなケースが存在することが知られていて, 以下の実装ではそれを回避するために Blend の最小添字規則を用いている.
コメントアウトしたものは最小添字規則に反するピボット選択であるが実際にはこっちのほうが高速に動く. たぶんコメントアウトの方を使うと有限回ステップで終了しないような恣意的なケースを作れてしまいそうなのであまりおすすめはしない(まあそんなケースそうそうないと思うけど).
verify のため $5,6$ 問 Simplex 法で問題を解いたが $1$ 度だけ $2$ 個目のコメントアウトの方を使うとループが永遠に続く場合がテストケースに混じっていたことがあった.
(誤差に強くない + もう少し verify しないといけないとは思っている.)
(コンストラクタ)
Simplex$(A,b,c)$ : $\max c^\top x$ s.t. $Ax \le b,\ x \ge 0$ で表される最適化問題を解く(ans(もしくは infinity, none) に答えが格納).
コード
#define EPS 1e-10
// max c * x s.t. A*x <= b, x >= 0
class Simplex {
private:
using Arr = vector<double>;
using Mat = vector<vector<double> >;
int* index;
double** a;
int row, column, L;
void Set(const Mat& A, const Arr& b, const Arr& c){
infinity = none = false;
row = A.size(),column = A[0].size() + 1;
index = new int[row + column];
int i, j;
for(i = 0; i < row + column; i++) index[i] = i;
L = row;
a = new double*[row + 2];
for(i = 0; i < row + 2; i++) a[i] = new double[column + 1];
for(i = 0; i < row; i++){
for(j = 0; j < column - 1; j++) a[i][j] = -A[i][j];
a[i][column-1] = 1, a[i][column] = b[i];
if(a[L][column] > a[i][column]) L = i;
}
for(j = 0; j < column - 1; j++) a[row][j] = c[j];
a[row+1][column-1] = -1;
}
void solve(){
int E, i, j;
int* ls = new int[column + 2];
for(E = column - 1;;){
if(L < row){
swap(index[E], index[L + column]);
a[L][E] = 1 / a[L][E];
int prv = column + 1;
for(j = 0; j < column + 1; j++){
if(j != E){
a[L][j] *= -a[L][E];
if(abs(a[L][j]) > EPS) ls[prv] = j, prv = j;
}
}
ls[prv] = column + 1;
for(i = 0; i < row + 2; i++){
if(abs(a[i][E]) < EPS || i == L) continue;
for(j = ls[column + 1]; j < column + 1; j = ls[j]){
a[i][j] += a[i][E] * a[L][j];
}
a[i][E] *= a[L][E];
}
}
E = -1;
// double pre = EPS;
for(j = 0; j < column; j++){
if(E < 0 || index[E] > index[j]){
if(a[row + 1][j] > EPS || (abs(a[row + 1][j]) < EPS && a[row][j] > EPS)) E = j;
// if(a[row + 1][j] > pre) E = j, pre = a[row + 1][j];
// else if(abs(a[row + 1][j]) < EPS && a[row][j] > pre) E = j, pre = a[row][j];
}
}
if(E < 0) break;
L = -1;
for(i = 0; i < row; i++){
if(a[i][E] < -EPS){
if(L < 0) L = i;
else if(a[L][column] / a[L][E] - a[i][column] / a[i][E] < -EPS) L = i;
else if(a[L][column] / a[L][E] - a[i][column] / a[i][E] < EPS && index[L] > index[i]) L = i;
// if(L < 0 || a[L][column] / a[L][E] - a[i][column] / a[i][E] < EPS) L = i;
}
}
if(L < 0){
infinity = true;
return;
}
}
if(a[row + 1][column] < -EPS){
none = true;
return;
}
x.assign(column - 1, 0);
for(i = 0; i < row; i++){
if(index[column + i] < column - 1) x[index[column + i]] = a[i][column];
}
ans = a[row][column];
}
public:
bool infinity, none;
double ans;
Arr x;
Simplex(const Mat& A, const Arr& b, const Arr& c){
Set(A,b,c);
solve();
}
};
verify 用の問題
Atcoder : Enlarge Circles 提出コード(非想定解)