差分约束系统

差分约束系统，一听这个名字，顿时觉得好高大上。其实并不是这样，那么什么是差分约束系统呢？简单来说就是，给你一堆诸如：x - y <= k这种形式的不等式组。问你有无可行解，有的话就给出可行解。现在大家应该能够知道什么是差分约束系统了吧。

1、差分约束系统

现在我们给出一个标准的查分约束的定义。

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

我们将上述的m个约束条件与n个变量，组成一个m * n矩阵A。那么矩阵A将会有如下性质：

每一行有且仅有包含一个1和-1，其他所有项均为0。

我们不妨令x为n维变量(包含n个变量)，b为m维向量(包含m个bk的集合，bk同上所述)。那么我们可以将差分约束系统表示为Ax <= b，其中每一个约束条件为：xj - xi <= bk(1 <= i, j <= n, i != j且 1 <= k <= m)。

例如我们考虑需找一个满足下列条件的5维向量x = (xi)的问题：

于是我们可以将问题转化为寻找下列的8个差分约束条件的变量x1,x2,x3,x4,x5的取值：

x1-x2≤0

x1-x5≤-1

x2-x5≤1

x3-x1≤5

x4-x1≤4

x4-x3≤-1

x5-x3≤-3

x5-x4≤-3

在这问题的一个可能答案是 x = (-5, -3, 0, -1, -4)，也可能是x = (0, 2, 5, 4, 1)。二者看似不同，其实也很大的联系。

引理：设x=(x1,x2,…,xn)是差分约束系统Ax≤b的一个解，d为任意常数。则x+d=(x1+d,x2+d,…,xn+d)也是该系统Ax≤b的一个解。

2、约束图

在说约束图之前，我们先来回忆一下有向图的关联矩阵表示：①每一列恰好有一个1和-1 ②第i行1的个数等于d+(vi),-1的个数等于d-(vi).且矩阵中所有1的个数等于-1的个数等于m。我们容易发现有向图的关联矩阵和我们的约束矩阵A的转置矩阵完全一致，那么现在我们考虑能否将差分约束系统用一副图来表示呢？

我们不妨这么想，将矩阵A看作一张由n个节点和m条边构成的有向图的关联矩阵的转置。

现在我们回过头来观察差分约束系统的约束条件：xj - xi <= bk。我们将其转化为：xj <= xi + bk。总感觉这个式子似曾相识有木有，现在我们看看最短路中的松弛操作。dist[vj] <= dist[vi] + w(vj, vi)。我们又一次的发现，两者有着高度的一致。于是我们将差分约束系统中的约束条件表示为矩阵A对应图的边的权值。

现在我们归纳下上述的想法，给定差分约束系统Ax ≤ b,其对应的约束图是一个带权值的有向图G = (V, E)。

V = {v0, v1, ... , vn}
E = {(vi, vj): xj - xi ≤ bk}∪{(v0, v1), ...,(v0,vn)}

我们发现约束图中多了一个节点v0，用来保证图中至少存在一个节点，从其出发到达其他所有的结点。由于v0的特殊性，我们定义所有从v0出发的边的权值都为0。

定理：给定一差分约束系统Ax ≤ b，设G=(V, E)为该差分约束系统对应的图。如果G不包含负权回路，那么x=( d(v0,v1) , d(v0,v2) , … , d(v0,vn) )是此系统的一可行解，其中d(v0,vi)是约束图中v0到vi的最短路径(i=1,2,…,n).如果G包含负权回路，那么此系统不存在可行解。(证明略，有兴趣的同学参见算法导论)

3、求解差分约束系统

我们既然已经将差分约束系统转化为最短路问题了，那么现在就很简单了。如果该约束图中存在负权环，那么我们就只能使用Bellman-Ford或者SPFA了。理解到这里我们已经能够解决简单的差分约束系统的问题了，我们先看看POJ 3169 Layout。这一题可以说直接套定义就能够AC了。然后依次可以尝试POJ 3159 Candies(这一题有坑，用Bellman-Ford超时没有疑问，SPFA这里不能使用stl的queue优化，否则超时，可以使用数组模拟队列或者使用stl中的stack。最后的最后，千万不要用vector建图，否则一直超时，我也不知道为什么，使用链式前向星就可以ac了)。然后有点难度的就是poj1364 King和poj1201 Intervals。

注意：

①如果是求最大值，那么我们就想办法将其变成x - y <= k。如果是x - y < k，那么我们就将他变为 x - y <= k-1。然后再建图。如果是求最小值，同理。

②如果权值为正，那么Dijkstra，bellman-ford， spfa都可以。如果权值为负，我们将不能使用Dijkstra。

最后附上poj 3169 Layout的ac代码：

