$\newcommand{\O}{\mathrm{O}}$ My Algorithm : kopricky アルゴリズムライブラリ

kopricky アルゴリズムライブラリ

Minimum Arborescence Algorithm

コードについての説明(個人的メモ)

有向グラフと点 $r$ が与えられたときに $r$ を根とする最小全域有向木(根からすべての点に有向辺の向きに沿って到達できるような木)を求めるアルゴリズム.
元論文は "Efficient algorithm for finding minimum spanning trees in undirected and directed graphs" [Gabow, Galil, Spencer, Tarjan 1985].
多重辺や自己ループがあっても問題なく動作する.

(関数)
solve$(root)$ : $root$ を根とする全域有向木が存在するならその最小コストを返し, $parent$ に最小全域有向木が格納される. 解が存在しないなら最大値を返す.

時間計算量: $\O(n \log n + m)$

コード

template<typename _Tp>
class LazyUnionFind {
private:
    const int sz;
    vector<_Tp> diff, lazy;
    vector<int> par, nrank;
    pair<int, _Tp> _find(int x){
        if(par[x] == x) return {x, (_Tp)0};
        auto res = _find(par[x]);
        par[x] = res.first, diff[x] += res.second, lazy[x] += res.second;
        return {res.first, lazy[x]};
    }
public:
    LazyUnionFind(const int node_size)
        : sz(node_size), diff(sz, (_Tp)0), lazy(sz, (_Tp)0), par(sz), nrank(sz, 0){
        iota(par.begin(), par.end(), 0);
    }
    int find(int x){ return _find(x).first; }
    _Tp find_value(int x){
        _find(x);
        return (par[x] == x) ? diff[x] : (diff[x] + lazy[par[x]]);
    }
    void change_value(_Tp value, int x){
        diff[x] += value, lazy[x] += value;
    }
    int unite(int x, int y){
        x = find(x), y = find(y);
        if(x == y) return -1;
    	if(nrank[x] < nrank[y]) swap(x,y);
        par[y] = x, diff[y] -= lazy[x], lazy[y] -= lazy[x];
        if(nrank[x] == nrank[y]) nrank[x]++;
        return x;
    }
};

template<class edge, typename _Tp>
class FibonacciHeap
{
public:
    using data_type = const edge*;
    class node
    {
    public:
        data_type _data;
        unsigned short int _child;
        bool _mark;
        node *_par, *_prev, *_next, *_ch_last;
        node(data_type data) : _data(data), _child(0), _mark(false),
            _par(nullptr), _prev(nullptr), _next(nullptr), _ch_last(nullptr){}
        void insert(node *cur){
            if(_ch_last) insert_impl(cur, _ch_last);
            else _ch_last = cur, _ch_last->_prev = _ch_last->_next = _ch_last;
            ++_child, cur->_par = this;
        }
        void erase(node *cur){
            if(cur == cur->_prev){
                _ch_last = nullptr;
            }else{
                erase_impl(cur);
                if(cur == _ch_last) _ch_last = cur->_prev;
            }
            --_child, cur->_par = nullptr;
        }
    };

private:
    node *_minimum;
    vector<node*> rank;
    LazyUnionFind<_Tp>& uf;
    vector<FibonacciHeap*>& fh;
    vector<int>& heap;

