wolffd@0: function C = reachability_graph(G)
wolffd@0: % REACHABILITY_GRAPH C(i,j) = 1 iff there is a path from i to j in DAG G
wolffd@0: % C = reachability_graph(G)
wolffd@0: 
wolffd@0: if 1
wolffd@0:   % expm(G) = I + G + G^2 / 2! + G^3 / 3! + ...
wolffd@0:   M = expm(double(full(G))) - eye(length(G));
wolffd@0:   C = (M>0);
wolffd@0: else
wolffd@0:   % This computes C = G + G^2 + ... + G^{n-1}
wolffd@0:   n = length(G);
wolffd@0:   A = G;
wolffd@0:   C = zeros(n);
wolffd@0:   for i=1:n-1
wolffd@0:     C = C + A;
wolffd@0:     A = A * G;
wolffd@0:   end
wolffd@0:   C = (C > 0);
wolffd@0: end