##############################################################################
#                                                                            #
#  matrix.pm - Perl module for doing simple matrix math.                     #
#                                                                            #
#  Revision:  1.0  1/20/98                                                   #
#                                                                            #
#  Author:    David Harris     <dharris@drh.net>                             #
#                                                                            #
#  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;
}