线性规划单纯形模板

#include<bits/stdc++.h>
using namespace std;
const int Maxn=110,Maxm=59;
class Simplex{
/*
功能：
接受有n个约束，m个基本变量的方程组a[0~n][0~m]
a[0][]存放需要最大化的目标函数，a[][0]存放常数
Base[]存放基本变量的id,初始为1~m
Rest[]存放松弛变量的id,初始为m+1~m+n
返回此线性规划的最小值ans
要求方案的话，Base[]中的变量值为0,Rest[]中的变量值为相应行的[0]
如果solve
返回1,说明运行正常ans是它的最大值
返回0,说明无可行解
返回-1,说明解没有最大值
测试：
m=2,n=3
double a[4][3]={
{0,1,3},
{8,-1,1},
{-3,1,1},
{2,1,-4}
};
solve=1,ans=64/3;
注意ac不了可能是eps的问题
*/
public:
static const double Inf;
static const double eps;
int n,m;
double a[Maxn][Maxm];
int Base[Maxm],Rest[Maxn];
double val[Maxm];
double ans;
void pt(){
for(int i=0;i<=n;i++){
for(int j=0;j<=m;j++)printf("%.2f ",a[i][j]);
puts("");
}
}
void pivot(int x,int y){//将第x个非基本变量和第y个基本变量调换
swap(Rest[x],Base[y]);
double tmp=-1./a[x][y];
a[x][y]=-1.;
for(int j=0;j<=m;j++)a[x][j]*=tmp;
for(int i=0;i<=n;i++){
if(i==x||fabs(a[i][y])<eps)continue;
tmp=a[i][y];
a[i][y]=0;
for(int j=0;j<=m;j++)a[i][j]+=tmp*a[x][j];
}
}
bool opt(){
while(1){
int csi=0;
for(int i=1;i<=m;i++)if(a[0][i]>eps&&(!csi||Base[i]<Base[csi]))csi=i;
if(!csi)break;
int csj=0;
double cur;
for(int j=1;j<=n;j++){
if(a[j][csi]>-eps)continue;
double tmp=-a[j][0]/a[j][csi];
if(!csj||tmp+eps<cur||(fabs(tmp-cur)<eps&&Rest[j]<Rest[csj]))csj=j,cur=tmp;
}
if(!csj)return 0;
pivot(csj,csi);
}
ans=a[0][0];
return 1;
}
bool init(){
ans=0;
for(int i=1;i<=m;i++)Base[i]=i;
for(int i=1;i<=n;i++)Rest[i]=m+i;
int cs=1;
for(int i=2;i<=n;i++)if(a[i][0]<a[cs][0])cs=i;
if(a[cs][0]>=-eps)return 1;
static double tmp[Maxm];
for(int i=0;i<=m;i++)tmp[i]=a[0][i],a[0][i]=0;
for(int i=1;i<=n;i++)a[i][m+1]=1.;
a[0][m+1]=-1.;Base[m+1]=m+n+1;
pivot(cs,++m);
opt();
m--;
if(a[0][0]<-eps)return 0;
cs=-1;
for(int i=1;i<=n;i++){
if(Rest[i]>m+n){
cs=i;
break;
}
}
if(cs>=1){
int nxt=-1;
m++;
for(int i=1;i<=m;i++)if(a[cs][i]>eps||a[cs][i]<-eps){nxt=i;break;}
pivot(cs,nxt);
m--;
}
for(int i=1;i<=m;i++){
if(Base[i]>m+n){
swap(Base[i],Base[m+1]);
for(int j=0;j<=n;j++)a[j][i]=a[j][m+1];
break;
}
}
for(int i=1;i<=m;i++)a[0][i]=0;a[0][0]=tmp[0];
for(int i=1;i<=m;i++)if(Base[i]<=m)a[0][i]=tmp[Base[i]];
for(int i=1;i<=n;i++){
if(Rest[i]<=m){
for(int j=0;j<=m;j++)a[0][j]+=tmp[Rest[i]]*a[i][j];
}
}
return 1;
}
void getval(){
for(int i=1;i<=m;i++)val[i]=0;
for(int i=1;i<=n;i++)if(Rest[i]<=m)val[Rest[i]]=a[i][0];
//for(int i=1;i<=m;i++)printf("%.2f ",val[i]);puts("");
}
int solve(){
if(!init())return 0;
if(!opt())return -1;
getval();
return 1;
}
}solver;
const double Simplex:: Inf=1e80;
const double Simplex:: eps=1e-8;
int main(){
int m,n,type;
scanf("%d%d%d",&m,&n,&type);
solver.a[0][0]=0;
for(int i=1;i<=m;i++)scanf("%lf",&solver.a[0][i]);
for(int i=1;i<=n;i++){
for(int j=1;j<=m+1;j++){
if(j==m+1)scanf("%lf",&solver.a[i][0]);
else {
scanf("%lf",&solver.a[i][j]);
solver.a[i][j]=-solver.a[i][j];
}
}
}
solver.m=m,solver.n=n;
int rep=solver.solve();
if(rep==0)puts("Infeasible");
else if(rep==-1)puts("Unbounded");
else {
printf("%.12f\n",solver.ans);
if(type==1){
for(int i=1;i<=m;i++)printf("%.12f%c",solver.val[i],i==m?'\n':' ');
}
}
}