with(ArrayTools):

# Helper method to filter a set S using f : S -> {true,false}
filter := proc(f, S) map(x -> if f(x) then x else  end if, S) end proc:

# Inputs a digraph G. Outputs the adjacency list of G having length |V|.
toAdjList := proc(G)
local adj, v, e;
    adj := Array(1 .. 0);
    for v to G[1] do Append(adj, Array(1 .. 0), inplace = true) end do;
    for e in G[2] do Append(adj[e[1]], e[2], inplace = true) end do;
    adj
end proc:

# Inputs a rooted digraph G := [V, E, S] and performs depth-first search
# starting from the root, S. Outputs an array having length |V| whose ith entry
# is 1 if vertex i was reached, and 0 otherwise.
DFS := proc(G)
local adj, rec, visited;
    adj := toAdjList(G);
    rec := proc(v)
        local w;
            visited[v] := 1;
            for w in adj[v] do if visited[w] = 0 then rec(w) end if end do;
            visited
        end proc;
    visited := Array(1 .. G[1], fill = 0);
    rec(G[3])
end proc:

# Inputs a rooted digraph G := [V, E, S], where every vertex is be reachable
# from S. Outputs an array idom where idom[i] is the immediate dominator of
# vertex i, i.e. the parent of the ith node in the dominator tree. Note that
# idom[S] is not well-defined since S is the root, so idom[S] = 0.
constructDominator := proc(G)
local n, E, r, V, domBy, v, S, i, dom, w, idom, u;
    n := G[1];
    E := G[2];
    r := G[3];
    V := DFS(G);
    domBy := Array(1 .. n);
    for v to n do
        S := DFS([n, filter(e -> not v in e, E), r]);
        S := [seq(V[i] - S[i], i = 1 .. n)];
        domBy[v] := filter(i -> S[i] = 1, {seq(i, i = 1 .. n)}) minus {v}
    end do;
    dom := Array(1 .. n, fill = {});
    for v to n do for w in domBy[v] do dom[w] := dom[w] union {v} end do end do;
    idom := Array(1 .. n, fill = 0);
    for w to n do
        for u in dom[w] while idom[w] = 0 do
            idom[w] := u;
            for v in dom[w] do
                if v <> u and not v in dom[u] then idom[w] := 0 end if
            end do
        end do
    end do;
    idom
    # [n, {seq([idom(v), v], v = 1 .. n)} minus {[0, r]}, r]
end proc:

# Clever version of constructDominator by Lengauer and Tarjan (1979).
# Inputs a rooted digraph G := [V, E, S], where every vertex is be reachable
# from S. Outputs an array idom where idom[i] is the immediate dominator of
# vertex i, i.e. the parent of the ith node in the dominator tree. Note that
# idom[S] is not well-defined since S is the root, so idom[S] = 0.
constructDominatorC := proc(G)
local n, E, S, adj, e, u, v, i, i1, j, T, w, radj, tin, rtin, dfs, deg, rdeg,
treePar, par, label, sdom, merge, find, bucket, idom, idom_t;
    n := G[1];
    E := G[2];
    S := G[3];
    adj := toAdjList(G);
    radj := Array(1 .. n);
    for i to n do radj[i] := Array(1 .. 0) end do;
    tin := Array(1 .. n, fill = 0);
    rtin := Array(1 .. n, fill = 0);
    treePar := Array(1 .. n, fill = -1);
    T := 1;
    dfs := proc(u)
        local v;
            tin[u] := T;
            rtin[T] := u;
            T := T + 1;
            for v in adj[u] do
                if tin[v] = 0 then dfs(v); treePar[tin[v]] := tin[u] end if;
                Append(radj[tin[v]], tin[u], inplace = true)
            end do
        end proc;
    dfs(S);
    for i to n do
        if tin[i] = 0 then return "Graph is not entirely reachable from S"
        end if
    end do;
    par := Array(1 .. n, fill = 0);
    label := Array(1 .. n, fill = 0);
    sdom := Array(1 .. n, fill = 0);
    for i to n do par[i] := i; label[i] := i; sdom[i] := i end do;
    merge := proc(u, v) par[v] := u end proc;
    find := proc(u, x::integer := 0)
        local v;
            if u = par[u] then if x = 0 then return u end if; return -1
            end if;
            v := find(par[u], x + 1);
            if v = -1 then return u end if;
            if sdom[label[par[u]]] < sdom[label[u]] then
                label[u] := label[par[u]]
            end if;
            par[u] := v;
            if x = 0 then return label[u] end if;
            return v
        end proc;
    idom := Array(1 .. n, fill = 0);
    bucket := Array(1 .. n, fill = 0);
    for i to n do bucket[i] := Array(1 .. 0) end do;
    for i1 to n do
        i := n - i1 + 1;
        for v in radj[i] do
            u := find(v); if sdom[u] < sdom[i] then sdom[i] := sdom[u] end if
        end do;
        if 1 < i then Append(bucket[sdom[i]], i, inplace = true) end if;
        for w in bucket[i] do
            v := find(w);
            if sdom[v] = sdom[w] then idom[w] := sdom[w]
            else idom[w] := v
            end if
        end do;
        if 1 < i then merge(treePar[i], i) end if
    end do;
    idom_t := Array(1 .. n, fill = 0);
    idom_t[S] := 0;
    for i from 2 to n do
        if idom[i] <> sdom[i] then idom[i] := idom[idom[i]] end if
    end do;
    for i from 2 to n do idom_t[rtin[i]] := rtin[idom[i]] end do;
    return idom_t
end proc:

# Generate a random connected digraph with n vertices.
genRandomDigraph := proc(n)
local i, j, G, E, S, ra;
    G := [n, {}, 1];
    ra := rand(1 .. 2);
    for i to n do
        for j to n do
            if i <> j and ra() = 1 then G[2] := G[2] union {[i, j]} end if
        end do
    end do;
    for i from 2 to n do
        ra := rand(1 .. i - 1); G[2] := G[2] union {[ra(), i]}
    end do;
    G
end proc:

# Perform K trials with random digraphs of size n.
# Check if constructDominator and constructDominatorC give the same result.
# If not, print the graph where the error occured and the 2 different results.
checkClever := proc(a, b, n)
local i;
    for i to n do if a[i] <> b[i] then return false end if end do;
    return true
end proc:

# Check that constructDominator and constructDominatorC agree on K random
# digraphs on n vertices.
# If not, print the graph where the error occured and the 2 different results.
stressTest := proc(n, K)
local i, G, idom1, idom2;
    for i to K do
        G := genRandomDigraph(n);
        idom1 := constructDominatorC(G);
        idom2 := constructDominatorC(G);
        if not checkClever(idom1, idom2, n) then
            print("Error in graph: ", G);
            print("Brute Force Output: ", idom1);
            print("Fast Output: ", idom2);
            return
        end if
    end do;
    print("All tests passed")
end proc:

# Fig. 1 from paper (R=1, A=2, B=3, etc.)
Fig1 := [13,{[1,2],[1,3],[1,4],[2,5],[3,2],[3,5],[3,6],[4,7],[4,8],[5,13],[6,9],[7,10],[8,10],[8,11],[9,6],[9,12],[10,12],[11,10],[12,10],[12,1],[13,9]},1]: