Chaos Slover 补档计划 - 线性规划

2015.12.01 09:13 Tue| 4 visits oi_2016| ChaosSlover补档计划| Text

对于线性规划问题,我们常常有两种解决方式:单纯形算法和网络流(随机增量法?大雾 ><)。

单纯形算法可真是神奇啊!单次转轴操作是 O(NM) 的,并且时间有可能会被卡到指数级别,然而可以用来求解 1000*10000 的题 0.0 真是强大呢。无初始可行解的做法我还不会,wyf 的无初始解的代码比初始解都是 0 的代码长了一倍啊!

BZOJ3550 [ONTAK2010]Vacation

网络流解法

空灰冰魂的题解

像他写的那样列出来 2*n 个等式,建出如下图所示的网络。

图中每一条蓝色的边都代表着一个取值为 0 或 1 的变量 a,紫色的边代表辅助变量 t,而除源点和汇点外的每一个点则代表着一个流量平衡的等式。在跑网络流的过程中,经过一条边的流量就代表着一个变量的取值。这样我们就可以用边的容量来限制变量的取值了。选蓝色的边的时候,会产生一定的收益,而我们就是要让收益最大,因此跑费用流即可。

 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
93
#include <bits/stdc++.h>
using namespace std;

#define N 405
#define M 2010
#define inf 0x3f3f3f3f

int n, k;
int s, t, head[N], d[N], v[N], inc[N], from[N];

inline int read()
{
    int x = 0, f = 1; char ch = getchar();
    while (ch < '0' || ch > '9') { if (ch == '-') f = -1; ch = getchar(); }
    while (ch >= '0' && ch <= '9') { x = x * 10 + ch - '0'; ch = getchar(); }
    return x * f;
}

struct graph
{
    int next, to, val, cost;
    graph() {}
    graph(int _next, int _to, int _val, int _cost)
    : next(_next), to(_to), val(_val), cost(_cost) {}
} edge[M];

inline void add(int x, int y, int z, int w)
{
    static int cnt = 1;
    edge[++cnt] = graph(head[x], y, z, w);
    head[x] = cnt;
    edge[++cnt] = graph(head[y], x, 0, -w);
    head[y] = cnt;
}

queue<int> q;

bool spfa()
{
    memset(v, 0, sizeof v);
    memset(d, 0x3f, sizeof d);
    d[s] = 0; v[s] = 1; inc[s] = inf;
    q.push(s);
    while (!q.empty())
    {
        int x = q.front(), y; q.pop(), v[x] = 0;
        for (int i = head[x]; i; i = edge[i].next)
            if (edge[i].val && d[y = edge[i].to] > d[x] + edge[i].cost)
            {
                d[y] = d[x] + edge[i].cost;
                inc[y] = min(inc[x], edge[i].val); from[y] = i;
                if (!v[y]) v[y] = 1, q.push(y);
            }
    }
    return d[t] != inf;
}

int update()
{
    int x = t;
    while (x != s)
    {
        edge[from[x]].val -= inc[t];
        edge[from[x] ^ 1].val += inc[t];
        x = edge[from[x] ^ 1].to;
    }
    return d[t] * inc[t];
}

int edmonds_karp()
{
    int re = 0;
    while (spfa())
        re += update();
    return re;
}

int main()
{
    cin >> n >> k;
    s = n * 2 + 2; t = n * 2 + 3;
    for (int i = 1; i <= n; ++i)
        add(0, i, 1, -read());
    for (int i = 1; i <= n; ++i)
        add(i, n + i, 1, -read());
    for (int i = 1; i <= n; ++i)
        add(n + i, n * 2 + 1, 1, -read());
    for (int i = 0; i <= n * 2; ++i)
        add(i, i + 1, inf, 0);
    add(s, 0, k, 0); add(n * 2 + 1, t, k, 0);
    cout << -edmonds_karp() << endl;
    return 0;
}

单纯形解法

限制有如下两种:

  • 相邻 N 个 ai 之和小于等于 k。
  • ai 小于等于 1。

然后我们就可以列出来 5N+1 个限制条件,直接求解线性规划即可。

// 比网络流慢得可不是一点半点,大概是由于所有的计算都是 double 精度的吧。于是尝试全改成 float 就被卡精了。

 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
#include <bits/stdc++.h>
using namespace std;

#define M 1005
#define N 605
#define eps 1e-8
#define inf 1e10