#include <cstdio>
#include <algorithm>
#define MAX_ML  (10000)
#define MAX_MD  (10000)
#define MAX_N   (1000)
#define INF     (10000001)
using namespace std;
///输入
int N, ML, MD;
int AL[MAX_ML], BL[MAX_ML], DL[MAX_ML];
int AD[MAX_MD], BD[MAX_MD], DD[MAX_MD];
///最短距离
int d[MAX_N];

void solve(){
    fill(d, d+N, INF);
    d[0] = 0;
    for(int k = 0; k < N; k++){<span style="white-space:pre">	</span>//bellman-ford求解
        for(int i = 0; i + 1 < N; i++){
            if(d[i + 1] < INF){
                d[i] = min(d[i], d[i+1]);
            }
        }
        for(int i = 0; i < ML; i++){
            if(d[AL[i]-1] < INF){
                d[BL[i]-1] = min(d[BL[i]-1], d[AL[i]-1]+DL[i]);
            }
        }
        for(int i = 0; i < MD; i++){
            if(d[BD[i]-1] < INF){
                d[AD[i]-1] = min(d[AD[i]-1], d[BD[i]-1]-DD[i]);
            }
        }
    }
    int res = d[N-1];
    if(d[0] < 0){
        res = -1;
    }
    else if(res == INF){
        res = -2;
    }
    printf("%d\n", res);
}

int main(){
    scanf("%d%d%d", &N, &ML, &MD);
    for(int i = 0; i < ML; i++){
        scanf("%d%d%d", &AL[i], &BL[i], &DL[i]);
    }
    for(int i = 0; i < MD; i++){
        scanf("%d%d%d", &AD[i], &BD[i], &DD[i]);
    }
    solve();
    return 0;
}

最后的最后在附上poj 3159 Candies的ac代码（原因是它实在是太坑了）:

spfa版的：

#include <cstdio>
#include <algorithm>
#include <climits>
#include <stack>
#define MAX_M 150000 + 5
#define MAX_N 30000+10
using namespace std;
int n, m;
struct edge{int to, cost, next;};
int head[MAX_N];
edge graph[MAX_M];
int d[MAX_N], tot, vis[MAX_N];
void addedge(int from, int to, int cost){
    graph[tot].to = to;
    graph[tot].cost = cost;
    graph[tot].next = head[from];
    head[from] = tot++;
}
//spfa
void solve(){
    fill(d, d+n, INT_MAX);
    d[0] = 0;
    stack<int> s;
    s.push(0);
    vis[0] = 1;
    while(!s.empty()){
        int v = s.top(); s.pop();
        vis[v] = 0;
        for(int u = head[v]; u != -1; u = graph[u].next){
            edge &e = graph[u];
            if(d[e.to] > d[v] + e.cost){
                d[e.to] = d[v] + e.cost;
                if(!vis[e.to]){
                    s.push(e.to);
                    vis[e.to] = 1;
                }
            }
        }
    }
    printf("%d\n", d[n-1]);
}

int main(){
    int a, b, c;
    scanf("%d%d", &n, &m);
    fill(vis, vis + MAX_N, 0);
    fill(head, head+MAX_N, -1);
    tot = 0;
    for(int i = 0; i < m; i++){
        scanf("%d%d%d", &a, &b, &c);
        addedge(a-1, b-1, c);
    }
    solve();
    return 0;
}

Dijkstra版的：

#include <cstdio>
#include <algorithm>
#include <climits>
#include <queue>
#define MAX_M 150000 + 5
#define MAX_N 30000+10
using namespace std;
int n, m;
struct edge{
    int to, cost, next;
    edge(int tt = 0, int tc = 0, int tn = -1):
        to(tt), cost(tc), next(tn){}
    bool operator < (edge & a) const{
        return this->cost > a.cost;
    }
};
typedef pair<int, int> P;
int head[MAX_N];
edge graph[MAX_M];
int d[MAX_N], tot, vis[MAX_N];
void addedge(int from, int to, int cost){
    graph[tot] = edge(to, cost, head[from]);
    head[from] = tot++;
}
//优先队列优化的Dijkstra
void solve(){
    fill(d, d+n, INT_MAX);
    d[0] = 0;
    priority_queue<P, vector<P>, greater<P> > Q;
    Q.push(P(0, 0));
    while(!Q.empty()){
        P p = Q.top(); Q.pop();
        int v = p.second;
        if(d[v] < p.first) continue;
        for(int i = head[v]; i != -1; i = graph[i].next){
            int to = graph[i].to;
            int cost = graph[i].cost;
            if(d[to] > d[v] + cost){
                d[to] = d[v] + cost;
                Q.push(P(d[to], to));
            }
        }
    }
    printf("%d\n", d[n-1]);
}

int main(){
    int a, b, c;
    scanf("%d%d", &n, &m);
    fill(vis, vis + MAX_N, 0);
    fill(head, head+MAX_N, -1);
    tot = 0;
    for(int i = 0; i < m; i++){
        scanf("%d%d%d", &a, &b, &c);
        addedge(a-1, b-1, c);
    }
    solve();
    return 0;
}