$\newcommand{\O}{\mathrm{O}}$

有向グラフと点 $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;
}
};
AOJ : Minimum-Cost Arborescence
提出コード
yosupo さんの library checker : Directed MST
提出コード