struct lin_pro
{
    int n, m;
    double a[M][N], b[M], c[N], v;
    // a: 限制条件的系数矩阵 b: 限制条件的常数列向量
    // c: 目标函数的系数向量 v: 目标函数中的常数
    void pivot(int l, int e)
    {
        b[l] /= a[l][e];
        for (int i = 1; i <= n; ++i)
            if (i != e)
                a[l][i] /= a[l][e];
        a[l][e] = 1 / a[l][e];
        for (int i = 1; i <= m; ++i)
            if (i != l && fabs(a[i][e]) > eps)
            {
                b[i] -= a[i][e] * b[l];
                for (int j = 1; j <= n; ++j)
                    if (j != e)
                        a[i][j] -= a[i][e] * a[l][j];
                a[i][e] = -a[i][e] * a[l][e];
            }
        v += c[e] * b[l];
        for (int i = 1; i <= n; ++i)
            if (i != e)
                c[i] -= c[e] * a[l][i];
        c[e] = -c[e] * a[l][e];
    }

    double simplex()
    {
        while (true)
        {
            int e = 0;
            for (int i = 1; i <= n; ++i)
                if (c[i] > eps)
                    { e = i; break; }
            if (!e) return v;
            double temp = inf; int l = 0;
            for (int i = 1; i <= m; ++i)
                if (a[i][e] > eps && b[i] / a[i][e] < temp)
                    temp = b[i] / a[i][e], l = i;
            if (temp == inf) return inf;
            pivot(l, e);
        }
    }
} lp;

inline int read()
{
    int x = 0, f = 1; char ch = getchar();
    while (ch < '0' || ch > '9') { if (ch == '-') f = -1; ch = getchar(); }
    while (ch >= '0' && ch <= '9') { x = x * 10 + ch - '0'; ch = getchar(); }
    return x * f;
}

int n, k;

int main()
{
    cin >> n >> k;
    lp.n = n * 3;
    for (int i = 1; i <= n * 3; ++i)
        lp.c[i] = read();
    for (int i = 1; i <= n * 2 + 1; ++i)
    {
        lp.b[++lp.m] = k;
        for (int j = 0; j < n; ++j)
            lp.a[lp.m][i + j] = 1;
    }
    for (int i = 1; i <= n * 3; ++i)
        lp.b[++lp.m] = 1, lp.a[lp.m][i] = 1;
    cout << fixed << setprecision(0);
    cout << lp.simplex() << endl;
    return 0;
}

BZOJ1061 [Noi2008]志愿者招募

在排行榜上的前两名可真是吓到我了。

从前,有一种神奇的东西叫做对偶原理。对这道题而言,我们如果对每一天列一个限制条件的话,会得到一个需要最小化的目标函数,每一个限制条件还都是大于等于号的,转化为标准型后无法直接找到一组初始可行解。

对于一个线性规划,我们如果将它的限制条件的系数矩阵转置,限制条件的常数向量转置后作为目标函数的系数向量,目标函数的系数向量转置后作为限制条件的常数向量,大于等于号变小于等于号,求最小值变成求最大值(总之变个利索),就可以得到一个新的线性规划,称其为原始线性规划的对偶线性规划。

对偶原理就是互相对偶的线性规划的最优值是相同的。

回到这道题上,我们就可以将一开始列出的线性规划转化成对偶线性规划直接求解了。

 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
#include <bits/stdc++.h>
using namespace std;

#define N 1005
#define M 10005
#define inf 1e10
#define eps 1e-8

int n, m;

struct lin_pro
{
    int n, m;
    double a[M][N], b[M], c[N], v;

    void pivot(int l, int e)
    {
        b[l] /= a[l][e];
        for (int i = 1; i <= n; ++i)
            if (i != e)
                a[l][i] /= a[l][e];
        a[l][e] = 1 / a[l][e];
        for (int i = 1; i <= m; ++i)
            if (i != l && fabs(a[i][e]) > eps)
            {
                b[i] -= a[i][e] * b[l];
                for (int j = 1; j <= n; ++j)
                    if (j != e)
                        a[i][j] -= a[i][e] * a[l][j];
                a[i][e] = -a[i][e] * a[l][e];
            }
        v += c[e] * b[l];
        for (int i = 1; i <= n; ++i)
            if (i != e)
                c[i] -= c[e] * a[l][i];
        c[e] = -c[e] * a[l][e]; 
    }

    double simplex()
    {
        while (true)
        {
            int e = 0;
            for (int i = 1; i <= n; ++i)
                if (c[i] > eps)
                    { e = i; break; }
            if (!e) return v;
            double temp = inf; int l = 0;
            for (int i = 1; i <= m; ++i)
                if (a[i][e] > eps && b[i] / a[i][e] < temp)
                    temp = b[i] / a[i][e], l = i;
            if (temp == inf) return inf;
            pivot(l, e);
        }
    }
} lp;

inline int read()
{
    int x = 0, f = 1; char ch = getchar();
    while (ch < '0' || ch > '9') { if (ch == '-') f = -1; ch = getchar(); }
    while (ch >= '0' && ch <= '9') { x = x * 10 + ch - '0'; ch = getchar(); }
    return x * f;
}

