#!/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_multiplications )."\n";
print "NM=$num_multiplications\n";