Poj 地震之后（朱-刘算法，亦称Edmond Algorithm） 解题报告

题目链接：http://dsalgo.openjudge.cn/graph/6/
这是我少数几次看见以中国人的名字命名的算法（当然印象最深刻的是中国剩余定理），将我的理解记录如下：
朱-刘算法（Edmond algorithm）
首先，说明定义：
$G=<V,E>$是一个有向赋权图。
树形图指的是一个能包含G中所有节点的有向有根树。
最小树形图指的是树形图中边的权值和最小的那个树形图。
缩点：将几个节点看成一个节点，原来与几个节点有关系的边都重新连接到新的结点。
f(G, r)指的是输出树形图的函数，r为树根。
算法描述：
在开始算法之前，需要进行一些预处理：
a. 去掉所有终止点为r的弧。（必须）
b. 将所有平行（首尾连接的是同一对节点）且同向的边简化为其中权值最小的边。（可选）
c. 去掉图中的自环（首尾相同的弧）。（必须）
该算法递归定义：
a. 首先，对于每一个节点，选出以该节点为末端的具有最小权值的弧，记录得到最小权值弧的集合，记为B，p(·)表示父节点。
b. 如果B不含用任何环，那么B即为最小树形图，转至第d步；若存在孤立点，则不存在最小树形图。
c. 否则，B中至少含有一个有向环，记为C，对其进行缩点，构造出新的图，按照弧的不同进行分类操作。
记u，v为不同的节点，(u, v)为有向弧，vC为缩成的新节点，在权重和中要加上本来的环上所有最小权值弧的权值的和。
i. 若u在C内，v不再v内，则(vC, v)的权值即等于(u, c)的权值。
ii. 若u，v都不在C内，则(u，v)在新图中保持不变。
iii. 若u不在C内，v在C内，则(u, vC)的权值为(u, v)的权值 - (p(v), v)的权值。这里减去(p(v), v)的权值的原因是在后续的展开缩点的过程中，会将环中的一条弧拆开，而这条弧的权值已经加入了权重和，新弧的选择决定拆开的是环中的哪一条弧，所以将新弧的权值减去拆开的弧的权值才能在下一次递归中正确选择最小权值弧。
回到第a步。
d. 将最小树形图中的缩点逐个展开成环，减去的弧的终点存在于最小树形图中。将其余弧加入最小树形图。得到最终结果。
实现细节：
a. 在最小权值集合B中寻找有向环：
对于节点v，从节点v开始找其父节点P(v)，若最后到了根节点则v不在任何一个有向环中，否则v在一个有向环中。要注意判断v是否已经在某个有向环中。
b. 重建图的细节：
要重新标记节点序号，从缩点开始标记。

复杂的分析：
有人设计了O(E+VlnV)的实现，但是我的实现复杂度是O(EV)，一个粗略的分析是找环要O(V)，缩点要O(E)，最多可能有O(V)个环。

流程图样例：

代码如下：

#include <iostream>
#include <iomanip>
#include <cmath>
#include <cstring>

using namespace std;

const int MAXN = 110;
const int INF = 1 << 20;

struct Point {
    double x, y;
};

struct Edge {
    int sta, end;
    double weight;
};

double dist(Point one, Point two) {
    return sqrt(pow((one.x - two.x),2) + pow((one.y - two.y),2));
}

Point vertices[MAXN];
Edge edges[MAXN*MAXN];


double Directed_MST(int root, int V, int E) {
    int Pre[V], ID[V];
    double In[V];
    double ret = 0;
    while(true) {
        //find the minimal edge
        for(int i = 0; i < V; ++i) {
            In[i] = INF;
        }
        for(int i = 0; i < E; ++i) {
            int end = edges[i].end;
            int sta = edges[i].sta;
            if(edges[i].weight < In[end]) {
                In[end] = edges[i].weight;
                Pre[end] = sta;
            }
        }
        //detect if there is a spanning tree
        for(int i = 0; i < V; ++i) {
            if(i == root) continue;
            if(In[i] == INF) return -1;
        }
        //find the cycles
        int cntVex = 0;
        memset(ID, -1, sizeof(ID));
        //root has no income edge
        In[root] = 0;
        //mark the new id for each vertex
        for(int i = 0; i < V; ++i) {
            ret += In[i];
            int beg = i;
            //find a cycle contians i
            while(Pre[beg] != i && ID[beg] == -1 && beg != root) {
                beg = Pre[beg];
            }
            //if there is a cycle, mark new ID of the vertex in the cycle
            if(beg != root && ID[beg] == -1) {
                for(int j = Pre[beg]; j != beg; ++j) {
                    ID[j] = cntVex;
                }
                ID[beg] = cntVex++;
            }
        }
        //no cycles
        if(cntVex == 0) break; 
        //mark new ID for the vertices outside the cycles
        for(int i = 0; i < V; ++i) {
            if(ID[i] == -1) {
                ID[i] = cntVex++;
            }
        }
        //contract the cycles and remark the graph
        for(int i = 0; i < E; ++i) {
            int end = edges[i].end;
            edges[i].sta = ID[edges[i].sta];
            edges[i].end = ID[edges[i].end];
            if(edges[i].sta != edges[i].end) {
                edges[i].weight -= In[end];
            }
        }
        V = cntVex;
        root = ID[root];
    }
    return ret;
}

int main() {
    int n, m, x, y, one, two;
    while(cin.peek() != EOF && cin >> n >> m) {
        for(int i = 0; i < n; ++i) {
            cin >> x >> y;
            vertices[i].x = x;
            vertices[i].y = y;
        }
        for(int i = 0; i < m; ++i) {
            cin >> one >> two;
            edges[i].sta = one-1;
            edges[i].end = two-1;
            if(one != two){
                edges[i].weight = dist(vertices[one-1], vertices[two-1]);
            } else {
                edges[i].weight = INF;
            }
        }
        double ans = Directed_MST(0, n, m);
        if(ans == -1) {
            cout << "NO\n";
        } else {
            cout << fixed << setprecision(2) << ans << endl;
        }
    }
    return 0;
}