#!/usr/bin/env perl

# author: bliako
# date: 02/02/2019
# for:  https://perlmonks.org/?node_id=1229234
#
# Usage: mm.pl power [number-of-trials]
# Matrix is hard-coded.
#
# It calculates a matrix m to power p as efficiently as the author man
+aged.
# The brute-force solution is : m x m x ... x m for p times.
# But, there are gains to harvest
# if we re-use results from intermediate multiplications.
# For example, m^14 = m^12 x m^2
# So, first calculate m^2 = m x m and cache it (that's 1 multiplicatio
+n).
# Then calculate m^4 = m^2 x m^2 and cache it (that's another one).
# Then calculate m^6 = m^4 x m^2 and cache it (thats' another one).
# Then calculate m^8 = m^4 x m^4 and cache it (thats' another one).
# Then calculate m^12 = m^8 x m^4 and cache it (thats' another one).
# Then calculate m^14 = m^12 x m^2 and cache it (thats' another one).
# ...
# Then multiply the cached results with 1 more multiplications.
# Done with 7 multiplications!
# (as opposed to p-1=13 required by the brute-force approach)
# Now, there is scope for even more optimisation using
# the decomposition m^14 = m^4 x m^4 x m^4 x m^2
# which calculates the m^14 using only 5 multiplications!
# So, there is a naive randomised algorithm which picks plans
# at random and choses the one with the fewer number of multiplication
+s
# (without doing the actual multiplications)
# That's the "plan" which is then returned as a code-ref which is
# executed to yield the matrix to the asked power using PDL's matrix m
+ult (x)
#
# NOTE: m^a x m^b = m^(a+b)
#
use strict;
use warnings;

use PDL;
use PDL::Matrix;
use List::Util qw/shuffle/;

my $p = $ARGV[0]; # power to raise
my $n = $ARGV[1] // 500; # number of trials
# the matrix
my $m = PDL::Matrix->pdl([[3,5],[2,4]]);
if( $p > 65 ){
    print STDERR "$0 : power is too large, so I cowardly refuse to do 
+any matrix multiplications, only the plan will be reported.\n";
    make_plan($p, $n);
} else {
    my $res = make_plan($p, $n)->($m);
    print "$0 : ". $m . " to the power of $p is".$res;
}
print "$0 : done.\n";

# Given a power (P), it tries different sets
# of even numbers whose sum equals P:
# e.g. P=14 yields 2,2,2,2,2,2,2 or 4,2,6,2 or ...
# Odd-P has just a 1 at the end.
# Then it tries each decomposition for how many multiplications would
# require and returns a code-ref to execute the plan given an input PD
+L
# matrix.
# No matrix multiplications below, just trying different decomposition
+s
# with a naive randomised algorithm (it shuffles the items and picks t
+he 1st 2)
# There is definetely scope for improvement, possibly with a genetic a
+lgorithm?
sub    make_plan {
    my $P = $_[0];
    my $trials = $_[1] // 500;

    my $best_plan = undef;
    my $bestTstr = undef;
    my $min_num_multiplications = 1000000;
    my ($sum, $nT);
    while( --$trials > 0 ){
        my @T = (2) x int($P/2);
        my $Tstr = join(",", @T);
        my $num_multiplications = 1; # the m^2 case

        my %cache = ();
        my $aplan = '$cache{2} = $m x $m;'."\n";
        while(1){
            @T = split(/,/, $Tstr);
            if( ($nT=scalar(@T)) == 1 ) { last }
            my @shuffledT = List::Util::shuffle(@T);
            my ($T0, $T1) = ($shuffledT[0], $shuffledT[1]);
            $sum = $T0 + $T1;
            $Tstr =~ s/${T0},${T1}/$sum/g;
            $Tstr =~ s/,+/,/g; $Tstr =~ s/^,//; $Tstr =~ s/,$//;
            if( ! exists $cache{$sum} ){
                $aplan .= '$cache{'.$sum.'} = $cache{'.$T1.'} x $cache
+{'.$T0.'};'."\n";
                $cache{$sum} = 1;
                $num_multiplications++;
            }
        }
        $num_multiplications += scalar(@T)-1;
        if( $num_multiplications < $min_num_multiplications ){
            $best_plan = $aplan; $bestTstr = $Tstr;
            print "make_plan() : found better plan with $num_multiplic
+ations multiplications: $Tstr\n";
            $min_num_multiplications = $num_multiplications;
        }
        print "make_plan() : trial $trials found $num_multiplications 
+multiplications, best is $min_num_multiplications ...\n" if $trials %
+ 10 == 0;
    }
    $best_plan .= 'my $result = ';
    $best_plan .= '$cache{'.$_.'} x ' foreach split /,/, $bestTstr;
    if( $P % 2 ){
        $min_num_multiplications++;
        $best_plan .= '$m;'."\n";
    } else { $best_plan =~ s/ x $/;\n/ }
    my @bestT = split /,/, $bestTstr;
    print "\nmake_plan() : final plan with $min_num_multiplications mu
+ltiplications.\n";
    print "make_plan() : worst case requires ".($P%2?(int($P/2)+1):int
+($P/2))." multiplications and the naive case ".($P-1).".\n";
    print "make_plan() : the plan:\n$best_plan\n";

    $best_plan = 'my $m=$_[0]; my %cache=(); '.$best_plan.'return $res
+ult';
    my $ret = eval("sub{ $best_plan }");
    die "ERROR in eval of\nbegin code>>>>\n$best_plan\n<<<end code\n\n
+" if $@;
    return $ret
}
1;
[download]