    static void insert_impl(node *cur, node *next){
        cur->_prev = next->_prev, cur->_next = next;
        cur->_prev->_next = cur, next->_prev = cur;
    }
    static void erase_impl(node *cur){
        cur->_prev->_next = cur->_next, cur->_next->_prev = cur->_prev;
    }
    inline const _Tp get_key(node *cur) const noexcept {
        return cur->_data->cost + uf.find_value(cur->_data->to);
    }
    inline FibonacciHeap *home_heap(node *cur) const noexcept {
        return fh[uf.find(heap[uf.find(cur->_data->from)])];
    }
    void root_insert(node *cur){
        if(_minimum){
            insert_impl(cur, _minimum);
        }else{
            _minimum = cur, _minimum->_prev = _minimum->_next = _minimum;
        }
    }
    void root_erase(node *cur){
        if(cur == cur->_prev) _minimum = nullptr;
        else{
            if(cur == _minimum) _minimum = cur->_prev;
            erase_impl(cur);
        }
    }
    void _delete(node *cur){
        root_erase(cur);
        // delete cur;
    }
    node *_push(data_type e){
        node *new_node = new node(e);
        root_insert(new_node);
        return new_node;
    }
    void detect_minimum(){
        node *cand = nullptr;
        _Tp _min = numeric_limits<_Tp>::max();
        for(node *cur = _minimum->_next;;cur = cur->_next){
            _Tp value = get_key(cur);
            if(_min > value) cand = cur, _min = value;
            if(cur == _minimum) break;
        }
        _minimum = cand;
    }
    data_type _pop(node *given_minimum = nullptr){
        if(!_minimum) return nullptr;
        if(given_minimum) _minimum = given_minimum;
        else detect_minimum();
        data_type res = _minimum->_data;
        if(_minimum->_ch_last){
            for(node *cur = _minimum->_ch_last->_next;;){
                node *next = cur->_next;
                _minimum->erase(cur);
                home_heap(cur)->root_insert(cur);
                if(!_minimum->_ch_last) break;
                cur = next;
            }
        }
        node *next_minimum = _minimum->_next;
        if(next_minimum == _minimum){
            _delete(_minimum);
            return res;
        }
        for(node*& cur : rank) cur = nullptr;
        for(node *cur = next_minimum; cur != _minimum;){
            if(get_key(cur) < get_key(next_minimum)) next_minimum = cur;
            node *next = cur->_next;
            unsigned int deg = cur->_child;
            if(rank.size() <= deg) rank.resize(deg + 1, nullptr);
            while(rank[deg]){
                if(get_key(cur) < get_key(rank[deg]) || cur == next_minimum){
                    root_erase(rank[deg]), cur->insert(rank[deg]);
                }else{
                    root_erase(cur), rank[deg]->insert(cur);
                    cur = rank[deg];
                }
                rank[deg++] = nullptr;
                if(rank.size() <= deg) rank.resize(deg + 1, nullptr);
            }
            rank[deg] = cur;
            cur = next;
        }
        _delete(_minimum);
        _minimum = next_minimum;
        return res;
    }
    void _to_root(node *cur){
        if(!cur->_par) return;
        while(true){
            node *next = cur->_par;
            next->erase(cur);
            home_heap(cur)->root_insert(cur);
            cur->_mark = false, cur = next;
            if(!cur->_par) break;
            if(!cur->_mark){
                cur->_mark = true;
                break;
            }
        }
    }
    void _delete_node(node *cur){
        _to_root(cur), home_heap(cur)->_pop(cur);
    }

public:
    FibonacciHeap(LazyUnionFind<_Tp>& _uf, vector<FibonacciHeap*>& _fh, vector<int>& _heap)
        noexcept : _minimum(nullptr), uf(_uf), fh(_fh), heap(_heap){}
    node *push(data_type e){ return _push(e); }
    data_type pop(){ return _pop(); }
    void to_root(node *cur){ _to_root(cur); }
    void delete_node(node *cur){ _delete_node(cur); }
    friend void move_node(FibonacciHeap *fh1, FibonacciHeap::node *cur, FibonacciHeap *fh2){
        if(!cur->_par){
            fh1->root_erase(cur), fh2->root_insert(cur);
            return;
        }
        bool _first = true;
        while(true){
            node *next = cur->_par;
            next->erase(cur);
            if(_first) fh2->root_insert(cur), _first = false;
            else fh1->home_heap(cur)->root_insert(cur);
            cur->_mark = false, cur = next;
            if(!cur->_par) break;
            if(!cur->_mark){
                cur->_mark = true;
                break;
            }
        }
    }
    friend FibonacciHeap *meld(FibonacciHeap *fh1, FibonacciHeap *fh2){
        if(!fh2->_minimum){
            // delete fh2;
            return fh1;
        }
        if(!fh1->_minimum){
            // delete fh1
            return fh2;
        }
        fh1->_minimum->_prev->_next = fh2->_minimum->_next;
        fh2->_minimum->_next->_prev = fh1->_minimum->_prev;
        fh2->_minimum->_next = fh1->_minimum;
        fh1->_minimum->_prev = fh2->_minimum;
        // delete fh2;
        return fh1;
    }
};

