Classic Linear Inequality ( O( 2 ^ ( N + M ) * ( M ^ 3 + M * N ) ) )
I have given up using Simplex Method....
Given an augmented matrix, where each row, say { a, b, c, .. z }, represents a * x1 + b * x2 + .. ≤ z, we will find the maximum value for some vector of coefficient for these variables. A brute force solution is to take advantage of the corollary that maximal / minimal value is always in one of the vertices that all the linear equations forms. That is, we can enumerate each subset of all the vectors, find a solution, and update our answer with that solution.
N : number of variables
M : number of constraints
const double EPS = 1e-9; const double DINF = 1e9; struct Gauss_Solve_Matrix{ vvd mat; vd ans; Gauss_Solve_Matrix( const vvd &m ){ mat = m; ans.resize( m[ 0 ].size() - 1 ); } void Gauss_elimination(){ for( int i = 0; i < min( mat.size(), mat[ 0 ].size() - 1 ); ++i ){ if( mat[ i ][ i ] == 0 ) for( int j = i + 1; j < mat.size(); ++j ) if( mat[ j ][ i ] ){ swap( mat[ i ], mat[ j ] ); break; } if( mat[ i ][ i ] == 0 ) continue; for( int j = mat[ i ].size() - 1; j >= i; --j ) mat[ i ][ j ] /= mat[ i ][ i ]; for( int j = i + 1; j < mat.size(); ++j ) for( int k = mat[ j ].size() - 1; k >= i; --k ) mat[ j ][ k ] -= mat[ j ][ i ] * mat[ i ][ k ]; } } void back_substitution(){ for( int i = mat.size() - 1; i >= 0; --i ){ int t; // first non zero coefficient for( t = 0; t + 1 < mat[ i ].size(); ++t ) if( mat[ i ][ t ] ) break; if( t >= ( int ) mat[ i ].size() - 1 ) continue; double s = 0; for( int j = t + 1; j + 1 < mat[ i ].size(); ++j ) s += mat[ i ][ j ] * ans[ j ]; ans[ t ] = ( mat[ i ].back() - s ) / mat[ i ][ t ]; } } }; bool valid( const vvd &mat, const vd &sol ){ for( int i = 0; i < mat.size(); ++i ){ double s = 0; for( int j = 0; j + 1 < mat[ i ].size(); ++j ) s += mat[ i ][ j ] * sol[ j ]; if( s > EPS + mat[ i ].back() ) return false; } return true; } void solve(){ int M, N; cin >> M >> N; // M constrictions, N variables vvd mat( N + M + 1, vd( N + 1 ) ); for( int i = 0; i < M; ++i ) for( int j = 0; j < N; ++j ) cin >> mat[ i ][ j ]; for( int i = 0; i < M; ++i ) cin >> mat[ i ][ N ]; for( int i = M; i < M + N; ++i ) mat[ i ][ i - M ] = -1; // -a[ 0 ~ N - 1 ] <= 0 -> a[ 0 ~ N - 1 ] >= 0 vd qry( N ); for( int i = 0; i < N; ++i ) cin >> qry[ i ]; for( int i = 0; i < N + 1; ++i ) mat[ M + N ][ i ] = ( i != N ? qry[ i ] : DINF ); vd sol; double maxv = -DINF; for( int S = 1; S < ( 1 << N + M + 1 ); ++S ){ vvd smat; for( int i = 0; i < N + M + 1; ++i ) if( S >> i & 1 ) smat.push_back( mat[ i ] ); Gauss_Solve_Matrix *gsm = new Gauss_Solve_Matrix( smat ); gsm->Gauss_elimination(); gsm->back_substitution(); if( not valid( mat, gsm->ans ) ) continue; double value = 0; for( int i = 0; i < gsm->ans.size(); ++i ) value += gsm->ans[ i ] * qry[ i ]; if( upmax( maxv, value ) ) sol = gsm->ans; } if( fabs( maxv - ( -DINF ) ) < EPS ){ cout << "No solution\n"; return; } if( fabs( maxv - DINF ) < EPS ){ cout << "Infinity\n"; return; } cout << "Bounded solution\n"; cout << fixed << setprecision( 16 ); for( int i = 0; i < sol.size(); ++i ) cout << sol[ i ] << " \n"[ i + 1 == sol.size() ]; } int main(){ ios::sync_with_stdio( false ); solve(); return 0; }