#!/usr/bin/perl -w

use PDL;  # http://www.perl.com/CPAN/modules/by-module/PDL

# Create an array of five matrices.
@matrices = (pdl ([[1,2],[3,4],[5,6],[7,8],[9,10],[11,12],[13,14]]),
             pdl ([[1,2,3],[4,5,6]]),
             pdl ([[1,2,3],[4,5,6],[7,8,9]]),
             pdl ([[1,2,3,4,5,6],[7,8,9,10,11,12],[13,14,15,16,17,18]]),
             pdl ([[1,2,3],[4,5,6],[7,8,9],[10,11,12],[13,14,15],
                  [16,17,18]]));

# Initialize the three auxiliary matrices that we'll use to
# store the costs (number of scalar multiplications),
# the parenthesization so far, and the dimensions of what the
# intermediate product would be if we were to compute it.

for ($i = 0; $i < @matrices; $i++) {
    $costs[$i][$i]  = 0;
    $parens[$i][$i] = '$matrices[' . $i . ']';
    $dims[$i][$i]   = [dims $matrices[$i]];
}

# Determine the costs of the pairs ($i == 1), then the triples
# ($i == 2), the quadruples, and finally all five matrices.

for ($i = 1; $i < @matrices; $i++) {

    # Loop through all of the entries on each diagonal.
    #
    for ($j = $i; $j < @matrices; $j++) { # column

        # Determine the best parenthesization for the entry
        # at row $j-$i and column $j.
        #
        for ($k = $j - $i; $k < $j; $k++) {
            ($col1, $row1) = @{$dims[$j-$i][$k]};
            ($col2, undef) = @{$dims[$k+1][$j]};

            # Compute the cost of this parenthesization.
            #
            $try = $costs[$j-$i][$k] + $costs[$k+1][$j] +
                       $row1 * $col1 * $col2;

            # If it's the lowest we've seen (or the first we've seen),
            # store the cost, the dimensions, and the parenthesization.
            #
            if (!defined $costs[$j-$i][$j] or $try < $costs[$j-$i][$j]) {
                $costs[$j-$i][$j] = $try;
                $dims[$j-$i][$j] = [$col2, $row1];
                $parens[$j-$i][$j] = "(" . $parens[$j-$i][$k] . "x" .
                    $parens[$k+1][$j] . ")";
            }
        }
    }
}

# At this point, all of the information we need has been propagated
# to the upper right corner of our master matrix: the parenthesizations
# and the number of scalar multiplications.

print "Evaluating:\n", $parens[0][$#matrices], "\n";
print "\tfor a total of $costs[0][$#matrices] scalar multiplications.\n";

# Evaluate the string and, finally, multiply our matrices!
print eval $parens[0][$#matrices];

__END__
Evaluating:
($matrices[0]x(($matrices[1]x$matrices[2])x($matrices[3]x$matrices[4])))
        for a total of 132 scalar multiplications.
[
 [ 341010  377460  413910]
 [ 743688  823176  902664]
 [1146366 1268892 1391418]
 [1549044 1714608 1880172]
 [1951722 2160324 2368926]
 [2354400 2606040 2857680]
 [2757078 3051756 3346434]
]