单纯形法

看了集训队答辩，感觉要学习的有杜教筛高级版、线性规划、FFT、仙人掌、高级版线段树

不出意外的话一个月内博客内都不会有别的东西了QAQ

首先是喜闻乐见的单纯形法解线性规划。

今年（2016年）和线性规划有关的集训队论文有两篇，大家可以自行翻一下集训队论文（当然如果你没有拿到你可以去UOJ群下载啊），下面的大部分内容都是参阅akf那篇

线性规划的标准型一般长得像这样：

一般我们拿到的都是像标准型这样的问题，例如网络流问题，我们就是要最大化流量并且让每个点流量守恒。

假设我们把汇点到源点连一条容量为+∞的边那么所有点就都流量守恒了。我们用c(u,v)表示(u,v)这条边的容量，f(u,v)表示这条边的流量。

那么我们不难得到这样的一个标准型线性规划

不难把它转化成标准型这样的模型。

而最小费用流则也是一个标准型，我们还是把汇点到源点加一条容量为+∞，费用为0的边。

那么也可以得到一个比较标准的线性规划：

（较论文中的配文有修改）

如果你要最小费用最大流的话只要加一个约束条件$\sum_{(u,v)}{f(u,v)}=maxflow$

额那不等号怎么办呢？

对于不等约束条件${\sum_{j=1}^n{a_{ij}x_j}}\leq{b_i}$

在松弛型等式左边的那些变量我们把它们叫做基变量，在上面的松弛型表示中就是x[n+1]…x[n+m]，非基变量就是在等式右边的那些，在上面的表示中是x[1]…x[n]。

显然当bi非负时令所有基变量为bi，非基变量为0，即可得到一个满足条件的初始解。（至于bi可能为负的情况下面在单独说）

单纯形法有一个重要的操作，叫转轴(pivot)操作。就是说我们可以把一个基变量x[b]和非基变量x[n]互换，用x[b]和其他非基变量代换这个x[n]，这样x[b]就成了非基变量，而x[n]成了基变量。

一开始我们知道。

那么简单代换一下会发现，然后我们也可以把其他约束条件中的x[n]这样代换掉，于是就完成了这个转轴操作。

然后有了这个转轴过程的帮助，我们就可以实现线性规划算法了。

首先为了最大化目标函数，我们考虑不停地找一个系数为正的非基变量，然后增大这个x。

具体的这个最优化过程如下：

（注意邹逍遥原论文在这里的做法是选择满足的，ysy告诉我这样是不科学的，akf的论文里写的也是ce>0，于是做了一些修改）

这是为什么呢？akf在论文中给出了一个例子，从中我们可以大概理解一下。

这个线性规划是这样的：

我们考虑第一步我们选择换正系数的x1，因为增大x1必然会增大目标函数。

我们分别考虑三个限制，显然第一个限制下x1<=30，第二个限制下x1<=12，第三个限制x1<=9，然后x1<=9的下界最紧，所以我们用x1代换x6，得到下面的新线性规划：

哇目标函数增大到了27，可喜可贺。

接下来假设我们来增大x3，类似地可以得到：

然后可以增大x2，可以得到：

由于所有非基变量系数均为负的，所以我们不可能再得到更小的解了。

咦，有没有注意到上面漏了什么东西说要下面讲啊…

对了，在bi为负的时候，如果你把所有基变量带0，而非基变量就不一定是一组合法的初始解。

这时候我们就需要一个初始化操作，初始化操作的基本思想是引入一个辅助线性规划：

如果我们求得了这个辅助线性规划的最优解，那么如果最优解中x0>0显然原线性规划无解。

如果最优解中x0=0，那么x[1]…x[n+m]就是原线性规划的一个可行解。

而我们容易构造出辅助线性规划的一个可行解。

我们只要把x0作为换入变量，找到bi的最小值bl，把x[i+l]用x0来替代。

例如：

我们引入辅助线性规划：

bi最小值是-4，那么我们把x4用x0来替代：

然后我们就可以得到一组x的初始解用来进行最优化操作。

讲道理：为什么bi这样之后一定为正？

首先bi的最小值bl一定为负（废话，否则就不用进行最优化操作了）

然后我们考虑第l个约束会变为：

