uoj#179 线性规划 单纯形法の模板

题目链接：uoj#179
题目大意：

这是一道模板题。
本题中你需要求解一个标准型线性规划：
有n个实数变量x1,x2,⋯,xn和m条约束，其中第i条约束形如aij*xj≤bi ,j∈(1,n),i∈(1,m)
此外这n个变量需要满足非负性限制，即xj≥0。
在满足上述所有条件的情况下，你需要指定每个变量xj的取值，使得目标函数F=cj*xj ,j∈(1,n)的值最大。

输入格式
第一行三个正整数 n,m,t。其中t∈{0,1}。
第二行有n个整数c1,c2,⋯,cn，整数间均用一个空格分隔。
接下来m行，每行代表一条约束，其中第i行有n+1个整数ai1,ai2,⋯,ain,bi，整数间均用一个空格分隔。
输出格式
如果不存在满足所有约束的解，仅输出一行”Infeasible”。
如果对于任意的M，都存在一组解使得目标函数的值大于M，仅输出一行”Unbounded”。
否则，第一行输出一个实数，表示目标函数的最大值F。当第一行与标准答案的相对误差或绝对误差不超过10−6，你的答案被判为正确。
如果t=1，那么你还需要输出第二行，用空格隔开的n个非负实数，表示此时x1,x2,⋯,xn的取值，如有多组方案请任意输出其中一个。
判断第二行是否合法时，我们首先检验F−cjxj,j∈(1,n)是否为0，再对于所有ii，检验min{0,bi−aijxj,j∈(1,n)}是否为0。检验时我们会将其中大于0的项和不大于0的项的绝对值分别相加得到S+和S−，如果S+和S−的相对误差或绝对误差不超过10−6，则判为正确。
如果t=0，或者出现Infeasible或Unbounded时，不需要输出第二行。

题解：
单纯形法
这个人讲得很清楚。
http://www.cnblogs.com/zzqsblog/p/5457091.html

-update 2017/3/20
对了，对于是否有初始可行解的话，initialize这个函数是一种方法，但是可能会有点慢。大多时候都利用对偶原理来转化模型。额奥爷爷说当c[]都是负的时候才用对偶(?)

#include<cstdio>
#include<cstdlib>
#include<cstring>
#include<ctime>
#include<iostream>
#include<algorithm>
using namespace std;
#define maxn 25
const double eps=0.00000001,inf=1e15;

int n,m,id[maxn*2];
double a[maxn][maxn];//a[i][j]:i表第几条约束 j表第几个元素
double myabs(double x) {return x>0?x:-x;}
//a[0][i] -> ci 目标函数中第i个元素系数
//a[i][0] -> bi 第i条约束中的常数
//a[i][j] -> Aij 第i条约束中第j个元素的系数
//最大化 sigma(ci*xi),i∈N
//约束 xj=bj-sigma(aji*xi) ,j∈B

//转轴
void pivot(int l,int e)
//替入变量xe∈非基本变量(1~n)  替出变量xl∈基本变量(n+1~n+m)
{
    int tt=id[n+l];id[n+l]=id[e];id[e]=tt;
    int i,j;double t=a[l][e];a[l][e]=1;
    for (j=0;j<=n;j++) a[l][j]/=t;
    //重写其它等式：
    for (i=0;i<=m;i++)
     if (i!=l && myabs(a[i][e])>eps)
     {
        t=a[i][e];a[i][e]=0;
        for (j=0;j<=n;j++)
          a[i][j]-=a[l][j]*t;
     }
}

//初始化
//方法一：引入一个辅助线性规划 要求最大化-x0
//约束为 xj=bj-sigma(aji*xi)+x0 ,j∈B然后用x0替换bj为负的约束
//下面的是方法二：
bool initialize()
{
    while (1)
    {
        int i,j,e=0,l=0;
        for (i=1;i<=m;i++) 
            if (a[i][0]<-eps && (!l || (rand()&1))) l=i;
        if (!l) break;
        for (j=1;j<=n;j++) 
          if (a[l][j]<-eps && (!e || (rand()&1))) e=j;
        if (!e) {printf("Infeasible\n");return 0;}
        pivot(l,e);
        //在bi为负的时候，把所有基变量设为0不是一组合法的初始解
        //所以选择一个bi为负的基变量x[i+n]
        //然后在该约束右边找一个系数为正(即原系数为负)的非基变量进行转轴操作
        //如果没有系数为正显然就无解了
    }return 1;
}

//最优化
bool simplex()
{
    int i,j;
    while (1)
    {
        int l=0,e=0;double minn=inf;
        for (j=1;j<=n;j++)
         if (a[0][j]>eps) {e=j;break;}
        if (!e) break;
        //如果目标变量ci都小于0 那么最优解就是非基变量都为0
        for (i=1;i<=m;i++)
         if (a[i][e]>eps && a[i][0]/a[i][e]<minn)
          minn=a[i][0]/a[i][e],l=i;
        //在所有的式子中找出包含当前选中项（系数不为0）且最紧的一项
        if (!l) {printf("Unbounded\n");return 0;}
        //如果所有的a[i][e]都小于0，说明最优值正无穷
        pivot(l,e);
    }return 1;
}
double ans[maxn];

int main()
{
    //freopen("a.in","r",stdin);
    //freopen("a.out","w",stdout);
    srand(time(0));int t,i,j;
    scanf("%d%d%d",&n,&m,&t);
    //n个变量 m条约束
    for (i=1;i<=n;i++) scanf("%lf",&a[0][i]);
    for (i=1;i<=m;i++)
    {
        for (j=1;j<=n;j++)
          scanf("%lf",&a[i][j]);
        scanf("%lf",&a[i][0]);
    }
    for (i=1;i<=n;i++) id[i]=i;
    if (initialize() && simplex())
    {
        printf("%.8lf\n",-a[0][0]);
        if (t)
        {
            for (i=1;i<=m;i++) ans[id[n+i]]=a[i][0];
            for (i=1;i<=n;i++) printf("%.8lf ",ans[i]);
        }
    }
    return 0;
}