差分约束系统

图论 最短路

---

简述

差分约束系统（system of difference constraints），是求解关于一组变数的特殊不等式组之方法。如果一个系统由n个变量和m个约束条件组成，其中每个约束条件形如 xj - xi <= bk(i,j∈[1,n],k∈[1,m]),则称其为差分约束系统(system of difference constraints)。亦即，差分约束系统是求解关于一组变量的特殊不等式组的方法。

算法

比如有这样一组不等式：

$$X1 - X2 <= 0$$ $$X1 - X5 <= -1$$ $$X2 - X5 <= 1$$ $$X3 - X1 <= 5$$ $$X4 - X1 <= 4$$ $$X4 - X3 <= -1$$ $$X5 - X3 <= -3$$ $$X5 - X4 <= -3$$

全都是两个未知数的差小于等于某个常数（大于等于也可以，因为左右乘以-1就可以化成小于等于）。这样的不等式组就称作差分约束系统。
这个不等式组要么无解，要么就有无数组解。因为如果有一组解{X1, X2, ..., Xn}的话，那么对于任何一个常数k，{X1 + k, X2 + k, ..., Xn + k}肯定也是一组解，因为任何两个数同时加一个数之后，它们的差是不变的，那么这个差分约束系统中的所有不等式都不会被破坏。

差分约束系统的解法利用到了单源最短路径问题中的三角形不等式。即对于任何一条边u -> v，都有：

$$d(v) <= d(u) + w(u, v)$$

其中d(u)和d(v)是从源点分别到点u和点v的最短路径的权值，w(u, v)是边u -> v1的权值。
显然以上不等式就是

$$d(v) - d(u) <= w(u, v)$$ .这个形式正好和差分约束系统中的不等式形式相同。于是我们就可以把一个差分约束系统转化成一张图，每个未知数Xi对应图中的一个顶点Vi，把所有不等式都化成图中的一条边。对于不等式 $$Xi - Xj <= c$$ ，把它化成三角形不等式： $$Xi <= Xj + c$$ ，就可以化成边Vj -> Vi，权值为c。最后，我们在这张图上求一次单源最短路径，这些三角形不等式就会全部都满足了，因为它是最短路径问题的基本性质嘛。
话说回来，所谓单源最短路径，当然要有一个源点，然后再求这个源点到其他所有点的最短路径。那么源点在哪呢？我们不妨自已造一个。以上面的不等式组为例，我们就再新加一个未知数X0。然后对原来的每个未知数都对X0随便加一个不等式（这个不等式当然也要和其它不等式形式相同，即两个未知数的差小于等于某个常数）。我们索性就全都写成 $$Xn - X0 <= 0$$ ，于是这个差分约束系统中就多出了下列不等式： $$X1 - X0 <= 0$$ $$X2 - X0 <= 0$$ $$X3 - X0 <= 0$$ $$X4 - X0 <= 0$$ $$X5 - X0 <= 0$$
对于这5个不等式，也在图中建出相应的边。最后形成的图如下：

图中的每一条边都代表差分约束系统中的一个不等式。现在以V0为源点，求单源最短路径。最终得到的V0到Vn的最短路径长度就是Xn的一个解。从图1中可以看到，这组解是{-5, -3, 0, -1, -4}。当然把每个数都加上10也是一组解：{5, 7, 10, 9, 6}。但是这组解只满足不等式组(1)，也就是原先的差分约束系统；而不满足不等式组(2)，也就是我们后来加上去的那些不等式。当然这是无关紧要的，因为X0本来就是个局外人，是我们后来加上去的，满不满足与X0有关的不等式我们并不在乎。
也有可能出现无解的情况，也就是从源点到某一个顶点不存在最短路径。也说是图中存在负权的圈。这一点我就不展开了，请自已参看最短路径问题的一些基本定理。
其实，对于图1来说，它代表的一组解其实是{0, -5, -3, 0, -1, -4}，也就是说X0的值也在这组解当中。但是X0的值是无可争议的，既然是以它作为源点求的最短路径，那么源点到它的最短路径长度当然是0了。因此，实际上我们解的这个差分约束系统无形中又存在一个条件： X0 = 0 也就是说在不等式组(1)、(2)组成的差分约束系统的前提下，再把其中的一个未知数的值定死。这样的情况在实际问题中是很常见的。比如一个问题表面上给出了一些不等式，但还隐藏着一些不等式，比如所有未知数都大于等于0或者都不能超过某个上限之类的。比如上面的不等式组(2)就规定了所有未知数都小于等于0。

