File Coverage

blib/lib/Math/Assistant.pm

Criterion	Covered	Total	%
statement	184	195	94.3
branch	59	82	71.9
condition	14	23	60.8
subroutine	12	12	100.0
pod	5	5	100.0
total	274	317	86.4

package Math::Assistant;

use 5.008004;
use strict;
use warnings;

use Carp;
use Math::BigInt qw( bgcd );

require Math::BigFloat;
require Exporter;

our @ISA = qw(Exporter);

our %EXPORT_TAGS = ( 'all' => [ qw(
			Rank Det Solve_Det Gaussian_elimination test_matrix
			) ],
		    'algebra' => [ qw( Rank Det Solve_Det ) ],
);

our @EXPORT_OK = ( @{ $EXPORT_TAGS{'all'} } );

our @EXPORT = qw( );

our $VERSION = '0.05';

use base qw(Exporter);


# Rank of any integer matrix
# input:
#	$Dimention
# return:
#	$rank (order)
sub Rank{
    my $M = &Gaussian_elimination; # triangular matrix
    my $rank = 0;

    for( @$M ){ # rows
	for( @{$_} ){ # element in row
	    $rank++, last if $_
	}
    }
    $rank;
}


# Method of Gaussian elimination
# input:
#	$Dimention
# return:
#	$elimination_Dimention
sub Gaussian_elimination{
    my $M_ini = shift;

    my $t = &test_matrix( $M_ini );
    if( $t > 3 ){
	croak("Use of uninitialized value in matrix");

    }elsif( $t > 1 ){
	croak("Bad matrix");
    }

    my $rows = $#{$M_ini}; # number of rows
    my $cols = $#{$M_ini->[0]}; # number of columns

    my %zr; # for rows with 0
    my %zc; # for columns with 0

    my $M; # copy
    # search of rows with 0
    for( my $i = 0; $i <= $rows; $i++ ){
	for( my $j = 0; $j <= $cols; $j++ ){
	    $M->[$i][$j] = $M_ini->[$i][$j];
	    unless($M->[$i][$j]){ # zero element
		$zr{$i}++;
		$zc{$j}++;
	    }
	}
    }

    # Check Float (Fraction)
    for my $i ( @$M ){
	my $max_frac = &_max_len_frac($i, 0);

	if($max_frac){
	    $_ *= 10**$max_frac for @$i;
	}
    }

    if(keys %zr){ # yes, zero rows (columns)

	# I supplement indexes of nonzero rows
	for( my $i = 0; $i <= $rows; $i++ ){
	    $zr{$i} += 0;
	}

	my $R; # temp copy of matrix M
	my $v = 0;
	# replacement of rows
	for my $i ( sort { $zr{$a} <=> $zr{$b} } keys %zr ){
	    for( my $j = 0; $j <= $cols; $j++ ){
		$R->[$v][$j] = $M->[$i][$j];
	    }
	    $v++;
	}

	# I supplement indexes of non zero columns
	for( my $j = 0; $j <= $cols; $j++ ){
	    $zc{$j} += 0;
	}

	$v = 0;
	# replacement of columns
	for my $j ( sort { $zc{$b} <=> $zc{$a} } keys %zc ){
	    for( my $i = 0; $i <= $rows; $i++ ){
		$M->[$i][$v] = $R->[$i][$j];
	    }
	    $v++;
	}
    }

    undef %zr;
    undef %zc;

    for(my $n = 0; $n < $rows; $n++){
	last if $n > $cols; # other zero rows

	# Replacement of zero diagonal element
        my $s = 0;
M_Gauss1:
	while( $M->[$n][$n] == 0 && $s < $rows - $n ){
	    $s++;
	    for( my $j = $n + 1; $j <= $cols; $j++ ){

		if( $M->[$n][$j] ){ # no zero element
		    # shift of columns $n <-> $j
		    for( my $i = 0; $i <= $rows; $i++ ){
			($M->[$i][$j], $M->[$i][$n]) = ($M->[$i][$n], $M->[$i][$j]);
		    }
		    last M_Gauss1;
		}
	    }

	    # all zero elements in row
	    for( my $j = $n; $j <= $cols; $j++ ){
		# shift of rows $n <-> $n+1
		($M->[$n][$j], $M->[$n+$s][$j]) = ($M->[$n+$s][$j], $M->[$n][$j]);
	    }
	}

	last unless $M->[$n][$n]; # zero elements of rows

	for( my $i = $n+1; $i <= $rows; $i++ ){

	    # Divisibility check totally
	    my($k,$d);
	    my $b = $M->[$i][$n] / $M->[$n][$n];
	    if($b == int($b)){
		$k = -$b;
		$d = 1;
	    }else{
		$k = -$M->[$i][$n];
		$d = $M->[$n][$n];
	    }

	    # column (element in row)
	    for( my $j = $n; $j <= $cols; $j++ ){
		$M->[$i][$j] = $M->[$i][$j]*$d + $k*$M->[$n][$j];
	    }
	}

	my $gcd = bgcd( @{$M->[$n]}[$n..$cols] );
	if($gcd > 1){
	    $_ = $_ / $gcd for @{$M->[$n]}[$n..$cols];
	}
    }
    $M;
}