template<typename _Tp>
class Arborescence {
private:
    struct edge {
        const int from, to;
        const _Tp cost;
        edge(int _from, int _to, _Tp _cost)
            : from(_from), to(_to), cost(_cost){}
    };
    struct cycle_edge {
        int from, super_from, super_to;
        const edge *e;
        cycle_edge(int _from, int _super_from, int _super_to, const edge *_e)
            : from(_from), super_from(_super_from), super_to(_super_to), e(_e){}
    };
    const int V;
    int super_id;
    vector<vector<edge> > revG;
    vector<forward_list<cycle_edge> > cycle;
    forward_list<forward_list<cycle_edge> > path;
    vector<forward_list<const edge*> > passive;
    vector<int> used, heap, par;
    LazyUnionFind<_Tp> uf;
    vector<FibonacciHeap<edge, _Tp>*> fh;
    vector<typename FibonacciHeap<edge, _Tp>::node*> nodes;
    void _move_node(int prev, int vertex, int next, const edge *e){
        move_node(fh[prev], nodes[vertex], fh[next]);
        heap[vertex] = next, nodes[vertex]->_data = e;
    }
    void grow_path(int cur, forward_list<int>& visit){
        visit.push_front(cur);
        fh[cur] = new FibonacciHeap<edge, _Tp>(uf, fh, heap);
        for(const edge& e : revG[cur]){
            const int from = uf.find(e.from);
            if(cur == from) continue;
            if(heap[from] < 0){
                heap[from] = cur, nodes[from] = fh[cur]->push(&e);
            }else if(heap[from] == cur){
                if(e.cost < nodes[from]->_data->cost){
                    nodes[from]->_data = &e;
                    fh[cur]->to_root(nodes[from]);
                }
            }else{
                const int hfrom = uf.find(heap[from]);
                passive[hfrom].emplace_front(nodes[from]->_data);
                _move_node(hfrom, from, cur, &e);
            }
        }
    }
    int contract_cycle(int cur, forward_list<cycle_edge>&& cur_path){
        int modify = -1;
        for(auto it = next(cur_path.begin(), 1);; ++it){
            const int v = (it == cur_path.end()) ? cur : it->from;
            while(!passive[v].empty()){
                const edge *e = passive[v].front();
                const int from = uf.find(e->from);
                if(nodes[from] &&
                    nodes[from]->_data->cost + uf.find_value(nodes[from]->_data->to) >
                        e->cost + uf.find_value(e->to)){
                    _move_node(uf.find(heap[from]), from, v, e);
                }
                passive[v].pop_front();
            }
            if(v == cur) break;
            par[it->super_from] = super_id;
            modify = it->super_to;
        }
        par[cur_path.begin()->super_from = modify] = super_id;
        int ver = -1;
        for(auto it = next(cur_path.begin(), 1);; ++it){
            const int v = (it == cur_path.end()) ? cur : it->from;
            if(heap[v] >= 0) fh[v]->delete_node(nodes[v]), heap[v] = -1, nodes[v] = nullptr;
            if(ver >= 0){
                if(uf.unite(ver, v) == ver) fh[ver] = meld(fh[ver], fh[v]), fh[v] = nullptr;
                else fh[v] = meld(fh[ver], fh[v]), fh[ver] = nullptr, ver = v;
            }else ver = v;
            if(v == cur){
                cycle[super_id].push_front(cur_path.front()), cur_path.pop_front();
                cycle[super_id].splice_after(cycle[super_id].begin(), move(cur_path),
                    cur_path.before_begin(), it);
                break;
            }
        }
        if(!cur_path.empty()) cur_path.begin()->from = ver, cur_path.begin()->super_from = super_id;
        return ver;
    }
    void cycle_dfs(int u, int par_cycle){
        while(u != par[u] && par[u] != par_cycle){
            const int p = par[u];
            for(auto it = cycle[p].begin(); it != cycle[p].end(); ++it){
                if(u == it->super_to) continue;
                const int w = it->e->to;
                parent[w] = it->e->from;
                if(p != par[w]) cycle_dfs(w, p);
            }
            u = p;
        }
    }
    void identify_tree(){
        for(forward_list<cycle_edge>& cur_path : path){
            while(!cur_path.empty()){
                int v = cur_path.begin()->e->to;
                parent[v] = cur_path.begin()->e->from;
                cur_path.pop_front();
                cycle_dfs(v, -1);
            }
        }
    }
public:
    _Tp ans;
    vector<int> parent;
    Arborescence(int node_size)
        : V(node_size), super_id(V), revG(V), cycle(2*V), passive(V),
            used(V, -1), heap(V, -1), par(2*V), uf(V), fh(V, nullptr), nodes(V, nullptr),
                ans((_Tp)0), parent(V, -1){
        iota(par.begin(), par.end(), 0);
    }
    // ~Arborescence(){
        // for(int i = 0; i < V; ++i) if(nodes[i]) delete nodes[i];
        // for(int i = 0; i < V; ++i) if(fh[i]) delete fh[i];
    // }
    void add_edge(int from, int to, _Tp cost){
        revG[to].emplace_back(from, to, cost);
    }
    _Tp solve(const int root){
        used[root] = 1;
        for(int i = 0; i < V; ++i){
            if(used[i] != -1) continue;
            int cur = i, super_cur;
            forward_list<cycle_edge> cur_path;
            forward_list<int> visit;
            while(used[cur] != 1){
                if(used[cur] == -1){
                    grow_path(cur, visit);
                    used[super_cur = cur] = 0;
                }else{
                    cur = contract_cycle(cur, move(cur_path));
                    super_cur = super_id++;
                }
                const edge *e = fh[cur]->pop();
                if(!e) return numeric_limits<_Tp>::max();
                ans += e->cost + uf.find_value(e->to);
                uf.change_value(-e->cost - uf.find_value(e->to), cur);
                const int next = uf.find(e->from);
                heap[next] = -1, nodes[next] = nullptr;
                cur_path.emplace_front(next, e->from, super_cur, e);
                cur = next;
            }
            path.push_front(move(cur_path));
            for(const int ver : visit) used[ver] = 1;
        }
        identify_tree();
        return ans;
    }
};

verify 用の問題

AOJ : Minimum-Cost Arborescence 提出コード
yosupo さんの library checker : Directed MST 提出コード