################################################################################
#
#  kPerl Logics/Mathematics Library Group
#  Library for Analysis Methods
#
#  copyright (c)
#    Karol Szafranski, 2005
#  author
#    Karol Szafranski, karol.szafranski@leibniz-fli.de
#
################################################################################
#
#  FUNCTIONS, DATA
#
#   @EXPORT
#
# - analytic geometry
#   &determinant
#
#
#  STD OPTIONS
#
#  -debug      print debug protocol to STDERR
#
################################################################################

package Math::Analysis;

# includes
#use strict; use warnings;  # OK 20050623

# symbol export
use Exporter qw(import);
our @EXPORT = qw (
  &determinant
  );


################################################################################
# analytic geometry
################################################################################


# random integer
#
# INTERFACE
# - argument 1: reference to 2D matrix == nested array
#               1st dim. := rows
#               2nd dim. := columns
# - options:
#   -debug      [STD]
#   -init       constructor value for determinant result. This is useful if
#               the matrix contains objects
#   -ncol       number of columns, default int(@{$pMat->[0]})
#   -nrow       number of rows, default int(@$pMat)
# - return val: matrix determinant
#
# DEBUG, CHANGES, ADDITIONS
# - There is no linear formula to compute the determinant beyond 3D space.
#
sub determinant {
  my ($pMat,%opt) = @_;
  my $debug = $opt{-debug};

  # establish working copy of the argument matrix
  my $nrow = $opt{-nrow} || int(@$pMat);
  my $ncol = $opt{-ncol} || int(@{$pMat->[0]});
  my $pWork = [ @$pMat ];
  for (my $i=0; $i<$nrow; ++$i) {
    $pWork->[$i] = [ @{$pMat->[$i]} ];

    # "circularize" matrix by addition of columns
    for (my $j=0; $j<$ncol-1; ++$j) {
      $pWork->[$i][$ncol+$j] = $pWork->[$i][$j] ||= 0;
    }
  }

  # calculate
  # I do not use mutating operators +=,*= since this would not work with over-
  # loaded operators in expressions ($scalar += $pObj)
  my @diagon;
  my $discr = $opt{-init} || 0;
  for (my $i=0; $i<$ncol; ++$i) {
    my $plus=1;
    for (my $j=0; $j<$nrow; ++$j) {
      $plus = $plus * $pWork->[$j][$i+$j];
    }
    $debug and push @diagon,$plus;
    $discr = $discr + $plus;
  }
  for (my $i=0; $i<$ncol; ++$i) {
    my $minus=1;
    for (my $j=0; $j<$nrow; ++$j) {
      $minus = $minus * $pWork->[$nrow-1-$j][$i+$j];
    }
    $debug and push @diagon,-$minus;
    $discr = $discr - $minus;
  }

  if ($debug) {
    require MainLib::Data;
    printf STDERR "%s. array of diagonals:\n", (caller(0))[3];
    &MainLib::Data::DataPrint(\@diagon,-handle=>\*STDERR,-space=>2,-NoAddr=>1);
  }
  return $discr;
}


1;
# $Id: Analysis.pm,v 1.6 2018/06/05 18:02:56 szafrans Exp $
