Hi, here is a modified version of your code to show the number of matrix multiplications too:
#!/usr/bin/perl
# [id://1229234]
use warnings;
use strict;
use PDL;
sub matrix_power {
my ($M,$n,$nm) = @_;
if($n<0) {
return matrix_power(1/$M, -$n, $nm);
} elsif ($n==0) {
return 1;
} elsif ($n==1) {
return $M;
} elsif (0 == $n % 2) {
$$nm = $$nm + 1;
print "(AxA)^".int($n/2)."\n";
return matrix_power($M x $M, $n/2, $nm);
} else {
$$nm = $$nm + 2;
print "Ax((AxA)^".int(($n-1)/2).")\n";
return $M x matrix_power($M x $M, ($n-1)/2, $nm);
}
}
my $x = pdl([0,1], [2,3]);
my $num_multiplications = 0;
my $power = 100;
print "M ** $power = " . matrix_power( $x, $power, \$num_multiplicatio
+ns )."\n";
print "NM=$num_multiplications\n";
bw, bliako |