int main()
{
    cin >> n >> m;
    lp.n = n; lp.m = m;
    for (int i = 1; i <= n; ++i)
        lp.c[i] = read();
    for (int i = 1; i <= m; ++i)
    {
        int l = read(), r = read();
        for (int j = l; j <= r; ++j)
            lp.a[i][j] = 1;
        lp.b[i] = read();
    }
    cout << fixed << setprecision(0);
    cout << lp.simplex() << endl;
    return 0;
}

BZOJ3112&JDFZOJ2000 [ZJOI2013]防守战线

和志愿者招募可以说基本一样吧。

 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
#include <bits/stdc++.h>
using namespace std;

#define N 10005
#define M 1005
#define inf 1e10
#define eps 1e-8

int n, m;

struct lin_pro
{
    int n, m;
    double a[M][N], b[M], c[N], v;

    void pivot(int l, int e)
    {
        b[l] /= a[l][e];
        for (int i = 1; i <= n; ++i)
            if (i != e)
                a[l][i] /= a[l][e];
        a[l][e] = 1 / a[l][e];
        for (int i = 1; i <= m; ++i)
            if (i != l && fabs(a[i][e]) > eps)
            {
                b[i] -= a[i][e] * b[l];
                for (int j = 1; j <= n; ++j)
                    if (j != e)
                        a[i][j] -= a[i][e] * a[l][j];
                a[i][e] = -a[i][e] * a[l][e];
            }
        v += c[e] * b[l];
        for (int i = 1; i <= n; ++i)
            if (i != e)
                c[i] -= c[e] * a[l][i];
        c[e] = -c[e] * a[l][e]; 
    }

    double simplex()
    {
        while (true)
        {
            int e = 0;
            for (int i = 1; i <= n; ++i)
                if (c[i] > eps)
                    { e = i; break; }
            if (!e) return v;
            double temp = inf; int l = 0;
            for (int i = 1; i <= m; ++i)
                if (a[i][e] > eps && b[i] / a[i][e] < temp)
                    temp = b[i] / a[i][e], l = i;
            if (temp == inf) return inf;
            pivot(l, e);
        }
    }
} lp;

inline int read()
{
    int x = 0, f = 1; char ch = getchar();
    while (ch < '0' || ch > '9') { if (ch == '-') f = -1; ch = getchar(); }
    while (ch >= '0' && ch <= '9') { x = x * 10 + ch - '0'; ch = getchar(); }
    return x * f;
}

int main()
{
    cin >> n >> m;
    lp.m = n; lp.n = m;
    for (int i = 1; i <= n; ++i)
        lp.b[i] = read();
    for (int i = 1; i <= m; ++i)
    {
        int l = read(), r = read();
        for (int j = l; j <= r; ++j)
            lp.a[j][i] = 1;
        lp.c[i] = read();
    }
    cout << fixed << setprecision(0);
    cout << lp.simplex() << endl;
    return 0;
}

BZOJ3265 志愿者招募加强版

这个加强版其实只是用来卡网络流的……

 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
#include <bits/stdc++.h>
using namespace std;

#define N 1005
#define M 10005
#define inf 1e10
#define eps 1e-8

int n, m;

struct lin_pro
{
    int n, m;
    double a[M][N], b[M], c[N], v;

    void pivot(int l, int e)
    {
        b[l] /= a[l][e];
        for (int i = 1; i <= n; ++i)
            if (i != e)
                a[l][i] /= a[l][e];
        a[l][e] = 1 / a[l][e];
        for (int i = 1; i <= m; ++i)
            if (i != l && fabs(a[i][e]) > eps)
            {
                b[i] -= a[i][e] * b[l];
                for (int j = 1; j <= n; ++j)
                    if (j != e)
                        a[i][j] -= a[i][e] * a[l][j];
                a[i][e] = -a[i][e] * a[l][e];
            }
        v += c[e] * b[l];
        for (int i = 1; i <= n; ++i)
            if (i != e)
                c[i] -= c[e] * a[l][i];
        c[e] = -c[e] * a[l][e]; 
    }

    double simplex()
    {
        while (true)
        {
            int e = 0;
            for (int i = 1; i <= n; ++i)
                if (c[i] > eps)
                    { e = i; break; }
            if (!e) return v;
            double temp = inf; int l = 0;
            for (int i = 1; i <= m; ++i)
                if (a[i][e] > eps && b[i] / a[i][e] < temp)
                    temp = b[i] / a[i][e], l = i;
            if (temp == inf) return inf;
            pivot(l, e);
        }
    }
} lp;

inline int read()
{
    int x = 0, f = 1; char ch = getchar();
    while (ch < '0' || ch > '9') { if (ch == '-') f = -1; ch = getchar(); }
    while (ch >= '0' && ch <= '9') { x = x * 10 + ch - '0'; ch = getchar(); }
    return x * f;
}

