############################################################################## # # # matrix.pm - Perl module for doing simple matrix math. # # # # Revision: 1.0 1/20/98 # # # # Author: David Harris # # # # History: None yet. # # # # # # COPYRIGHT (C) 1998 David R. Harris All Rights Reserved. # # # ############################################################################## package matrix; ############################################################################## ############################################################################## # DOCUMENTATION # # This module is designed to do matrix mathematics in perl. The matrices can # be of any scalar value that has addition and multiplication operations # supported, namely integers. You could also used an overloaded something. # # NOTE: What a methematician would call a 3x3 matrix, this module will call # a matrix with height = 2 and width = 2, because we programmser begin # counting at zero. BEGIN: { # # Package globals # %multiply_cache; } ############################################################################## # SUB: new($ref_to_list) # Takes a reference to a list of references to lists and returns a matrix # object. sub new { my ($ref_to_list) = @_; return { width => $#{$ref_to_list->[0]}, height => $#$ref_to_list, data => flatten($ref_to_list), } } ############################################################################## # SUB: identity($width) # Returns an identity matrix with width $width. sub identity { my ($width) = @_; $data = []; foreach $shift ( 0 .. ($width) ) { my @this_row = (0) x ( $width + 1 ); $this_row[$shift] = 1; push(@$data, @this_row); } return { width => $width, height => $width, data => $data, } } ############################################################################## # SUB: copy_matrix($m1) # Makes a copy of a matrix. sub copy_matrix { my ($m1) = @_; return { width => $m1->{width}, height => $m1->{height}, data => [ @{$m1->{data}} ], } } ############################################################################## # SUB: multiply($m1, $m2) # Multiplies two matrixes. Returns a matrix with the height of the first # argument and the width of the second. Optimized to do repedative # multiplication of matricies of the same size. sub multiply { my ($m1, $m2) = @_; # # Make sure we have the correct dimensions # die("[matrix.pm @{[__LINE__]}][Incompatable dimensions]") if ( $m1->{width} != $m2->{height} ); # # Do the multiply # my $ref_to_sub = get_multiply_sub( $m1->{width}, $m1->{height}, $m2->{width}, $m2->{height} ); return { width => $m2->{width}, height => $m1->{height}, data => &$ref_to_sub($m1->{data}, $m2->{data}), } } ############################################################################## # SUB: add_into_scalar_multiply($m1, $m2, $s) # Adds to one matrix into another. Same is $m1 += $m2 * $s; sub add_into_scalar_multiply { my ($m1, $m2, $s) = @_; die("[matrix.pm @{[__LINE__]}][Incompatable dimensions]") if ( $m1->{width} != $m2->{width} || $m1->{height} != $m2->{height} ); # # Do the addition # my $loop_to = ( $m1->{width} + 1 ) * ( $m1->{height} + 1 ) - 1; my $c; my $m1data = $m1->{data}; my $m2data = $m2->{data}; for ( $c = 0 ; $c <= $loop_to ; $c++ ) { $m1data->[$c] += $m2data->[$c] * $s; } } ############################################################################## # SUB: matrix_exponential_use_precompute($precompute, $s, $order) # Computes the matrix exponential exp( $m * $s ) where $m is a matrix and $s # is a scalar. The order of the internal taylor serise is given by $order, # therefor the number of terms is $order + 1. Get $precompute from the # fuction matrix_exponential_precompute. sub matrix_exponential_use_precompute { my ($precompute, $s, $order) = @_; # # Loop for each term in the tylor serise creating the sum # my $k; my $results_hash = {}; my $this_multiplier = 1; my $running_sum = identity($precompute->{1}{width}); for ( $k = 1 ; $k <= $order ; $k++ ) { $this_multiplier *= $s / $k; add_into_scalar_multiply($running_sum, $precompute->{$k}, $this_multiplier); $results_hash->{$k} = copy_matrix($running_sum); } # # Return the answer # return $results_hash; } ############################################################################## # SUB: matrix_exponential_do_precompute($m, $order) # Returns the precomputation requred for matrix_exponential_use_precompute. sub matrix_exponential_do_precompute { my ($m, $order) = @_; # # Loop for each term in the tylor serise creating the sum # my $k; my $m_tothe_k = identity($m->{width}); my $return_value; for ( $k = 1 ; $k <= $order ; $k++ ) { $m_tothe_k = multiply($m_tothe_k, $m); $precompute->{$k} = $m_tothe_k; } # # Return the answer # return $precompute; } ############################################################################## ############################################################################## ############################################################################## # P R I V A T E S U B S ############################################################################## # PRIVATE SUB: get_multiply_sub($width1, $height1, $width2, $height2) # Returns a reference to a sub wich will multiply two matricies with the # sizes given. We need $width1 == $height2 for this to be valid and # return a matrix $width2 by $height1 sub get_multiply_sub { my ($width1, $height1, $width2, $height2) = @_; # # See if this request is valid # $main::debug_string .= "\$width1 = $width1\n"; $main::debug_string .= "\$height1 = $height1\n"; $main::debug_string .= "\$width2 = $width2\n"; $main::debug_string .= "\$height2 = $height2\n"; main::my_die("[matrix.pm @{[__LINE__]}][Invalid multiply]") if ($width1 != $height2); # # See if we have this cached. If so, return it. # my $key_value = join('_', $width1, $height1, $width2, $height2); return $multiply_cache{$key_value} if ( defined $multiply_cache{$key_value} ); # # Generate a the code of the ref to sub # my ($i, $k, $j); my $code = 'return sub { my ($m1, $m2) = @_; return [' . "\n"; for ( $k = 0 ; $k <= $height1 ; $k++ ) { for( $i = 0 ; $i <= $width2 ; $i++ ) { for ( $j = 0 ; $j <= $width1 ; $j++ ) { $code .= " \$m1->[@{[ $j + ( $width1 + 1 ) * $k ]}] * \$m2->[@{[ $i + ( $width2 + 1 ) * $j ]}] +" } $code .= " 0,\n"; } } $code .= " ]; }"; # # Store and return the ref to sub # $main::debug_string .= "\$code = <<'EOT'\n$code\nEOT\n"; $multiply_cache{$key_value} = eval($code); die("[matrix.pm @{[__LINE__]}][Error in eval][\$\@ = $@]") if ($@ ne ''); return $multiply_cache{$key_value}; } ############################################################################## # PRIVATE SUB: flatten($ref_to_list) # Flattens a ref to a list to be only one level deep. Fore example, the list # [ 1, [2, 3] would become [1, 2, 3]. Only does one level of flattening. sub flatten { my ($ref_to_list) = @_; my $return_list = []; foreach $this_element ( @$ref_to_list ) { if ( ref($this_element) ne 'ARRAY' ) { push(@$return_list, $this_elemement); } else { push(@$return_list, @$this_element); } } return $return_list; } 1; # Required for packages ############################################################################## ############################################################################## ############################################################################## # U N U S E D S U B S __END__ ############################################################################## # SUB: scalar_multiply($scalar, $m1) # Multiplies a matrix by a scalar and returns a matrix. sub scalar_multiply { my ($m1) = @_; return { width => $m1->{width}, height => $m1->{height}, data => [ @{$m1->{data}} * $scalar ], } } ############################################################################## # SUB: add($m1, $m2) # Adds to matrices together and returns a new one. sub add { my ($m1, $m2) = @_; die("[matrix.pm @{[__LINE__]}][Incompatable dimensions]") if ( $m1->{width} =! $m2->{width} || $m1->{height} =! $m2->{height}} ); # # Generate new data # my $loop_to = ( $m1->{width} + 1 ) * ( $m1->{height} + 1 ) - 1; my $data = []; $#$data = $loop_to; my $c; my $m1data = $m1->{data}; my $m2data = $m2->{data}; for ( $c = 0 ; $c <= $loop_to ; $c++ ) { $data[$c] = $m1data[$c] + $m2data[$c]; } # # Make and return a new matrix # return { width => $m1->{width}, height => $m1->{height}, data => [ @{$m1->{data}} * $scalar ], } } ############################################################################## # SUB: add_into($m1, $m2) # Adds to one matrix into another. Same is $m1 += $m2; sub add_into { my ($m1, $m2) = @_; die("[matrix.pm @{[__LINE__]}][Incompatable dimensions]") if ( $m1->{width} =! $m2->{width} || $m1->{height} =! $m2->{height}} ); # # Do the addition # my $loop_to = ( $m1->{width} + 1 ) * ( $m1->{height} + 1 ) - 1; my $c; my $m1data = $m1->{data}; my $m2data = $m2->{data}; for ( $c = 0 ; $c <= $loop_to ; $c++ ) { $m1data[$c] += $m2data[$c]; } } ############################################################################## # SUB: empty_matrix($width, $height, $scalar) # Returns a matrix $width by $height filled with the value $scalar. sub empty_matrix { my ($width, $height, $scalar) = @_; return { width => $width, height => $height, data => [ ($scalar) x ( ( $width + 1 ) * ( $height + 1 ) - 1 ) ], } } ############################################################################## # SUB: matrix_exponential($m, $s, $order) # Computes the matrix exponential exp( $m * $s ) where $m is a matrix and $s # is a scalar. The order of the internal taylor serise is given by $order, # therefor the number of terms is $order + 1. sub matrix_exponential { my ($m, $s, $order) = @_; # # Loop for each term in the tylor serise creating the sum # my $k; my $m_tothe_k = identity_matrix($m->{width}); my $this_multiplier = 1; my $running_sum = identity($m->{width}); for ( $k = 1 ; $k <= $order ; $k++ ) { $m_tothe_k = multiply($m_tothe_k, $m); $this_multiplier *= $s / $k; add_into_scalar_multiply($running_sum, $m_tothe_k, $this_multiplier); } # # Return the answer # return $running_sum; }