100 人の生徒のうち、忘れ物 A, B, C, D をした人が 70 人、75 人、80 人、85 人いたとする。「1 種類のみを忘れた人数」の最大値を求めたい。
二段階単体法で求めた。整数解が出る保証はなかったが。30 人という結果となった。
実行結果
A = 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 0 1 0 1 0 1 0 1 0 1 0 1 0 1 0 1 0 0 1 1 0 0 1 1 0 0 1 1 0 0 1 1 0 0 0 0 1 1 1 1 0 0 0 0 1 1 1 1 0 0 0 0 0 0 0 0 1 1 1 1 1 1 1 1 (L228) b = 100 70 75 80 85 (L228) c = 0 -1 -1 0 -1 0 0 0 -1 0 0 0 0 0 0 0 (L228) opt = -30 (L234) res = 0 0 5 0 10 0 0 0 15 0 0 0 0 0 -2.66454e-15 70 (L235) state = 0 (L236)
つまり、
- (A, B, C, D) をすべて忘れた人：70 人
- (B) だけを忘れた人：5 人
- (C) だけを忘れた人：10 人
- (D) だけを忘れた人：15 人
ということで 30 人。
実験コード
#include <bits/stdc++.h> using namespace std; #define COUT(x) cout << #x << " = " << (x) << " (L" << __LINE__ << ")" << endl template<class T> inline bool chmax(T& a, T b) { if (a < b) { a = b; return 1; } return 0; } template<class T> inline bool chmin(T& a, T b) { if (a > b) { a = b; return 1; } return 0; } using pint = pair<int, int>; using pll = pair<long long, long long>; // 4-neighbor (or 8-neighbor) const vector<int> dx = {1, 0, -1, 0, 1, -1, 1, -1}; const vector<int> dy = {0, 1, 0, -1, 1, 1, -1, -1}; /*///////////////////////////////////////////////////////*/ // debug /*///////////////////////////////////////////////////////*/ #define COUT(x) cout << #x << " = " << (x) << " (L" << __LINE__ << ")" << endl template<class T1, class T2> ostream& operator << (ostream &s, pair<T1,T2> P) { return s << '<' << P.first << ", " << P.second << '>'; } template<class T> ostream& operator << (ostream &s, vector<T> P) { for (int i = 0; i < P.size(); ++i) { if (i > 0) { s << " "; } s << P[i]; } return s; } template<class T> ostream& operator << (ostream &s, deque<T> P) { for (int i = 0; i < P.size(); ++i) { if (i > 0) { s << " "; } s << P[i]; } return s; } template<class T> ostream& operator << (ostream &s, vector<vector<T> > P) { for (int i = 0; i < P.size(); ++i) { s << endl << P[i]; } return s << endl; } template<class T> ostream& operator << (ostream &s, set<T> P) { for (auto it : P) { s << "<" << it << "> "; } return s; } template<class T> ostream& operator << (ostream &s, multiset<T> P) { for (auto it : P) { s << "<" << it << "> "; } return s; } template<class T1, class T2> ostream& operator << (ostream &s, map<T1,T2> P) { for (auto it : P) { s << "<" << it.first << "->" << it.second << "> "; } return s; } // 単体法（最小添字規則によるナイーブな2段階単体法）, 有理数ライブラリにも適用可能 // A : n×m行列, n <= mでないと止まらない // min cx s.t. Ax = b, x >= 0 template<class T> struct Matrix { vector<vector<T> > val; Matrix(int n = 1, int m = 1) {val.clear(); val.resize(n, vector<T>(m));} Matrix(int n, int m, T x) {val.clear(); val.resize(n, vector<T>(m, x));} void init(int n, int m, T x = 0) {val.clear(); val.resize(n, vector<T>(m, x));} void resize(int n) {val.resize(n);} void resize(int n, int m, T x = 0) {val.resize(n); for (int i = 0; i < n; ++i) val[i].resize(m, x);} int size() {return val.size();} inline vector<T>& operator [] (int i) {return val[i];} friend ostream& operator << (ostream& s, Matrix<T> M) {s << endl; for (int i = 0; i < M.val.size(); ++i) s << M[i] << endl; return s;} }; const double INF = 1LL<<29; const double EPS = 1e-10; template<class T> int PreProcess(Matrix<T> &A, vector<T> &b) { int rank = 0; Matrix<T> M(A.size(), A[0].size()+1); for (int i = 0; i < A.size(); ++i) { for (int j = 0; j < A[0].size(); ++j) M[i][j] = A[i][j]; M[i][A[0].size()] = b[i]; } for (int i = 0; i < A[0].size(); ++i) { int pivot = rank; T max = 0; for (int j = rank; j < M.size(); ++j) { if (max < abs(M[j][i])) { pivot = j; max = abs(M[j][i]); } } if (max > EPS) { swap(M[pivot], M[rank]); T div = M[rank][i]; for (int j = 0; j < M[0].size(); ++j) { M[rank][j] = M[rank][j] / div; } for (int j = 0; j < M.size(); ++j) { if (j != rank && abs(M[j][i]) > EPS) { T fac = M[j][i]; for (int k = 0; k < M[0].size(); ++k) { M[j][k] = M[j][k] - M[rank][k] * fac; } } } ++rank; } } A.resize(rank); b.resize(rank); for (int i = 0; i < rank; ++i) { for (int j = 0; j < A[0].size(); ++j) A[i][j] = M[i][j]; b[i] = M[i].back(); } return rank; } enum STATE { OPTIMIZED, INFEASIBLE, UNBOUNDED }; template<class T> T Simplex(Matrix<T> A, vector<T> b, vector<T> c, vector<T> &res, STATE &state) { res = vector<T>(); PreProcess(A, b); // let A be row-fullrank // build tableau int n = A.size(), m = A[0].size(); for (int i = 0; i < n; ++i) if (b[i] < 0) {b[i] *= -1; for (int j = 0; j < m; ++j) A[i][j] *= -1; } vector<int> base(n), non(m); for (int i = 0; i < n; ++i) base[i] = m+i; for (int j = 0; j < m; ++j) non[j] = j; A.resize(n+2, n+m+1); for (int i = 0; i < n; ++i) A[i][m+i] = 1, A[i][n+m] = b[i]; for (int i = 0; i < n; ++i) { A[n][n+m] += A[i][n+m]; for (int j = 0; j < m; ++j) A[n][j] += A[i][j]; } for (int j = 0; j < m; ++j) A[n+1][j] = -c[j]; // start simplex for (int phase = 0; phase < 2; ++phase) { while (true) { int nn = -1, nb = -1; for (int i = 0; i < non.size(); ++i) { if (non[i] >= m) continue; // We cannot let slack value move to the base if (A[n][non[i]] > EPS) { if (nn == -1) nn = i; else if (non[i] < non[nn]) nn = i; // Bland's smallest subscript rule //if (chmax(FMax, A[n][non[i]])) nn = i; // max coefficient rule } } if (nn == -1) { if (phase == 1) break; // All done! if (A[n][A[0].size()-1] > EPS) { // No feasible solution! state = INFEASIBLE; return -1; } // detail postprocess of phase 0 bool ok = true; for (int i = 0; i < base.size(); ++i) { if (base[i] >= m) { ok = false; nb = i; break; // If base doesn't contain slack, go to phase 1 } } if (ok) break; for (int i = 0; i < non.size(); ++i) { if (non[i] < m && abs(A[nb][non[i]]) > EPS) { nn = i; break; // slack has base, continue phase 0 } } if (nn == -1) { break; // It can't be happened (A is row-fullrank) } } int col = non[nn]; if (nb == -1) { T min_ratio = INF; for (int i = 0; i < base.size(); ++i) if (A[i][col] > EPS) { T tmp = A[i][A[0].size()-1] / A[i][col]; if (min_ratio > tmp + EPS) { min_ratio = tmp, nb = i; } else if (min_ratio >= tmp - EPS) { if (nb == -1) nb = i; else if (base[i] < base[nb]) nb = i; // Bland's smallest subscript rule } } if (nb == -1) { // It cannot happen at the 1st stage state = UNBOUNDED; return -1; } } int row = nb; T piv = A[row][col]; for (int j = 0; j < A[0].size(); ++j) A[row][j] = A[row][j] / piv; for (int i = 0; i < A.size(); ++i) { if (i == row) continue; T pn = A[i][col]; for (int j = 0; j < A[0].size(); ++j) A[i][j] = A[i][j] - A[row][j] * pn; } swap(base[nb], non[nn]); } if (phase == 0) { swap(A[n], A[n+1]); A.val.pop_back(); for (int i = 0; i < n; ++i) for (int j = 0; j < A.size(); ++j) A[j].erase(A[j].begin() + m); /* for (int i = 0; i < base.size(); ++i) { // base slack can be removed if (base[i] >= m) { A.val.erase(A.val.begin() + i); int eraseID = base[i]; base.erase(base.begin() + i); --i; --n; } } */ for (int i = 0; i < non.size(); ++i) if (non[i] >= m) non.erase(non.begin() + i--); } } res = vector<T>(m, 0); for (int i = 0; i < base.size(); ++i) res[base[i]] = A[i].back(); state = OPTIMIZED; return A[n].back(); } int main() { // A : n×m行列, n <= mでないと止まらない // min cx s.t. Ax = b, x >= 0 Matrix<double> A(5, 16, 0); vector<double> b(5, 0), c(16, 0); for (int i = 0; i < 16; ++i) A[0][i] = 1; b = {100, 75, 80, 85, 90}; for (int i = 0; i < 4; ++i) { for (int bit = 0; bit < 16; ++bit) { if (bit & (1 << i)) A[i+1][bit] = 1; } } for (int bit = 0; bit < 16; ++bit) { if (__builtin_popcount(bit) == 1) c[bit] = -1; } COUT(A); COUT(b); COUT(c); vector<double> res; STATE state; double opt = Simplex(A, b, c, res, state); COUT(opt); COUT(res); COUT(state); // 70, 75, 80, 85 // (2)：5 // (3)：10 // (4)：15 // (1, 2, 3, 4)：70 }