而-bi显然为正的。

而其他约束会变为：

由于bl是最小值，所以bi>=bl，所以-bl+bi>=0。

而我们发现UOJ上似乎有一种更加神奇的初始化方法？

我们的目标显然是让所有系数bi都为非负的。

我们选择一个bi为负的基变量x[i+n]，然后我们在该约束右边找一个系数为正的非基变量，然后进行转轴操作，如果没有系数为正显然就无解了。

额其实这和这个初始化操作本质上是一样的。

所以这样就可以完成整个线性规划的过程。

我们来分析一下时间复杂度？

pivot操作的复杂度显然是O(NM)的，但是最优化操作中pivot操作的调用次数可能会成为指数级。

但是我们可以发现要达到这个指数级的调用次数，边权也必须是指数级的。所以在OI中往往跑得比谁都快。所以“能在1s内跑出范围为几百的数据”。

然而一般线性规划是可以在多项式范围内求解的，只不过…

椭球算法 O(n^6*m^2)

内点算法 O(n^3.5*m^2)

改进的内点算法 O(n^3.5*m)

在oi中一般并不能有什么卵用

所以单纯形法就这么讲完啦。

实现的时候可以发现，写单纯形并不要真的把松弛型建出来，只要假装第i个约束条件就对应第i个基变量，把系数取负就可以了。

科学的代码 http://uoj.ac/submission/61872 (worse)

我的代码（有一些细节写的比较丑

#include <iostream>
#include <stdio.h>
#include <stdlib.h>
#include <algorithm>
#include <string.h>
#include <vector>
#include <math.h>
#include <limits>
#include <set>
#include <map>
using namespace std;
typedef double ld;
#define SZ 233
int n,m,t,id[SZ];
ld a[SZ][SZ],vv[SZ];
const ld eps=1e-7;
int dcmp(ld x)
{
    if(x<-eps) return -1;
    if(x>eps) return 1;
    return 0;
}
void pivot(int r,int c) //基变量r和非基变量c 
{
    swap(id[r+n],id[c]);
    ld x=-a[r][c]; a[r][c]=-1;
    for(int i=0;i<=n;i++) a[r][i]/=x;
    for(int i=0;i<=m;i++)
    {
        if(dcmp(a[i][c])&&i!=r);else continue;
        x=a[i][c]; a[i][c]=0;
        for(int j=0;j<=n;j++) a[i][j]+=x*a[r][j];
    }
}
void solve()
{
    for(int i=1;i<=n;i++) id[i]=i;
    while(1) //init
    {
        int x=0,y=0;
        for(int i=1;i<=m;i++)
        {
            if(dcmp(a[i][0])<0&&(!x||(rand()&1))) x=i; 
        }
        if(!x) break;
        for(int i=1;i<=n;i++)
        {
            if(dcmp(a[x][i])>0&&(!y||(rand()&1))) {y=i; break;}
        }
        if(!y) {puts("Infeasible"); return;}
        pivot(x,y);
    }
    while(1) //simplex
    {
        int x=0,y=0;
        for(int i=1;i<=n;i++)
        {
            if(dcmp(a[0][i])>0&&(!x||(rand()&1))) x=i; 
        }
        if(!x) break;
        double w,t; bool f=1;
        for(int i=1;i<=m;i++)
        {
            if(dcmp(a[i][x])<0&&((t=-a[i][0]/a[i][x]),f||t<w)) w=t, y=i, f=0;
        }
        if(!y) {puts("Unbounded"); return;}
        pivot(y,x);
    }
    printf("%.9lf\n",a[0][0]);
    if(!t) return;
    for(int i=1;i<=n;i++) vv[i]=0;
    for(int i=n+1;i<=n+m;i++) vv[id[i]]=a[i-n][0];
    for(int i=1;i<=n;i++) printf("%.9lf ",vv[i]);
}
int main()
{
    scanf("%d%d%d",&n,&m,&t);
    for(int i=1;i<=n;i++) scanf("%lf",&a[0][i]);
    for(int i=1;i<=m;i++)
    {
        for(int j=1;j<=n;j++) scanf("%lf",&a[i][j]), a[i][j]*=-1;
        scanf("%lf",&a[i][0]);
    }
    solve();
}