int main()
{
    cin >> n >> m;
    lp.n = n; lp.m = m;
    for (int i = 1; i <= n; ++i)
        lp.c[i] = read();
    for (int i = 1; i <= m; ++i)
    {
        int k = read();
        while (k--)
        {
            int l = read(), r = read();
            for (int j = l; j <= r; ++j)
                lp.a[i][j] = 1;
        }
        lp.b[i] = read();
    }
    cout << fixed << setprecision(0);
    cout << lp.simplex() << endl;
    return 0;
}

BZOJ3118 Orz the MST

知道这是一道线性规划的问题之后,不看数据范围的话它是很好想的。

首先我们知道确定在最小生成树上的边边权只会减小,其余的边边权只会增大。设边 l 的边权的改变量为 xl,若边 j 可以连通最小生成树上的边 i 两边的点,那么我们需要保证 xi+xj>=wi-wj。最终我们要使得改变量和消耗费用的乘积之和最小。这就是一个很显然的线性规划模型了。

然而看一下数据范围,发现限制条件最多会有 210000 个啊!然而据说数组开 4000 就过了 233。

至于找约束的过程呢,只有 SB 才会写深搜记了一堆再LCA。天哪!他还找环!天哪!这何止比较蠢!天哪,orz 那个用 Splay 维护凸包求 LIS 的人。

明明只需要小小地 DFS 一遍然后枚举非树边就好了 233。

  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
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
#include <bits/stdc++.h>
using namespace std;

#define N 4005
#define M 1005
#define eps 1e-8
#define inf 1e10

int n, m, u[M], v[M], w[M], ff[M], a[M], b[M];
int head[305], col[305];

inline int read()
{
    int x = 0, f = 1; char ch = getchar();
    while (ch < '0' || ch > '9') { if (ch == '-') f = -1; ch = getchar(); }
    while (ch >= '0' && ch <= '9') { x = x * 10 + ch - '0'; ch = getchar(); }
    return x * f;
}

struct lin_pro
{
    int n, m;
    double a[M][N], b[M], c[N], v;

    void pivot(int l, int e)
    {
        b[l] /= a[l][e];
        for (int i = 1; i <= n; ++i)
            if (i != e)
                a[l][i] /= a[l][e];
        a[l][e] = 1 / a[l][e];
        for (int i = 1; i <= m; ++i)
            if (i != l && fabs(a[i][e]) > eps)
            {
                b[i] -= a[i][e] * b[l];
                for (int j = 1; j <= n; ++j)
                    if (j != e)
                        a[i][j] -= a[i][e] * a[l][j];
                a[i][e] = -a[i][e] * a[l][e];
            }
        v += c[e] * b[l];
        for (int i = 1; i <= n; ++i)
            if (i != e)
                c[i] -= c[e] * a[l][i];
        c[e] = -c[e] * a[l][e];
    }

    double simplex()
    {
        while (true)
        {
            int e = 0;
            for (int i = 1; i <= n; ++i)
                if (c[i] > eps)
                    { e = i; break; }
            if (!e) return v;
            double temp = inf; int l = 0;
            for (int i = 1; i <= m; ++i)
                if (a[i][e] > eps && b[i] / a[i][e] < temp)
                    temp = b[i] / a[i][e], l = i;
            if (temp == inf) return inf;
            pivot(l, e);
        }
    }
} lp;

struct graph
{
    int next, to;
    graph() {}
    graph(int _next, int _to)
    : next(_next), to(_to) {}
} edge[M];

inline void add(int x, int y)
{
    static int cnt = 1;
    edge[++cnt] = graph(head[x], y);
    head[x] = cnt;
}

void dfs(int x, int y, int c)
{
    col[x] = c;
    for (int i = head[x]; i; i = edge[i].next)
        if (edge[i].to != y && col[edge[i].to] != c)
            dfs(edge[i].to, y, c);
}

int main()
{
    cin >> n >> m;
    lp.m = m;
    for (int i = 1; i <= m; ++i)
    {
        u[i] = read(), v[i] = read(), w[i] = read();
        ff[i] = read(), a[i] = read(), b[i] = read();
        if (ff[i]) add(u[i], v[i]), add(v[i], u[i]), lp.b[i] = b[i];
        else lp.b[i] = a[i];
    }
    for (int i = 1; i <= m; ++i)
        if (ff[i])
        {
            dfs(u[i], v[i], i);
            for (int j = 1; j <= m; ++j)
                if (!ff[j] &&
                    ((col[u[j]] == i && col[v[j]] != i) ||
                     (col[u[j]] != i && col[v[j]] == i)))
                    if (w[i] > w[j])
                        lp.c[++lp.n] = w[i] - w[j],
                        lp.a[i][lp.n] = lp.a[j][lp.n] = 1;
        }
    cout << fixed << setprecision(0);
    cout << lp.simplex() << endl;
    return 0;
}