对于这种有一个未知数定死的差分约束系统，还有一个有趣的性质，那就是通过最短路径算法求出来的一组解当中，所有未知数都达到最大值。下面我来粗略地证明一下，这个证明过程要结合Bellman-Ford算法的过程来说明。
假设X0是定死的；X1到Xn在满足所有约束的情况下可以取到的最大值分别为M1、M2、……、Mn（当然我们不知道它们的值是多少）；解出的源点到每个点的最短路径长度为D1、D2、……、Dn。
基本的Bellman-Ford算法是一开始初始化D1到Dn都是无穷大。然后检查所有的边对应的三角形不等式，一但发现有不满足三角形不等式的情况，则更新对应的D值。最后求出来的D1到Dn就是源点到每个点的最短路径长度。
如果我们一开始初始化D1、D2、……、Dn的值分别为M1、M2、……、Mn，则由于它们全都满足三角形不等式（我们刚才已经假设M1到Mn是一组合法的解），则Bellman-Ford算法不会再更新任合D值，则最后得出的解就是M1、M2、……、Mn。
好了，现在知道了，初始值无穷大时，算出来的是D1、D2、……、Dn；初始值比较小的时候算出来的则是M1、M2、……、Mn。大家用的是同样的算法，同样的计算过程，总不可能初始值大的算出来的结果反而小吧。所以D1、D2、……、Dn就是M1、M2、……、Mn。

那么如果在一个未知数定死的情况下，要求其它所有未知数的最小值怎么办？只要反过来求最长路径就可以了。最长路径中的三角不等式与最短路径中相反：

$$d(v) >= d(u) + w(u, v)$$ 也就是 $$d(v) - d(u) >= w(u, v)$$
所以建图的时候要先把所有不等式化成大于等于号的。其它各种过程，包括证明为什么解出的是最小值的证法，都完全类似。

总结

• 第一： 感觉难点在于建图
• 第二：
  ①：对于差分不等式，a - b <= c ，建一条 b 到 a 的权值为 c 的边，求的是最短路，得到的是最大值
  ②：对于不等式 a - b >= c ，建一条 b 到 a 的权值为 c 的边，求的是最长路，得到的是最小值
  ③：存在负环的话是无解
  ④：求不出最短路（dist[ ]没有得到更新）的话是任意解
• 第三： 一种建图方法：
  设x[i]是第i位置（或时刻）的值（跟所求值的属性一样），那么把x[i]看成数列，前n项和为s[n]，则x[i] = s[i] - s[i-1]；
  那么这样就可以最起码建立起类似这样的一个关系：0 <= s[i] - s[i-1] <= 1;
  其他关系就要去题目探索了

一些例题

POJ 1201

• 题意： 第一行输入n，下面输入n个限制条件，条件的格式为 ai bi ci, 0<=ai<=bi<=50000,1<=ci<=bi-ai+1.表示在线段[ai,bi]上至少选ci个点，使被选出的点的个数最少而且满足所有的限制条件，输出这个最小值。
• 分析：设s[i]指的是从所给所有集合并集的左端点到i的最少要取多少元素, 则对于每组给定的 ai bi ci, 有不等式 $$s[bi + 1] - s[ai] >= ci$$ 成立， 整理可得 $$s[ai] - s[bi + 1] <= -ci$$ 于是，可以建一条 ai 到 bi + 1 权值为 -ci 的边，最后对于每一个相邻的整点，i, i + 1，我们可以想到这里面可选的点最多有1个最少有0个，所以加边
  i, i + 1, 1
  i + 1, i, 0
  到此，建图完成。
• 代码如下

