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