【模板】线性规划单纯形算法

110 阅读 0 评论 73 点赞

我是靠谱客的博主传统冥王星，这篇文章主要介绍【模板】线性规划单纯形算法，现在分享给大家，希望可以做个参考。

原题为[UOJ #179]

Code

复制代码

#include <bits/stdc++.h>
#define fo(i,a,b) for(int i=a;i<=b;++i)
#define fod(i,a,b) for(int i=a;i>=b;--i)
#define N 50
using namespace std;
const double eps=1e-9;
double a[N][N],ans[N];
int n,m,t,id[N],d[N];
int rd(int k)
{
return rand()%k;
}
void pivot(int l,int e)
{
swap(id[l+n],id[e]);
double t=a[l][e];a[l][e]=1;
fo(i,0,n) a[l][i]/=t;
fo(i,0,m)
{
if(i!=l&&abs(a[i][e])>eps)
{
double t=a[i][e];a[i][e]=0;
fo(j,0,n) a[i][j]-=t*a[l][j];
}
}
}
bool prp()
{
while(1)
{
int l=0,e=0;
d[0]=0;
fo(i,1,m) if(a[i][0]<-eps) d[++d[0]]=i;
if(!d[0]) return 1;
l=d[rd(d[0])+1],d[0]=0;
fo(i,1,n) if(a[l][i]<-eps) d[++d[0]]=i;
if(!d[0]) return 0;
e=d[rd(d[0])+1];
pivot(l,e);
}
}
bool solve()
{
while(1)
{
int l=0,e=0;double mi=1e18;
fo(i,1,n)
if(a[0][i]>eps&&(!e||a[0][i]>a[0][e])) e=i;
if(!e) return 1;
fo(i,1,m)
{
if(a[i][e]>eps&&a[i][0]/a[i][e]<mi) mi=a[i][0]/a[i][e],l=i;
}
if(!l) return 0;
pivot(l,e);
}
}
int main()
{
srand(9982443);
cin>>n>>m>>t;
fo(i,1,n)
{
int x;
scanf("%d",&x);
a[0][i]=x;
}
fo(i,1,m)
{
int x;
fo(j,1,n) scanf("%d",&x),a[i][j]=x;
scanf("%d",&x),a[i][0]=x;
}
fo(i,1,n+m) id[i]=i;
if(!prp())
{
printf("Infeasiblen");
return 0;
}
if(!solve())
{
printf("Unboundedn");
return 0;
}
printf("%.10lfn",-a[0][0]);
if(t==1)
{
fo(i,1,m) ans[id[i+n]]=a[i][0];
fo(i,1,n) printf("%.10lf ",ans[i]);
}
}

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
#include <bits/stdc++.h>
#define fo(i,a,b) for(int i=a;i<=b;++i)
#define fod(i,a,b) for(int i=a;i>=b;--i)
#define N 50
using namespace std;
const double eps=1e-9;
double a[N][N],ans[N];
int n,m,t,id[N],d[N];
int rd(int k)
{
return rand()%k;
}
void pivot(int l,int e)
{
swap(id[l+n],id[e]);
double t=a[l][e];a[l][e]=1;
fo(i,0,n) a[l][i]/=t;
fo(i,0,m)
{
if(i!=l&&abs(a[i][e])>eps)
{
double t=a[i][e];a[i][e]=0;
fo(j,0,n) a[i][j]-=t*a[l][j];
}
}
}
bool prp()
{
while(1)
{
int l=0,e=0;
d[0]=0;
fo(i,1,m) if(a[i][0]<-eps) d[++d[0]]=i;
if(!d[0]) return 1;
l=d[rd(d[0])+1],d[0]=0;
fo(i,1,n) if(a[l][i]<-eps) d[++d[0]]=i;
if(!d[0]) return 0;
e=d[rd(d[0])+1];
pivot(l,e);
}
}
bool solve()
{
while(1)
{
int l=0,e=0;double mi=1e18;
fo(i,1,n)
if(a[0][i]>eps&&(!e||a[0][i]>a[0][e])) e=i;
if(!e) return 1;
fo(i,1,m)
{
if(a[i][e]>eps&&a[i][0]/a[i][e]<mi) mi=a[i][0]/a[i][e],l=i;
}
if(!l) return 0;
pivot(l,e);
}
}
int main()
{
srand(9982443);
cin>>n>>m>>t;
fo(i,1,n)
{
int x;
scanf("%d",&x);
a[0][i]=x;
}
fo(i,1,m)
{
int x;
fo(j,1,n) scanf("%d",&x),a[i][j]=x;
scanf("%d",&x),a[i][0]=x;
}
fo(i,1,n+m) id[i]=i;
if(!prp())
{
printf("Infeasiblen");
return 0;
}
if(!solve())
{
printf("Unboundedn");
return 0;
}
printf("%.10lfn",-a[0][0]);
if(t==1)
{
fo(i,1,m) ans[id[i+n]]=a[i][0];
fo(i,1,n) printf("%.10lf ",ans[i]);
}
}