95 lines
2.6 KiB
Perl
95 lines
2.6 KiB
Perl
use strict;
|
|
package SquareMatrix;
|
|
use Carp; # standard, "it's not my fault" module
|
|
|
|
use overload (
|
|
'""' => \&_string, # overload string operator so we can just print
|
|
'*' => \&_mult, # multiplication, needed for expo
|
|
'*=' => \&_mult, # ditto, explicitly defined to trigger copy
|
|
'**' => \&_expo, # overload exponentiation
|
|
'=' => \&_copy, # copy operator
|
|
);
|
|
|
|
sub make {
|
|
my $cls = shift;
|
|
my $n = @_;
|
|
for (@_) {
|
|
# verify each row given is the right length
|
|
confess "Bad data @$_: matrix must be square "
|
|
if @$_ != $n;
|
|
}
|
|
|
|
bless [ map [@$_], @_ ] # important: actually copy all the rows
|
|
}
|
|
|
|
sub identity {
|
|
my $self = shift;
|
|
my $n = @$self - 1;
|
|
my @rows = map [ (0) x $_, 1, (0) x ($n - $_) ], 0 .. $n;
|
|
bless \@rows
|
|
}
|
|
|
|
sub zero {
|
|
my $self = shift;
|
|
my $n = @$self;
|
|
bless [ map [ (0) x $n ], 1 .. $n ]
|
|
}
|
|
|
|
sub _string {
|
|
"[ ".join("\n " =>
|
|
map join(" " => map(sprintf("%12.6g", $_), @$_)), @{+shift}
|
|
)." ]\n";
|
|
}
|
|
|
|
sub _mult {
|
|
my ($a, $b) = @_;
|
|
my $x = $a->zero;
|
|
my @idx = (0 .. $#$x);
|
|
for my $j (@idx) {
|
|
my @col = map($a->[$_][$j], @idx);
|
|
for my $i (@idx) {
|
|
my $row = $b->[$i];
|
|
$x->[$i][$j] += $row->[$_] * $col[$_] for @idx;
|
|
}
|
|
}
|
|
$x
|
|
}
|
|
|
|
sub _expo {
|
|
my ($self, $n) = @_;
|
|
confess "matrix **: must be non-negative integer power"
|
|
unless $n >= 0 && $n == int($n);
|
|
|
|
my ($tmp, $out) = ($self, $self->identity);
|
|
do {
|
|
$out *= $tmp if $n & 1;
|
|
$tmp *= $tmp;
|
|
} while $n >>= 1;
|
|
|
|
$out
|
|
}
|
|
|
|
sub _copy { bless [ map [ @$_ ], @{+shift} ] }
|
|
|
|
# now use our matrix class
|
|
package main;
|
|
|
|
my $m = SquareMatrix->make(
|
|
[1, 2, 0],
|
|
[0, 3, 1],
|
|
[1, 0, 0] );
|
|
print "### Order $_\n", $m ** $_ for 0 .. 10;
|
|
|
|
$m = SquareMatrix->make(
|
|
[ 1.0001, 0, 0, 1 ],
|
|
[ 0, 1.001, 0, 0 ],
|
|
[ 0, 0, 1, 0.99998 ],
|
|
[ 1e-8, 0, 0, 1.0002 ]);
|
|
|
|
print "\n### Matrix is now\n", $m;
|
|
print "\n### Big power:\n", $m ** 100_000;
|
|
print "\n### Too big:\n", $m ** 1_000_000;
|
|
print "\n### WAY too big:\n", $m ** 1_000_000_000_000;
|
|
print "\n### But identity matrix can handle that\n",
|
|
$m->identity ** 1_000_000_000_000;
|