# Interger determinant for quadratic matrix
# input:
#	$Dimention
#	facultative parameter
# return:
#	determinant
#	undef
sub Det{
    my( $M_ini, $opt ) = @_;

    my $dm = $#{ $M_ini }; # dimension of matrix

    return $M_ini->[0][0] if $dm < 1; # dim matrix = 1

    # Check Float (Fraction)
    my $fraction = 0;
    unless( exists $opt->{'int'} && $opt->{'int'} ){
	for my $i ( @$M_ini ){
	    $fraction = &_max_len_frac($i, $fraction);
	}
    }

    my $M = $fraction ? [ map{ [ map{ $_ * 10**$fraction } @$_ ] } @$M_ini ] :
			[ map{ [ @$_ ] } @$M_ini ]; # copy

    my $exch = 0; # number of permutations
    my $denom = 1;

    for( my $i = 0; $i < $dm; $i++ ){ # take all the matrix rows except 1st
	my $minV = abs( $M->[$i][$i] );

	if( ! $minV && $i < $dm - 1 ){
	    my $minN = $i;

	    # Search of row with abs minimal element on the diagonal
	    for( my $j = $i + 1; $j <= $dm; $j++ ){
		my $v = abs( $M->[$j][$i] );

		if( $v && ($v < $minV || ! $minV) ){
		    $minN = $j;
		    $minV = $v;
		}
	    }

	    return 0 unless $minV; # determinant = 0

	    ( $M->[$i], $M->[$minN] ) = ( $M->[$minN], $M->[$i] );
	    $exch++;
	}

	my $v1 = $M->[$i][$i];

	for( my $j = $i + 1; $j <= $dm; $j++ ){
	    my $v2 = $M->[$j][$i];

	    for( my $k = $i + 1; $k <= $dm; $k++ ){
		$M->[$j][$k] = ( $M->[$j][$k] * $v1 - $M->[$i][$k] * $v2 ) / $denom;
	    }
	}
	$denom = $v1;
    }

    for( $M->[$dm][$dm] ){
	if( $_ < 0 ){
	    $_ = abs;
	    $exch++;
	}
	$_ = 1 + int if abs($_ - int ) >= 0.5;	# Rounding
	$_ *= -1 if $exch%2;

	return $fraction ? $_ / 10**($fraction * ($dm + 1)) : $_;
    }
}


sub Solve_Det{
    my $M = shift || croak("Missing matrix");
    my $B = shift || croak("Missing vector");
    my $opts = shift;

    my $rows = $#{ $M };	# number of rows
    my $cols = $#{ $M->[0] };	# number of columns

    my $t = &test_matrix( $M );
    if( $t > 3 ){
	croak("Use of uninitialized value in matrix");

    }elsif( $t ){
	croak("Matrix is not quadratic");
    }

    croak("Vector doesn't correspond to a matrix") if $rows != $#{$B};
    croak("Use of uninitialized value in vector") if scalar( grep ! defined $_, @$B );

    if( defined $opts ){
	if( exists $opts->{'eqs'} ){
	    die "Unknown parameter \'$opts->{'eqs'}\'!\n"
		unless $opts->{'eqs'}=~/^(?:row|column)/i;
	}else{
	    $opts->{'eqs'} = 'row';
	}
    }else{
	$opts->{'eqs'} = 'row';
    }

    my $solve;

    # main determinant
    my $det_main = &Det( $M, $opts ) || return undef; # no one solution

    for( my $v = 0; $v <= $cols; $v++ ){

	my $R; # copy of matrix M
	for( my $i = 0; $i <= $rows; $i++ ){

	    if($opts->{'eqs'}=~/^col/i){
		$R->[$i] = $v == $i ? $B : $M->[$i];

	    }else{
		for( my $j = 0; $j <= $cols; $j++ ){
		    $R->[$i][$j] = $v == $j ? $B->[$i] : $M->[$i][$j];
		}
	    }

	}

	my $det = &Det( $R, $opts );
	my $dm = $det_main;

	if( $det ){
	    if( $dm < 0 ){
		$dm *= -1;
		$det *= -1;
	    }

	    # Check Float (Fraction)
	    my $max_frac = &_max_len_frac([$dm, $det], 0);

	    if($max_frac){
		$_ *= 10**$max_frac for($dm, $det);
	    }

	    my $gcd = bgcd(abs($det), abs($dm));
	    if($gcd > 1){
		$det = $det / $gcd;
		$dm = $dm / $gcd;
	    }

	    $solve->[$v] = $dm == 1 ? $det : "$det/$dm";

	}else{
	    $solve->[$v] = 0;
	}
    }
    $solve;
}


sub test_matrix{
    my $M = shift;

    return 4 if scalar( grep{ grep ! defined $_, @{$_} } @$M );

    my $res = scalar( grep $#{$_} != $#{ $M }, @$M ) ? 1 : 0;	# quadra
    $res += 2 if scalar( grep $#{$_} != $#{ $M->[0] }, @$M );	# reqtan;

    $res;
}


sub _max_len_frac{
    my($M, $max_frac) = @_;

    for( @$M ){
	next if Math::BigInt::is_int($_);
	my(undef, $frac) = Math::BigFloat->new($_)->length();
	$max_frac = $frac if $frac > $max_frac;
    }
    $max_frac;
}

1;
__END__