#include <queue>
#include <cstdio>
#include <iostream>
#include <algorithm>
using namespace std;
const int MAXN = 51010;
const int MAXE = 160010;
struct Edge{
    int v, w, nxt;
};
int n, m, ecnt;
bool vis[MAXN];
Edge edge[MAXE];
long long dis[MAXN];
int head[MAXN], val[MAXN], cnt[MAXN];
void init(){
    ecnt = 0;
    memset(cnt, 0, sizeof(cnt));
    memset(edge, 0, sizeof(edge));
    memset(head, -1, sizeof(head));
    fill(dis, dis + MAXN, (1 << 29));
}
void addEdge(int u, int v, int w){
    edge[ecnt].v = v;
    edge[ecnt].w = w;
    edge[ecnt].nxt = head[u];
    head[u] = ecnt++;
}
bool SPFA(int src){
    queue <int> que;
    que.push(src);
    vis[src] = true;
    dis[src] = 0;
    int u, v;
    while(!que.empty()){
        u = que.front();
        que.pop();
        vis[u] = false;
        cnt[u]++;
        if(cnt[u] > n) return false;
        for(int i = head[u]; i + 1; i = edge[i].nxt){
            v = edge[i].v;
            if(dis[v] > dis[u] + edge[i].w){
                dis[v] = dis[u] + edge[i].w;
                if(!vis[v]){
                    que.push(v);
                    vis[v] = true;
                }
            }
        }
    }
    return true;
}
int main(){
    while(scanf("%d", &n) != EOF){
        init();
        int l, r, w, s = 51050, t = 0;
        for(int i = 0; i < n; i++){
            scanf("%d%d%d", &l, &r, &w);
            l++, r++;
            addEdge(r, l - 1, -w);
            s = min(s, l);
            t = max(t, r);
        }
        for(int i = s; i <= t; i++){
            addEdge(i - 1, i, 1);
            addEdge(i, i - 1, 0);
        }
        if(SPFA(t)) printf("%lld\n", -dis[s - 1]);
        else printf("0\n");
    }
    return 0;
}

HDU 3666

• 题意：给你一个N*M的矩阵，求是否存在两列数a1,a2,a3...an 和 b1,b2.....bm使得对矩阵中的每个数进行下面的操作之后的值在[L,U]之间，操作为：a[i] * map[i][j] / b[j]。 N,M<=400
• 分析：要求是否满足 $$L <= a[i] * map[i][j] / b[j] <= R$$ 现将该式整理，使其符合正常的差分约束的等式，两边同时取 log 得到 $$log(L/mp[i][j])<=log(a[i])-log(b[j])<=log(R/mp[i][j])$$ 将 a[i] 编号置为 1--n, b[j] 的编号置为 n + 1--n + 2, 分解上面的公式得到 $$log(b[j]) - log(a[i]) <= log(mp[i][j]/L)$$ $$log(a[i]) - log(b[j]) <= log(R / mp[i][j])$$ 有此可以得到建边方法：
  i, j + n, log(mp[i][j] / L)
  j + n, i, log(R / mp[i][j])
  到此改题就被解决了
• 代码如下

#include <cmath>
#include <stack>
#include <cstdio>
#include <iostream>
#include <algorithm>
using namespace std;
const int MAXN = 840;
const int MAXE = 321000;
struct Edge{
    int v, nxt;
    double w;
};
int n, m, ecnt;
int cnt[MAXN];
bool vis[MAXN];
int head[MAXN];
Edge edge[MAXE];
double dis[MAXN];
void init(){
    ecnt = 0;
    memset(cnt, 0, sizeof(cnt));
    memset(vis, 0, sizeof(vis));
    memset(edge, 0, sizeof(edge));
    memset(head, -1, sizeof(head));
    //fill(dis, dis + MAXN, (1 << 29));
}
void addEdge(int u, int v, double w){
    edge[ecnt].v = v;
    edge[ecnt].w = w;
    edge[ecnt].nxt = head[u];
    head[u] = ecnt++;
}
bool SPFA(int src){
    stack <int> sta;
    int num = sqrt((double)(n + m + 1));
    for(src = 1; src <= n + m; src++){
        vis[src] = true;
        sta.push(src);
        dis[src] = 0;
        cnt[src]++;
    }
    while(!sta.empty()){
        int u = sta.top();
        sta.pop();
        vis[u] = false;
        for(int i = head[u]; i + 1; i = edge[i].nxt){
            int v = edge[i].v;
            if(dis[v] > dis[u] + edge[i].w){
                dis[v] = dis[u] + edge[i].w;
                if(!vis[v]){
                    sta.push(v);
                    cnt[v]++;
                    if(cnt[v] > num) return false;
                    vis[v] = true;
                }
            }
        }
    }
    return true;
}
int main(){
    double l, r;
    while(scanf("%d%d%lf%lf", &n, &m, &l, &r) != EOF){
        init();
        int v;
        for(int i = 1; i <= n; i++)
            for(int j = 1; j <= m; j++){
                scanf("%d", &v);
                addEdge(j + n, i, log(r / v));
                addEdge(i, j + n, log(v / l));
                //addEdge(j + n, i, log(v / l));
                //addEdge(i, j + n, log(r / v));
            }
        if(SPFA(1)) printf("YES\n");
        else printf("NO\n");
    }
    return 0;
}