use strict;

## just any old directed graph...

my $adj = [ [0,1,1,0,1],
            [1,0,0,0,1],
            [0,1,0,0,0],
            [1,1,1,0,0],
            [0,0,1,1,0] ];
my $dim = 4; ## size of the graph
my $N = 3;   ## length of paths (# of edges, not # of vertices).

## convert each entry in the adjacency matrix into a list of paths
    
for my $i (0 .. $dim) {
  for my $j (0 .. $dim) {
    $adj->[$i][$j] = $adj->[$i][$j] ? [$j] : [];
  }   
}

## compute the $N-th power of the adjacency matrix with our modified
## multiplication

my $result = $adj;
for (2 .. $N) {
   print_paths($result);
   print "========\n";
   $result = matrix_mult($result, $adj);
}

print_paths($result);

## the i,j entry of the matrix is a list of all the paths from i to j, but
## without "i," at the beginning, so we must add it

sub print_paths {
    my $M = shift;
    my @paths;
    for my $i (0 .. $dim) {
      for my $j (0 .. $dim) {
        push @paths, map { "$i,$_" } @{ $M->[$i][$j] };
      }
    }
    print map "$_\n", sort @paths;
}

## modified matrix multiplication. instead of multiplication, we
## combine paths from i->k and k->j to get paths from i->j (this is why
## we include the endpoint in the path, but not the starting point).
## then instead of addition, we union all these paths from i->j

sub matrix_mult {
    my ($A, $B) = @_;
    my $result;
    for my $i (0 .. $dim) {
      for my $j (0 .. $dim) {
        my @result;
        for my $k (0 .. $dim) {
          push @result, combine_paths( $A->[$i][$k], $B->[$k][$j] );
        }
        $result->[$i][$j] = \@result;
      }
    }
    $result;
}

## the sub to combine i->k paths and k->j paths into i->j paths --
## we simply concatenate with a comma in between.

sub combine_paths {
    my ($x, $y) = @_;
    my @result; 
    
    for my $i (@$x) {
      for my $j (@$y) {
        push @result, "$i,$j";
      } 
    }   
    @result;
}