标签：
看完《算法导论》肯定会写单纯形
因为单纯形不仅好写而且《算法导论》里讲的很清楚
附赠uoj179的模板一个
1 #include<iostream> 2 #include<cstdio> 3 #include<algorithm> 4 #include<cmath> 5 #include<cstring> 6 #include<stdlib.h> 7 8 using namespace std; 9 const double eps=1e-9; 10 double a[50][50]; 11 int b[50],u[50],n,m,ty; 12 13 int read() 14 { 15 int x=0,f=1;char ch=getchar(); 16 while(ch<‘0‘||ch>‘9‘){if(ch==‘-‘)f=-1;ch=getchar();} 17 while(ch>=‘0‘&&ch<=‘9‘){x=x*10+ch-‘0‘;ch=getchar();} 18 return x*f; 19 } 20 int dcmp(double x) 21 { 22 if (fabs(x)<=eps) return 0; 23 if (x>0) return 1; else return -1; 24 } 25 26 void pivot(int x,int y) 27 { 28 swap(b[x],u[y]); 29 double k=a[x][y]; a[x][y]=1; 30 for (int i=0; i<=n; i++) a[x][i]/=k; 31 for (int i=0; i<=m; i++) 32 if (i!=x&&dcmp(a[i][y])!=0) 33 { 34 k=a[i][y]; a[i][y]=0; 35 a[i][0]+=(i?-1:1)*k*a[x][0]; 36 for (int j=1; j<=n; j++) 37 a[i][j]-=k*a[x][j]; 38 } 39 } 40 bool initial() 41 { 42 for (int i=1; i<=n; i++) u[i]=i; 43 for (int i=1; i<=m; i++) b[i]=n+i; 44 while (1) 45 { 46 int x=0,y=0; 47 for (int i=1; i<=m; i++) 48 if (dcmp(a[i][0])<0) x=i; //加break会TLE？ 49 if (!x) return 1; 50 for (int i=1; i<=n; i++) 51 if (dcmp(a[x][i])<0) y=i;//加break会TLE？ 52 if (!y) return 0; 53 pivot(x,y); 54 } 55 } 56 57 int simplex() 58 { 59 if (!initial()) return 0; 60 while (1) 61 { 62 int x=0,y=0; 63 for (int i=1; i<=n; i++) 64 if (dcmp(a[0][i])>0) {y=i; break;} 65 if (!y) return 1; 66 double mi=1e15; 67 for (int i=1; i<=m; i++) 68 if (dcmp(a[i][y])>0&&(!x||a[i][0]/a[i][y]<mi)) {mi=a[i][0]/a[i][y];x=i;} 69 if (!x) return -1; 70 pivot(x,y); 71 } 72 } 73 int main() 74 { 75 n=read(); m=read(); ty=read(); 76 for (int i=1; i<=n; i++) a[0][i]=read(); 77 for (int i=1; i<=m; i++) 78 { 79 for (int j=1; j<=n; j++) a[i][j]=read(); 80 a[i][0]=read(); 81 } 82 switch (simplex()) 83 { 84 case 1:{ 85 printf("%.8lf\n",a[0][0]); 86 if (ty) 87 { 88 for (int i=1; i<=n; i++) u[i]=0; 89 for (int i=1; i<=m; i++) if (b[i]<=n) u[b[i]]=i; 90 for (int i=1; i<=n; i++) printf("%.8lf ",u[i]?a[u[i]][0]:double(0)); 91 } 92 break; 93 } 94 case 0:puts("Infeasible");break; 95 case -1:puts("Unbounded");break; 96 } 97 system("pause"); 98 return 0; 99 } 100
我唯一不解的是时间复杂度问题，据说会被卡但实际上基本跑得飞起（下面会用实例证明）
还有模板的initialization上有一处我很不解就是为什么在48行和51行处的循环加break和不加break在时间会带来明显的差距……
（uoj上这两处加上break会TLE……）求神犇指教……
下面是实际应用，凡是最大流，费用流的问题大概都能用线性规划解决，而且会很快变裸题……
比如这个1061，很明显把每天的要求人数bi作为约束，每种志愿者雇佣数量做变量xi
也就是求最小化∑ci*xi i=1..m
并且x要满足约数条件 ∑aij*xi>=bi i=1..n, j=1..m, ai,j=(sj<=i<=ti)?1:0，且xj非负
但是这与标准线性规划刚好相反，标准应该是全是<=且最大化
当然如果你看过《算法导论》关于单纯形最优解证明的话，你就会知道在对偶线性规划下这很简单
对于标准型线性规划：最大化∑ci*xi i=1..n ∑aij*xi>=bi i=1..m, j=1..n，x非负
我们很容易转化成对偶情况：最小化∑bi*xi i=1..m, ∑aij*xi>=bi i=1..n, j=1..m，x非负
这两种线性规划最优值相等（不禁联想到最大流和最小割的关系）
于是我们直接上单纯形即可
但是也许有疑虑，n<=1000,m<=10000能过吗&……
然而答案是c++版本的单纯形跑了1212ms，而我以前pascal版本的费用流跑了1764ms，
所以单纯形还是非常非常厉害的，大胆的使用吧……
1 #include<iostream> 2 #include<cstdio> 3 #include<cstring> 4 #include<cmath> 5 #include<stdlib.h> 6 7 using namespace std; 8 const double eps=1e-9; 9 int b[11100],u[11100],n,m; 10 double a[10010][1010]; 11 12 int dcmp(double x) 13 { 14 if (fabs(x)<=eps) return 0; 15 if (x>0) return 1; else return -1; 16 } 17 18 void pivot(int x,int y) 19 { 20 swap(b[x],u[y]); 21 double k=a[x][y];a[x][y]=1; 22 for (int i=0; i<=n; i++) a[x][i]/=k; 23 for (int i=0; i<=m; i++) 24 if (i!=x&&dcmp(a[i][y])!=0) 25 { 26 k=a[i][y]; a[i][y]=0; 27 a[i][0]+=(i?-1:1)*k*a[x][0]; 28 for (int j=1; j<=n; j++) a[i][j]-=k*a[x][j]; 29 } 30 } 31 32 void simplex() 33 { 34 for (int i=1; i<=n; i++) u[i]=i; 35 for (int i=1; i<=m; i++) b[i]=i+n; 36 while (1) 37 { 38 int x=0,y=0; 39 for (int i=1; i<=n; i++) 40 if (dcmp(a[0][i])>0) {y=i;break;} 41 if (!y) break; 42 double mi=1e20; 43 for (int i=1; i<=m; i++) 44 if (dcmp(a[i][y])>0&&(!x||mi>a[i][0]/a[i][y])) {x=i; mi=a[i][0]/a[i][y];} 45 if (!x) break; 46 pivot(x,y); 47 } 48 } 49 50 int main() 51 { 52 scanf("%d%d",&n,&m); 53 for (int i=1; i<=n; i++) scanf("%lf",&a[0][i]); 54 for (int i=1; i<=m; i++) 55 { 56 int s,t,c; 57 scanf("%d%d%d",&s,&t,&c); 58 for (int j=s; j<=t; j++) a[i][j]=1; 59 a[i][0]=c; 60 } 61 simplex(); 62 printf("%d\n",(int)a[0][0]); 63 system("pause"); 64 return 0; 65 }
标签：
原文地址：http://www.cnblogs.com/phile/p/5658603.html
