Caml-list

A couple weeks ago, Kim Quyen Ly asked a question about algorithms for strongly connected components.

I answered the best known algorithms were Mehlhorn-Gabow's and Tarjan's, both linear in number of arcs.

I wrote an implementation in Caml but was unsatisfied with it because in theory these algorithms use 3 data structures in total (2 arrays + 1 stack or 1 array + 2 stacks).

However I had an extra stack because of the recursion, and couldn't figure out how to merge the call-stack with the open-node stack of the dfs.

I finally checked in the original Mehlhorn paper (Algorithmica 1996) and Sedgewick implementations, to find out that not only did they use a recursive function as well but they had MORE arrays and stacks than theoretically required. I guess I will have to wait Knuth reaches the corresponding TAOCP volume to uncompile his assembler code and finally know the truth.

So here are my implementations. The example graph was built to show the case where the call-stack and open-nodes stack differ (node 4).

I recommend potential users to prove correction before using the implementation, I am a lousy coder.

(* Make matrix from list *)

let to_matrix = function list ->

  let n = 1 + List.fold_left (fun current (i, j) -> max current (max i j)) 0 list in

  let matrix = Array.make_matrix n n 0 in

  let rec add_arc = function

    | [] -> matrix

    | (i, j) :: tail -> matrix.(i).(j) <- 1; add_arc tail

  in add_arc list

(* Example built to show the open-node stack / dfs call-stack difference *)

let example = to_matrix [(0, 1); (1, 2); (2, 3); (3, 4); (4, 2); (2, 1); (3, 5); (5, 6); (6, 5)]

(* Mehlhorn Gabow scc *)

let cmg_scc = function matrix ->

  let n = Array.length matrix in

  let 

      visited_at_depth = Array.make n max_int and

      roots = Stack.create () and

      open_nodes = Stack.create ()

  in

  let rec unstack_until = function i ->

    match Stack.pop open_nodes with

      | n when n = i -> [i]

      | n -> n :: unstack_until i

  in

  let rec dfs depth = function i ->

    let result = ref [] in

    

    (* mark *)

    Stack.push depth roots;

    Stack.push i open_nodes;

    visited_at_depth.(i) <- depth;

    

    (* dive *)

    for j = 0 to n - 1 do

      if (matrix.(i).(j) = 1) && (visited_at_depth.(j) = max_int) then 

result := dfs (depth + 1) j @ !result

    done;

    (* process reverse-arcs *)

    for j = 0 to n - 1 do

      if (matrix.(i).(j) = 1) && (visited_at_depth.(j) < depth) then 

let scc_returns_to_depth = visited_at_depth.(j) in

while Stack.top roots > scc_returns_to_depth do ignore (Stack.pop roots) done

    done;    

    

    (* emit connected component if current node is root *)

    if depth = Stack.top roots then

      (

ignore (Stack.pop roots);

unstack_until i :: !result

      )

    else

      !result

  in

  let result = ref [] in

  for i = 0 to n - 1 do 

    if (visited_at_depth.(i) = max_int) then result := (dfs 0 i) @ !result

  done;

  !result

     

(* Tarjan scc *)

let tarjan_scc = function matrix ->

  

  let n = Array.length matrix in

  let 

      visited_at_depth = Array.make n max_int and

      scc_root = Array.make n max_int and

      open_nodes = Stack.create () and

      result = ref []

  in

  

  let rec unstack_until = function i ->

    match Stack.pop open_nodes with

      | n when n = i -> [i]

      | n -> n :: unstack_until i

  in

  let rec dfs depth = function i ->

    

    (* mark *)

    visited_at_depth.(i) <- depth;

    scc_root.(i) <- depth;

    Stack.push i open_nodes;

    for j = 0 to n - 1 do

      if matrix.(i).(j) = 1 then

if visited_at_depth.(j) = max_int then

 scc_root.(i) <- min scc_root.(i) (dfs (depth + 1) j) (* dive *)  

else 

 scc_root.(i) <- min scc_root.(i) visited_at_depth.(j) (* reverse-arc *)

    done;

    

    (* emit connected component if current node is root *)

    if scc_root.(i) = visited_at_depth.(i) then 

      result := unstack_until i :: !result;

    scc_root.(i)

  in

  

  for i = 0 to n - 1 do 

    if (visited_at_depth.(i) = max_int) then ignore (dfs 0 i)

  done;

  !result

         Diego Olivier